cas-bioinf/covid19retrospective
What the Package Does (One Line, Title Case)

Package index
Search the cas-bioinf/covid19retrospective package
Vignettes
Functions
110
Source code
30
Man pages
3
Home
/
GitHub
/
cas-bioinf/covid19retrospective
/

In cas-bioinf/covid19retrospective: What the Package Does (One Line, Title Case) 

knitr::opts_chunk$set(echo = FALSE, fig.width=8,out.width = "100%", cache=TRUE)

devtools::load_all()
library(tidyverse)
library(patchwork)
library(tsiMisc) # devtools::install_github("tsieger/tsiMisc")
theme_set(cowplot::theme_cowplot())
data <- read_data_for_analysis()
wide <- data$marker_data_wide

n_patients <- nrow(data$patient_data)
n_hospitals <- length(unique(data$patient_data$hospital_id))

epidemiology_cz <- read_csv(here::here("public_data", "nakazeni-vyleceni-umrti-testy.csv"), col_types = cols(
  datum = col_date(format = ""),
  kumulativni_pocet_nakazenych = col_double(),
  kumulativni_pocet_vylecenych = col_double(),
  kumulativni_pocet_umrti = col_double(),
  kumulativni_pocet_testu = col_double(),
  kumulativni_pocet_ag_testu = col_double()
))

epidemiology_cz_october <- epidemiology_cz %>% filter(datum == "2020-10-01")

# hospitalized_cz <- read_csv(here::here("public_data", "hospitalizace.csv"), col_types = cols(
#   datum = col_date(format = ""),
#   pacient_prvni_zaznam = col_double(),
#   kum_pacient_prvni_zaznam = col_double(),
#   pocet_hosp = col_double(),
#   stav_bez_priznaku = col_double(),
#   stav_lehky = col_double(),
#   stav_stredni = col_double(),
#   stav_tezky = col_double(),
#   jip = col_double(),
#   kyslik = col_double(),
#   hfno = col_double(),
#   upv = col_double(),
#   ecmo = col_double(),
#   tezky_upv_ecmo = col_double(),
#   umrti = col_double(),
#   kum_umrti = col_double()
# ))
# 
# hospitalized_cz_october <- hospitalized_cz %>% filter(datum == "2020-10-01")
# hospitalized_cz_october$kum_pacient_prvni_zaznam

# epidemiology_cz_october$kumulativni_pocet_nakazenych
# epidemiology_cz_october$kumulativni_pocet_umrti

Introduction

The Covid-19 pandemic caused by severe acute respiratory syndrome coronavirus (SARS-CoV-2) has, as of June 2021, led to over 172 million cases and over 3.7 million deaths. The present study was designed and conducted during March - October 2020, when Czech Republic experienced a relatively mild first wave of the pandemic due to early and strict lockdowns. Low numbers of cases continued throughout the summer but during September and October, after most of the data collection for this study concluded, the number of cases was raising again. On October 1st 2020, Czech Republic had accumulated r as.character(epidemiology_cz_october$kumulativni_pocet_nakazenych) total confirmed Covid-19 cases and r epidemiology_cz_october$kumulativni_pocet_umrti confirmed Covid-19 related deaths [@http://zotero.org/users/5567156/items/SWNTWFLU].

At the time the study was conducted, the proposed treatments included antivirals approved for other indications (chloroquine, hydroxychloroquine, lopinavir/ritonavir, remdesivir, favipiravir, umifenovir), azithromycin, corticosteroids, immunoglobulins, tocilizumab and convalescent plasma [@http://zotero.org/users/5567156/items/8XWTUMUI; @http://zotero.org/users/5567156/items/6Q8Z8VJR]. Notably, the anti-malarial and anti-rheumatic drug hydroxychloroquine and the macrolide antibiotic azithromycin showed promise in early data and were broadly available and thus were frequently used in the early stages of the pandemic. Remdesivir, previously designed and approved for Ebola, SARS and MERS, also reported good initial results. However, during spring and summer 2020 remdesivir was available in Czech Republic only in limited amounts via compassionate use programme. The RECOVERY trial reported positive results of coroticosteriod dexamethasone for severe cases in June 2020 [@http://zotero.org/users/5567156/items/W8KK3TX2], but at this point the number of Covid-19 patients hospitalized in Czech Republic was low and dexamethasone thus did not see wider use until later in the pandemic.

Our understanding of the efficacy of Covid-19 treatments has improved substantially since the present study was conducted. As of April 2021, the pharmacological treatments that were deemed to be beneficial for at least one outcome in a systematic review of randomized trials were the corticosteroid dexamethasone (mortality, mechanical ventilation), colchicine (mortality, length of hospital stay), the antiviral remdesivir (mechanical ventilation), Janus kinase inhibitors (mechanical ventilation, duration of ventilation), IL-6 inhibitors (mechanical ventilation, length of hospital stay), the antiviral favipiravir (length of hospital stay, resolution of symptoms) and the anti-androgen proxalutimide (admission to hospital). Hydroxychloroquine, interferon beta, lopinavir-ritonavir, azithromycin, vitamin C, vitamin D, anticoagulants and ACE inhibitors were considered to not be better than standard of care and lopinvair-ritonavir showed evidence of harm, although most of the conclusions were considered to be of low certainty [@http://zotero.org/users/5567156/items/YP62P6XJ]. Interestingly, in observational studies, hydroxychloroquine was often found to be associated with better outcomes [e. g., @http://zotero.org/users/5567156/items/WZPNBRYN; @http://zotero.org/users/5567156/items/4HR7MUJ6; @http://zotero.org/users/5567156/items/NXC9CDI9]. No benefit was also observed in a meta-analysis of randomized trials of convalescent plasma treatment [@http://zotero.org/users/5567156/items/JCV6R3RC].

High IL-6, D-dimer values were observed to be associated with worse outcome and increased disease severity [@http://zotero.org/users/5567156/items/K3N3VG2X]. Large study of electronic health records [@http://zotero.org/users/5567156/items/ARTB8HYZ] showed an increase in C-reactive protein in early disease and increase of D-dimer and white blood cell count in later stages of the disease.

An ongoing challenge in evaluating Covid-19 treatments is that the analysis and interpretation of the data is often inappropriate or misleading, most notably interpreting lack of evidence due to small sample size as evidence of no effect [see e.g., @http://zotero.org/users/5567156/items/5PZEEPJV; @http://zotero.org/users/5567156/items/AZDZHLJV].

Additionally, many methods for predicting disease severity of Covid-19 were published, but the methods are at high risk of bias and lack external validation [@http://zotero.org/users/5567156/items/TS57HJ3T].

The present study aims to describe the outcomes and disease course of hospitalized patients with mild to severe clinical presentation in a multicentric Czech cohort during the early stages of the pandemic, explore the association between the outcomes and pharmacological interventions and to provide external validation to previously published prognostic models for Covid-19 severity.

Methods

Patients and study design

ethics_data <- readxl::read_excel(here::here("manuscript", "ethics.xlsx"), sheet = "Sheet1", na = "NA")
ethics_string <- ethics_data %>% pull(Site) %>% paste0(collapse = ", ")

n_ecmo <- data$marker_data_wide %>% filter(breathing == "ECMO") %>% pull(patient_id) %>% unique() %>% length()
n_nippv <- data$marker_data_wide %>% filter(breathing == "NIPPV") %>% pull(patient_id) %>% unique() %>% length()

n_second_wave <- data$patient_data %>% filter(!first_wave) %>% nrow()

A convenience sample of patients from r n_hospitals sites was collected. The study sites span the whole spectrum of sizes from large university hospitals in major cities with multiple dedicated Covid-19 wards (Thomayer University Hospital in Prague, Motol University Hospital in Prague, Kralovské Vinohrady University Hospital in Prague, General University Hospital in Prague, University Hospital in Pilsen) through major regional/specialized hospitals (Na Homolce Hospital in Prague, Military Hospital Olomouc) as well as smaller hospitals caring for just several Covid-19 patients at a time (AGEL Hospital Nový Jičín, Hořovice Hospital, Třebíč Hospital). The sites were chosen based on availability and willingness of the personnel to participate in data collection. None of the study sites was exclusively dedicated to treating Covid-19 patients. For each site, the dataset contains all patients hospitalized in the participating wards over the data collection period. The data collection started at the onset of the Covid-19 pandemic in March 2020 (except for one site where some older records were inaccessible), but the end date for collection differs between sites due to time constraints of the participating physicians. Three sites included total of r n_second_wave patients that could be considered part of "second wave" (admitted after September 1st). Last patient included in the dataset was admitted on October 12th. See S1 Fig for per-site data collection periods. Patients over the age of 18 were included if they had PCR-confirmed infection of SARS-CoV-2 and were not participating in a clinical trial of any Covid-19 pharmacotherapy.

Not all patients developed pneumonia or other symptoms of Covid-19. All patients received the standard of care which could include supplemental oxygen and ventilation and antibiotics for bacterial superinfections, as determined by the attending physician. Some patients were not indicated for all treatment modalities (especially mechanical ventilation) based on decision of the attending physician and underlying patient condition. We note that the participating sites were not homogeneous in either patient population or treatment protocols. The choice of pharmacological treatment was based on the decision of the attending clinician and its availability.

The study was approved by the Ethical committees of r ethics_string, all data were collected in fully anonymized form. Data was collected between June and October 2020 for patients that were treated between March and October 2020.

Data collection

We collected data on comorbidities and information about disease progression on daily resolution including breathing support required, oxygen flow rate, experimental anti-Covid-19 and antimicrobial drugs taken and several laboratory markers (PCR positivity for SARS-CoV-2, C-reactive protein, D-dimer, Interleukin 6, Ferritin, lymphocyte count). Full protocol for data collection is attached in S1 file and the data collection table in S2 file. Due to very low number of patients using extra-corporeal membrane oxygenation (N = r n_ecmo) or non-invasive positive pressure ventilation (N = r n_nippv) in our sample, we merged those categories with mechanical ventilation.

Statistical analysis

The character of the convenience sample does not allow for a proper assessment of the association between treatments and patient outcomes, because the treatments had not been assigned to patients at random but were only observed retrospectively. This can be partially remedied by adjusting for patient characteristics in the analysis, but such adjustments will always be imperfect and the analysis needs to be treated as exploratory and interpreted cautiously.

Since many details of analysis may influence the conclusions made, we performed multiverse analysis [@http://zotero.org/users/5567156/items/DXWM259V] and report results for all the hypothesis tested across multiple different models using both frequentist and Bayesian paradigms. For each model class we worked with several possible sets of adjustments. All analyses were performed in the R language [@http://zotero.org/users/5567156/items/67XKBC2W], visualization and data cleaning was run with the tidyverse package [@http://zotero.org/users/5567156/items/8EWY8TKM]. 

 # and - where applicable - with multiple estimands of the outcome as different tasks require different estimands [@http://zotero.org/users/5567156/items/MH797GNC]

First class of models are frequentist survival and multistate models under the proportional hazards assumption as implemented in the coxph function from the survival package [@http://zotero.org/users/5567156/items/I9UVTUJF]. We primarily use a model with competing risks for death and discharge from hospital (see Fig \@ref(fig:modelstates)a).

Second class of models are Bayesian hidden Markov models (HMM) of disease progression implemented via a custom extension to the brms package [@http://zotero.org/users/5567156/items/7ER3ZSER]. The parametrization of the HMM is inspired by Williams et al.[@http://zotero.org/users/5567156/items/F9K2L5Q8]: the actual process of disease is assumed to be continuous and allow only for transitions between neighboring states (as shown in Fig \@ref(fig:modelstates)b, c). The total probability of transition between any two states over the period of a day is then computed as the total probability of transition across all possible paths. This class of models does not satisfy the proportional hazards assumption, instead, it is assumed the process has the Markov property - i.e., that the (potentially unobserved) state and the covariates at a given day contain all the information to determine probabilities of the states on the next day. We use two versions of such models, one working solely with the observed breathing support and one assuming a hidden improving/worsening distinction. All of the hidden Markov models take into account whether best supportive care was initiated and a patient was thus not indicated to progress to more intensive treatment modalities.

knitr::include_graphics(here::here("manuscript","static_figs","model_states.png"), auto_pdf = TRUE)

Finally, we used a set of Bayesian regression models implemented with the brms package [@http://zotero.org/users/5567156/items/7ER3ZSER]. Those included overall survival, state at day 7 or 28 as either binary or categorical outcome and a Bayesian version of the Cox proportional-hazards model.

#Since `stan_jm` does not support competing risks, we imputed the second event...

Except for age, sex and comorbidities, all covariates are treated as time-varying, e.g., the effect of taking a drug is only included for the days after the drug was taken. More details on the exact model formulations can be found in the supplementary statistical analysis in S3 file.

#Sex differences are to be expected for both sociological and biological [@http://zotero.org/users/5567156/items/5KQZ3XEE] reasons.

Evaluating prognostic models

We searched the living systematic review of Covid-19 prognostic models [@http://zotero.org/users/5567156/items/TS57HJ3T] for those that could be applied to our dataset (i.e., where we have gathered all the input features). We primarily focused on the Area Under Receiver Operating Characteristic Curve (AUC), and its bootstrapped 95% confidence intervals which we computed using the pROC package [@http://zotero.org/users/5567156/items/UN6BIC67]. When there were multiple reasonable ways to evaluate the outcome or a predictor in our dataset, we computed and reported all of those options. We used two simple scores with age or the decade of age as the sole predictor to have a baseline to compare the scores against.

Complete code for all analyses is available at https://github.com/cas-bioinf/covid19retrospective/.

Results

Baseline characteristics

In total, we were able to gather data for r n_patients patients, see Table \@ref(tab:summarytable) for the overall characteristics of the patient sample and several subgroups we used in the analysis, including treatments taken. Counts of all treatment combinations are shown in S2 Fig and S3 Fig shows outcomes by study site, demonstrating quite large hospital-specific differences. The dataset includes 19 patients already reported in a study of inflammatory signatures of Covid-19 [@http://zotero.org/users/5567156/items/IA4FQCF3]. 

data_table_1 <- rbind(
  data$patient_data %>% mutate(group = "All"),
  data$patient_data %>% mutate(group = if_else(ever_hcq, "HCQ", "No HCQ")),
  data$patient_data %>% filter(ever_favipiravir) %>% mutate(group = "Favipiravir")
  # data$patient_data %>% filter(any_d_dimer) %>% mutate(group = "Measured D-dimer"),
  # data$patient_data %>% filter(any_IL_6) %>% mutate(group = "Measured IL-6")
  )


data_table_1 %>% group_by(group) %>% 
  summarise(
    N = n(),
    `Distinct sites` = length(unique(hospital_id)),
    Male = sum_percent(sex == "M"),
    `Age (mean, IQR)` = mean_iqr(age),
    `Admitted for Covid` = sum_percent(admitted_for_covid),
    `Took hydroxychloroquine` = sum_percent(ever_hcq),
    `Took azithromycin` = sum_percent(ever_az),
    `Took dexamethasone` = sum_percent(ever_dexamethasone),
    `Took favipiravir` = sum_percent(ever_favipiravir),
    `Took remdesivir` = sum_percent(ever_remdesivir),
    `Convalescent plasma` = sum_percent(ever_convalescent_plasma),
    `Ischemic Heart Disease` = sum_percent(ischemic_heart_disease),
    `Takes antihypertensives` = sum_percent(has_hypertension_drugs),    
    `Heart Failure` = sum_percent(heart_failure),
    COPD = sum_percent(COPD),
    Asthma = sum_percent(asthma),
    `Other lung disease` = sum_percent(other_lung_disease),
    Diabetes = sum_percent(diabetes),
    `Renal Disease` = sum_percent(renal_disease),
    `Liver Disease` = sum_percent(liver_disease),
    `Smoking` = sum_percent(smoking),
    `BMI (mean, IQR)` = mean_iqr(BMI),
    `Best supportive care` = sum_percent(!is.na(best_supportive_care_from)),
    Death = sum_percent(outcome == "Death"),
    Discharged = sum_percent(outcome == "Discharged"),
    .groups = "drop"
  ) %>%
  mutate(group = factor(group, levels= c("All", "HCQ", "No HCQ", "Favipiravir"), ordered = TRUE)) %>%
  arrange(group) %>%
  as.data.frame() %>%
  column_to_rownames("group") %>%
  t() %>% 
  knitr::kable(caption = "Patient characteristics for the overall sample and treatment subgroups. Note that the favipiravir subgroup is not exclusive with either the HCQ or No HCQ group. IQR = interquartile range, COPD = Chronic obstructive pulmonary disease, BMI = body-mass index, Best supportive care = patient was not indicated to undergo more intensive treatment modality.", booktabs = TRUE)

# Note that the 'Measured D-dimer' and 'Measured IL-6' are neither exclusive nor exhaustive.

median_ratios <-  data$marker_data %>% 
  filter(marker %in% c("lymphocyte_count", "d_dimer", "CRP", "IL_6")) %>%
  group_by(patient_id, marker) %>%
  filter(n() > 2) %>%
  summarise(ratio = max(value) / min(value), .groups = "drop") %>%
  group_by(marker) %>%
  summarise(med_ratio = median(ratio), .groups = "drop") %>%
  column_to_rownames("marker") %>% as.matrix() %>% round(1)

Disease progression

In Fig \@ref(fig:progression) we show the overall disease progression for all patients and in Fig \@ref(fig:selectedmarkers) we show the time-course of a subset of the markers we have measured. The data show some interesting patterns: patients with low Interleukin-6 or D-dimer values are overrepresented among patients with better outcomes, most patients had high CRP upon admission and for many the CRP levels stayed elevated over the whole hospitalization. However, the limited nature of the data does not allow for any statistically robust conclusions. We also see that the marker levels were not substantially stratified by study site. Those patterns should however be interpreted with care due to systematic missingness of the data - in particular, patients that fared worse were probably more likely to have the markers measured. However, we believe this kind of patient-level view is useful to appreciate the extent of both between-patient and within-patient variability.

The between-patient variability is notable even across outcomes - when ordering the patients by the highest CRP levels experienced throughout the hospital stay, the top 20% of patients that breathed ambient air for the whole hospitalization experienced higher levels than the bottom 20% of patients that required ventilation or died. This overlap is even larger when comparing only against the patients that died and D-dimer, Interleukin-6 and lymphocyte count also show a notably larger overlap than CRP (S4 Fig). 

my_paste <- function(x) {
  paste0(as.integer(x), collapse = "")
}



data_for_fig_progression <- wide  %>% 
  filter(!is.na(breathing_s), day >= 0) %>%
  arrange(hospital_id, patient_id, day) %>%
  mutate(patient_id = factor(patient_id), patient_id = fct_reorder(patient_id, breathing_s,  .fun = my_paste, .desc = FALSE)) %>%
  mutate(first_half = as.integer(patient_id) < length(levels(patient_id)) / 2 + 1, 
         is_outcome =  breathing_s %in% c("Discharged", "Death"),
         breathing_s = fct_recode(breathing_s, `Ambient air` = "AA")
         ) 

fig_progression <- data_for_fig_progression %>% 
  mutate(patient_id = fct_reorder(patient_id, breathing_s,  .fun = my_paste, .desc = TRUE)) %>%
  filter(!is_outcome) %>%
  ggplot(aes(x = patient_id, y = day, fill = breathing_s)) + 
  geom_tile(width = 0.7) + 
  geom_point(data = data_for_fig_progression %>% filter(is_outcome) %>% mutate(outcome = factor(breathing_s, ordered = FALSE)),
             inherit.aes = FALSE, aes(x = patient_id, y = day, shape = outcome), size = 0.9) +
  scale_fill_viridis_d("", direction = -1) +
  scale_y_continuous("Days since hospitalization") +
  # scale_color_manual("", values = c("black", "darkred"))+
  scale_shape_manual("", values = c(1, 15), guide = guide_legend(override.aes = list(size = 4))) +
  theme(axis.text.x = element_blank(), strip.background = element_blank(), strip.text = element_blank(),
        legend.direction = "horizontal", legend.position = "bottom") +
  theme(axis.title.x = element_blank(), axis.ticks.x = element_blank(), axis.line.x = element_blank())

fig_progression
ggsave(fig_progression, filename = "supp_figs_production/fig2.tiff", width = 7.5)

#data$marker_data %>%
# TODO: use SpO2 in AA as SpO2 native (and vice versa)
plot_marker_data <- function(markers, log_transform, smooth = FALSE, max_patients = NULL, max_day = NULL) {

  if(length(markers) > 1) {
    facet = facet_grid(marker ~ worst_breathing_s, scales = "free_y")
    scale_y_name = "Value"
  } else {
    facet = facet_wrap( ~ worst_breathing_s, nrow = 1)
    unit <- data$marker_data %>% filter(marker == names(markers)) %>% pull(unit) %>% unique()
    scale_y_name <- markers
  }

  data_for_plot <- data$marker_data %>% 
    filter(marker %in% names(markers), day >= 0) %>% 
    group_by(marker, patient_id) %>%
    filter(n() > 1) %>%
    ungroup() %>%
    mutate(marker = markers[marker]) 

  if(!is.null(max_patients)) {
    patient_markers_to_use <- data_for_plot %>% 
      select(marker, patient_id) %>%
      distinct() %>%
      group_by(marker) %>%
      slice_sample(n = max_patients) 

    data_for_plot <- data_for_plot %>%
      semi_join(patient_markers_to_use, by = c("patient_id", "marker"))
  }

  if(log_transform) {
    scale_y <- scale_y_log10(scale_y_name)
    data_for_plot <- data_for_plot %>% 
      group_by(marker) %>%
      mutate(value = if_else(value == 0, min(value[value > 0]) / 2, value)) %>%
      ungroup()
  } else {
    scale_y <- scale_y_continuous(scale_y_name)
  }

  if(smooth) {
    main_geom <- geom_smooth(aes(group = outcome))
  } else {
    main_geom <- geom_line(alpha = 0.5)
  }

  if(!is.null(max_day)) {
    expand <- expand_limits(x = max_day)
  } else {
    expand <- NULL
  }

  data_for_plot %>%
    inner_join(data$patient_data, by = c("hospital_id", "patient_id")) %>%
    mutate(outcome = fct_relevel(fct_collapse(outcome, Hospitalized = c("Hospitalized", "Transferred")),
            "Discharged", "Hospitalized", "Death"),
           worst_breathing_s =fct_recode(worst_breathing_s, `Ambient Air` = "AA")) %>%
    ggplot(aes(x = day, y = value,  group = patient_id, color = hospital_id)) + 
    main_geom +
    facet + 
    expand +
    #scale_color_brewer(type = "qual", palette = 2, guide = guide_legend(override.aes = list(alpha = 1, size = 2))) + 
    scale_y +
    scale_x_continuous("Days since hospitalization") +
    scale_color_discrete(guide = FALSE)
    #scale_color_viridis_d(guide = guide_legend(override.aes = list(alpha = 1, size = 2))) + scale_y
}

max_day <- data$marker_data_wide %>% filter(!is.na(CRP)) %>% pull(day) %>% max()
hide_axis_theme <- theme(axis.text.x = element_blank(),
                         axis.line.x = element_blank(),
                         axis.title.x = element_blank(),
                         axis.ticks.x = element_blank())
hide_strip_theme <- theme(strip.background = element_blank(), strip.text = element_blank())
title_theme <- theme(axis.title.y = element_text(size = 10))

fig_markers <- ((plot_marker_data(c("CRP" = "CRP"), log_transform = TRUE, max_day = max_day) + hide_axis_theme + title_theme) /
(plot_marker_data(c("d_dimer" = "D-dimer"), log_transform = TRUE, max_day = max_day) + hide_axis_theme + title_theme + hide_strip_theme) /
(plot_marker_data(c("lymphocyte_count" = "Ly"), log_transform = TRUE, max_day = max_day) + hide_axis_theme + title_theme + hide_strip_theme) /
(plot_marker_data(c("IL_6" = "IL-6"), log_transform = TRUE, max_day = max_day) + title_theme + hide_strip_theme)) 

fig_markers

ggsave(fig_markers, filename = "supp_figs_production/fig3.tiff", width = 7.5)
  #plot_marker_data(c("d_dimer" = "D-dimer", "IL_6" = "IL 6", "CRP" = "CRP"), log_transform = TRUE)
#plot_marker_data(c("d_dimer", "IL_6", "CRP"), log_transform = FALSE)

# plot_marker_data(c("CRP"), log_transform = TRUE)
# plot_marker_data(c("CRP"), log_transform = TRUE, max_patients = 30)
# plot_marker_data(c("CRP"), log_transform = TRUE, max_patients = 30)

Association between patients' characteristics and treatments

As noted above, the nature of the convenience sample did not enforce random assignment of treatments to patients. In fact, patients with worse baseline characteristics, which lead to worse outcomes, were less likely to receive hydroxychloroquine (see Fig \@ref(fig:hcqassoc)). This clearly creates a bias towards a positive effect of hydroxychloroquine on the outcome (and potentially for other treatments as well - most were used in combination with hydroxychloroquine), which, however, could be false. 

tmp<-data$patient_data[,c('age','sex','days_from_symptom_onset','admitted_for_covid','ischemic_heart_disease','has_hypertension_drugs','heart_failure',
  'COPD','asthma','other_lung_disease','diabetes','renal_disease','liver_disease','smoking','BMI','obesity','last_record','high_creatinin','high_pt_inr','low_albumin',
  'heart_problems','outcome','ever_hcq','ever_az','worst_breathing','had_invasive'),]
tmp$outcome[tmp$outcome=='Transferred']<-'Hospitalized'
tmp$outcome<-dropLevels(tmp$outcome)
tmp2<-tmp[tmp$outcome%in%c('Death','Discharged'),]
tmp2$outcome<-dropLevels(tmp2$outcome)
tmp2$outcome<-tmp2$outcome=='Death'

m.ihd<-glm(outcome~ischemic_heart_disease,tmp2,family='binomial')
m.hd<-glm(outcome~has_hypertension_drugs,tmp2,family='binomial')
m.hf<-glm(outcome~heart_failure,tmp2,family='binomial')
m.copd<-glm(outcome~COPD,tmp2,family='binomial')
m.asthma<-glm(outcome~asthma,tmp2,family='binomial')
m.old<-glm(outcome~other_lung_disease,tmp2,family='binomial')
m.dia<-glm(outcome~diabetes,tmp2,family='binomial')
m.rd<-glm(outcome~renal_disease,tmp2,family='binomial')
m.ld<-glm(outcome~liver_disease,tmp2,family='binomial')
m.s<-glm(outcome~smoking,tmp2,family='binomial')
m.ob<-glm(outcome~obesity,tmp2,family='binomial')
m.hc<-glm(outcome~high_creatinin,tmp2,family='binomial')
m.hi<-glm(outcome~high_pt_inr,tmp2,family='binomial')
m.la<-glm(outcome~low_albumin,tmp2,family='binomial')

m2.ihd<-glm(ever_hcq~ischemic_heart_disease,tmp2,family='binomial')
m2.hd<-glm(ever_hcq~has_hypertension_drugs,tmp2,family='binomial')
m2.hf<-glm(ever_hcq~heart_failure,tmp2,family='binomial')
m2.copd<-glm(ever_hcq~COPD,tmp2,family='binomial')
m2.asthma<-glm(ever_hcq~asthma,tmp2,family='binomial')
m2.old<-glm(ever_hcq~other_lung_disease,tmp2,family='binomial')
m2.dia<-glm(ever_hcq~diabetes,tmp2,family='binomial')
m2.rd<-glm(ever_hcq~renal_disease,tmp2,family='binomial')
m2.ld<-glm(ever_hcq~liver_disease,tmp2,family='binomial')
m2.s<-glm(ever_hcq~smoking,tmp2,family='binomial')
m2.ob<-glm(outcome~obesity,tmp2,family='binomial')
m2.hc<-glm(ever_hcq~high_creatinin,tmp2,family='binomial')
m2.hi<-glm(ever_hcq~high_pt_inr,tmp2,family='binomial')
m2.la<-glm(ever_hcq~low_albumin,tmp2,family='binomial')

models<-list(m.ihd,m.hd,m.hf,m.copd,m.asthma,m.old,m.dia,m.rd,m.ld,m.s,m.ob,m.hc,m.hi,m.la)
models2<-list(m2.ihd,m2.hd,m2.hf,m2.copd,m2.asthma,m2.old,m2.dia,m2.rd,m2.ld,m2.s,m2.ob,m2.hc,m2.hi,m2.la)
modelNames<-c('IHD','HD','HF','COPD','Asthma','LungD','Dia','RenalD','LiverD','Smoking','Obesity','HighCr','HighInr','LowAlb')
gtmp<-c()
for (i in rev(seq(along=models))) {
  ii<-length(models)-i+1
  m<-models[[i]]
  m2<-models2[[i]]
  ci<-confint(m)[2,]
  ci2<-confint(m2)[2,]
  gtmp<-rbind(gtmp,data.frame(est=coef(m)[2],low=ci[1],high=ci[2],name=modelNames[i],type='outcome'))
  gtmp<-rbind(gtmp,data.frame(est=-coef(m2)[2],low=-ci2[1],high=-ci2[2],name=modelNames[i],type='hcq'))
}
fig_assoc <- ggplot(gtmp,aes(y=name,x=est,xmin=low,xmax=high,color=type))+
  geom_vline(xintercept=0,color='blue')+
  geom_point(position=position_dodge(width=0.5))+
  geom_errorbarh(position=position_dodge(width=0.5),height= 0.5)+
  scale_x_continuous('Log odds ratio')+
  scale_y_discrete('')+
  scale_color_manual('Model of:',values=c('black','red'),labels=c('Death','No HCQ'))+
  theme(legend.background = element_rect(fill="gray95"),legend.position=c(.83,.7))

fig_assoc 
ggsave(fig_assoc, filename = "supp_figs_production/fig4.tiff", width = 7, height = 4.5)

# condition on the sum of comorbidities
pd<-data$patient_data
pdd<-within(pd,cnd<-ischemic_heart_disease+has_hypertension_drugs+heart_failure+COPD+    other_lung_disease+    renal_disease+  high_creatinin)
m.hcq_treatment2<-glm(ever_hcq~cnd,pdd,family='binomial')
m.hcq_treatment0<-glm(ever_hcq~1,pdd,subset=!is.na(cnd),family='binomial')

Taken quantitatively, the comorbidities known upon hospitalization were informative with respect to the future hydroxychloroquine treatment: the score representing the cumulative presence of ischemic heart disease, hypertension drugs, former heart failure, COPD, other lung diseases, renal disease, or high creatinine was associated with a lower chance of taking hydroxychloroquine over the course of the hospitalization (the chance was only r round(100*exp(coef(m.hcq_treatment2)[2]),1)%, 95% confidence interval (r round(100*exp(confint(m.hcq_treatment2)[2,]),1))%, Chi-square test in the logistic regression model, $\chi^2$=r round(abs(anova(m.hcq_treatment2,m.hcq_treatment0,test='Chisq'))[2,4],2), df=1, P=r round(anova(m.hcq_treatment2,m.hcq_treatment0,test='Chisq')[2,5],3)).

Association between treatments and outcomes

all_hypo_drugs <- read_all_multiverse_res()  %>%
  filter(model_check == "OK", !adj.first_wave, estimand_group == "")

first_fig_multiverse_drugs <- 6
last_fig_multiverse_drugs <- 8
supplementary_figure_markers <- last_fig_multiverse_drugs + 1 
supplementary_figures_multiverse <-  paste0(first_fig_multiverse_drugs, "-", last_fig_multiverse_drugs)

Here, we focus on hydroxychloroquine and azithromycin as those are the only treatments with larger sample size. We also investigate favipiravir as it is less well reported in the literature. Hydroxychloroquine was dosed almost exclusively in a 5-day regime starting with a loading dose of 800mg on the first day and followed by 400mg. Majority of patients complemented hydroxychloroquine with azithromycin while azithromycin was rarely used alone (see Table \@ref(tab:summarytable)). Azithromycin was most frequently dosed 250 or 500mg/day, but doses ranging from 100mg/day to 1500 mg/day were observed. Favipiravir was used only at one site with a loading dose of 3600mg on the first day, followed by at most 9 days with a 1600mg dose. All but one of the patients receiving favipiravir also received hydroxychloroquine. Treatment was initiated mostly within two days of admission (see S5 Fig).

The results of the multiverse analysis for association between hydroxychloroquine, azithromycin and favipiravir usage and death is shown in Fig \@ref(fig:multiverseresultsdrugs) - here, we only show models that were not found to have immediate problems representing the data well or computational issues (see S3 file for details). Results for all models we tested are reported in Sr supplementary_figures_multiverse Figs, with additional details in S3 file. The results do not change noticeably when only patients from the first wave are included (Sr supplementary_figures_multiverse Figs).

Most models report that using hydroxychloroquine is associated with lower risk of death. We must however bear in mind the potential bias noted in the previous section. Also, we see that for the HMM models, as we add adjustments the credible intervals do not widen but instead shift towards zero. This is a weak indication that further adjustments could drive the effect towards zero. We did not attempt to model additional adjustments as the models became computationally unstable. The case of hydroxychloroquine serves as a "control group" for our other results - since randomized trials give us high confidence that hydroxychloroquine does not substantially reduce mortality, we can be quite certain the associations we observe for hydroxychloroquine are just a measure of bias in the data. Additionally, our models either cannot determine the sign of association between azithromycin and risk of death or even show an increase in risk of death. This serves as a weak evidence that a substantial improvement in mortality from azithromycin is unlikely.

Most models exclude very strong association between increased risk of death and using favipiravir, but our data are necessarily quite limited, which is reflected in the very wide uncertainties around estimates. We also cannot put strict bounds on the association between favipiravir and length of hospitalization.

We also examined the association between treatments and length of hospital stay for all the patients that survived. Almost all models cannot discern the sign of the association for all treatments examined (Sr supplementary_figures_multiverse Figs). Similarly, we studied the association between D-dimer and Interleukin 6 and outcomes, with unconclusive results as well (Sr supplementary_figure_markers Fig) 

#In the current dataset, we see a possibly strong negative association between using hydroxychloroquine and risk of death across all models, although for several models we cannot rule out that the association is of negligible strength. Across all models, the data offer no definitive conclusions on the association between taking hydroxychloroquine and length of hospital stay (represented as risk for being discharged). After adjusting for hydroxychloroquine usage, we also cannot make any strong claims on the association of using azithromycin with either outcome. In all cases, the results seem to not be strongly affected by model choice.


my_lims <- c(-5.5,3.5)
multiverse_hcq1 <- plot_multiverse(all_hypo_drugs %>% filter(hypothesis == "hcq_death"), x_limits = my_lims, adjustments_to_hide = "first_wave")
multiverse_az1 <- plot_multiverse(all_hypo_drugs %>% filter(hypothesis == "az_death"), x_limits = my_lims, adjustments_to_hide = "first_wave")

# multiverse_fpv1 <- plot_multiverse(all_hypo_drugs %>% filter(hypothesis == "favipiravir_death"))
# multiverse_fpv2 <- plot_multiverse(all_hypo_drugs %>% filter(hypothesis == "favipiravir_hospital"))



# ((multiverse_hcq1 & hide_axis_theme) / (multiverse_hcq2 & hide_axis_theme) / (multiverse_az1 & hide_axis_theme) / (multiverse_az2 & hide_axis_theme) / (multiverse_fpv1 & hide_axis_theme) / (multiverse_fpv2)) + plot_layout(heights = c(1, 0.65,0.8,0.4, 0.8,0.3)) & theme(axis.title.y = element_blank())

fpv_lims <- c(-12,2)

multiverse_fpv1 <- plot_multiverse(all_hypo_drugs %>% filter(hypothesis == "favipiravir_death"), x_limits = fpv_lims, adjustments_to_hide = "FPV")

fig_drugs <- ((multiverse_hcq1 & hide_axis_theme) / (multiverse_az1 & theme(axis.title.x = element_blank())) / multiverse_fpv1)  + plot_layout(heights = c(0.8, 0.7,0.70)) & theme(axis.title.y = element_blank())

fig_drugs

ggsave(fig_drugs, filename = "supp_figs_production/fig5.tiff", width = 7.5, height = 7 * (7.5/8))

Published prognostic models are not better than using age as the sole predictor of outcome

Following Wynants et al.[@http://zotero.org/users/5567156/items/TS57HJ3T] we found five prediction models we were able to recompute: Li et al. report the ACP index [@http://zotero.org/users/5567156/items/EN4GICJ6] combining CRP and age to form 3 grades, Chen & Liu [@http://zotero.org/users/5567156/items/4DDUVEKS] report a continuous score using age, CRP, D-dimer and lymphocyte count, Shi et al. [@http://zotero.org/users/5567156/items/X4GWJIP5] use age, sex and hypertension to form 4 grades, Caramelo et al. use age, sex, hypertension, diabetes, cardiac disease, chronic respiratory disease and cancer to form a continuous score [@http://zotero.org/users/5567156/items/3A2RN4ZS], Bello-Chavolla et al. [@http://zotero.org/users/5567156/items/INZVABQI] use age, diabetes, obesity, pneumonia, chronic kidney disease, COPD and immunosuppression to build a score ranging from -6 to 22. For the latter two scores we had to impute some of the predictors as they had no immediate equivalent in our dataset. The outcomes present in the studies were: 12-day mortality, 30-day mortality and mortality without any further details, here we report results for both 12-day and 30-day mortality. Full details on the scores and how we used our dataset to compute them is given in the S3 file.

All prognostic models we tested performed similarly to or notably worse than using age as the only predictor and also worse than originally reported (Fig \@ref(fig:riskscores)). Additionally, some publications did not provide enough detail to unambiguously reconstruct how the score and/or outcome was assessed. We thus concur with Wynants et al. [@http://zotero.org/users/5567156/items/TS57HJ3T] that reported prediction scores are at high risk of bias and need additional careful evaluation.

col_types_auc <- cols(
  auc = col_double(),
  auc_low = col_double(),
  auc_high = col_double(),
  outcome = col_character(),
  score = col_character(),
  note = col_character(),
  subgroup = col_character()
)
all_auc <- read_csv(here::here("manuscript", "all_auc.csv"), col_types = col_types_auc, na = "NA") %>% filter(subgroup != "(admission)")
all_auc_orig <- read_csv(here::here("manuscript", "all_auc_orig.csv"), col_types = col_types_auc, na = "NA") %>% filter(subgroup != "(admission)")

fig_risk <- compare_auc_to_orig(all_auc_orig,
                    all_auc, show_outcome = FALSE)

fig_risk
ggsave(fig_risk, filename = "supp_figs_production/fig6.tiff", width = 7.5)

Discussion

Our data show the extent of between-patient variability in progression of the disease in terms of both length of hospital stay, duration of various types of breathing support and basic markers. A direct comparison with other studies is hard to perform as almost always only summaries of measurements are reported.

For multiple candidate Covid-19 treatments, observational data repeatedly contradicted results of randomized controlled trials (contrast e.g. Catteau et al. [@http://zotero.org/users/5567156/items/WZPNBRYN] to the RECOVERY trial [@http://zotero.org/users/5567156/items/9X3UH8ZM] for hydroxychloroquine and Liu et al. [@http://zotero.org/users/5567156/items/ES7EXX6D] to Agarwal et al. [@http://zotero.org/users/5567156/items/UEXLDJCB] for convalescent plasma). Our results for hydroxychloroquine also fit into this pattern. This should make us wary about over-interpreting the results of this study for azithromycin and favipiravir, although some higher-quality evidence that suggests clinical benefit of favipiravir has been reported [@http://zotero.org/users/5567156/items/YP62P6XJ].

The current (April 2021) Covid-19 treatment guidelines in Czech Republic recommend monoclonal antibodies and in some cases convalescent plasma or favipiravir as early treatment for high-risk patients with mild or no symptoms. For more severe cases dexamethason and anticoagulants are recommended while remdesivir is recommended only for patients that have severe disease but do not require mechanical ventilation [@http://zotero.org/users/5567156/items/YA3UAMVY]. This is similar to recommendations from the National Institute of Health in the USA who additionally recommend tocilizumab in some cases while not recommending convalescent plasma and favipiravir [@http://zotero.org/users/5567156/items/PCJBPPW7]. We do not believe our results should directly inform clinical decision, though we see some potential for inclusion of our results in future meta-analyses.

Regarding methodology, there are multiple approaches that are - at least in principle - capable of deriving causal conclusions from observational data, most notably the DAG framework [@http://zotero.org/users/5567156/items/IERDIX65; @http://zotero.org/users/5567156/items/HUSZBEA6] and target trial framework [@http://zotero.org/users/5567156/items/XRTF2VCR; @http://zotero.org/users/5567156/items/WK8XQTZQ]. In all approaches - and also in some randomized designs - it is common that substantial uncertainty about the best statistical model for the task at hand remains and cannot be eliminated. Nevertheless, most published papers present results only from a single statistical model. We believe that this uncertainty about model choice should not be ignored, rather we should embrace the uncertainty and employ a multiverse analysis or other forms of robustness checks to explore how our conclusions would differ had different assumptions been adopted. In this work we tried to show how such an analysis could be performed and reported in practice. We note that modelling choices that are often made semi-arbitrarily, e.g., logistic regression for survival at 28 days vs. a Cox proportional hazards model for time to event, did in our case lead to substantially different results.

The hidden-Markov models (HMMs) used in this study are of some interest because they fitted the data well and allowed for inclusion of a larger number of predictors than the Cox proportional hazards model without making the posterior uncertainty unreasonably large. We believe this is because such models make better use of the available detailed data. Additionally, HMMs can be used even when the outcomes are observed only indirectly or noisily --- as in the original study that motivated our models which concerns the progression of Alzheimer's disease [@http://zotero.org/users/5567156/items/F9K2L5Q8]. Noisily observed outcomes can pose problems for the proportional hazards model and require some special care [@http://zotero.org/users/5567156/items/YM2LGA3Y; @http://zotero.org/users/5567156/items/RTERDC3X]. We should however note that the Markov property assumed in HMMs is likely to be a reasonable approximation in fewer settings than the proportional hazards assumption of the Cox model.

Common problems with prognostic models in medicine are small sample sizes used to develop the models, weak or problematic statistical methods and lack of external validation on independent datasets [@http://zotero.org/users/5567156/items/E76N3Y5U]. Those problems are prevalent also in prediction models for Covid-19 [@http://zotero.org/users/5567156/items/TS57HJ3T]. We used our dataset to validate several models and observed very poor performance for four out of the five models tested. In fact, the simple assumption that older patients tend to have worse outcomes provides better or similar results to all of the models we were able to implement. This is despite all of the scores including age as a predictor. There seem to be two causes - three of the models dichotomize the age into just two groups which is known to lose information [@http://zotero.org/users/5567156/items/59L7MJPE; @http://zotero.org/users/5567156/items/N4552Y8Q]. Of the other two models Chen & Liu [@http://zotero.org/users/5567156/items/4DDUVEKS] use stepwise variable selection which is known to be a problematic approch [@http://zotero.org/users/5567156/items/J5VWRSEJ]. The resulting model puts largest relative weight on laboratory markers and deemphasizes age. Caramelo et al. [@http://zotero.org/users/5567156/items/3A2RN4ZS] take the decade of age as a very strong predictor and perform the best on our data. Still their results are not distinguishable from just using age. Our findings are almost the same as in a similar but larger validation study using 22 models and 411 patients from the United Kingdom where no tested model provided better prediction for mortality than age alone [@http://zotero.org/users/5567156/items/9XLSLVWJ].

Conclusions

We provide very weak observational evidence against a substantial beneficiary effect of using azithromycin (both with or without hydroxychloroquine) and against substantial negative effect of using favipiravir in hospitalized Covid-19 patients. We also observed better outcomes associated with taking hydroxychloroquine, which is likely linked to substantial confounding by indication. Where our results contradict randomized trials, the most likely explanation is systematic bias in our dataset.

A lesson from our analysis is that the assessment of treatment efficacy from observational data is sensitive to modelling assumptions while it is usually almost impossible to determine which of the models is more likely to reflect reality (if any). We believe that using multiverse analysis is an appropriate way to explore data in such contexts as it lets us be transparent about this sensitivity. We further believe that using hidden Markov models is a promising complement to the standard Cox proportional hazards analysis when detailed information on disease progression is available, particularly because it lets us impose additional structure on the model and thus make inferences with more disease states than would be possible to handle in the Cox framework, making better use of the available data.

Additionally, our experience indicates that a substantial fraction of published prognostic models will perform much worse on new patients than on the datasets they were built for and that external validation is crucial. We suggest that comparing the prognostic models against simple baselines (e.g., decade of age as the single predictor) should be a first step in validation. Furthermore, some of the published scores lack enough information to let others implement the score in the same way.

Supporting information

References



cas-bioinf/covid19retrospective documentation built on Sept. 7, 2021, 6:19 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close