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

knitr::opts_chunk$set(echo = FALSE)
devtools::load_all()
library(tidyverse)
library(survival)

theme_set(cowplot::theme_cowplot())
data <- read_data_for_analysis()
wide <- data$marker_data_wide


hypothesis_res_list <- list()

cache_dir <- here::here("local_temp_data")
if(!dir.exists(cache_dir)) {
  dir.create(cache_dir)
}

Competing risks proportional hazards models

In all of those models we use the coxph function from the Survival package and treat death and discharged as competing risks. The status column in the data corresponds to the three-valued state "Hospitalized" (1) / "Discharged" (2) / "Death" (3). I.e. in all model outputs below, there are two transitions "1:2" is the transition to "Discharged" and "1:3" is the transition to "Death". We evaluate the hypotheses about treatment effects or marker associations directly via the estimated log hazard ratio of the treatment for the relevant transition. Notably, the sign of the effect on the transition to "Discharged" is flipped for the multiverse analysis because it assumes that positive effect means increased hospital stay.

We do not work with models with more than about 5 predictors, as those become impossible to estimate due to insufficient amount of events.

data_surv <- wide %>%  
  filter(day > 0) %>%
  mutate(daym1 = day - 1,
    status = fct_collapse(breathing, Hospitalized = breathing_levels) %>% factor(ordered = FALSE) %>% fct_relevel("Hospitalized", "Discharged", "Death") %>% fct_explicit_na("Hospitalized")
  ) %>%
  inner_join(data$patient_data, by = c("hospital_id", "patient_id")) 

Treatment effects

Age + sex + hydroxychloroquine

The model uses the formula:

f_competing_hcq <- Surv(daym1, day, status) ~ age + sex + took_hcq
f_competing_hcq

Fitted coefficients below:

fit_competing_hcq <- coxph(f_competing_hcq, data = data_surv, id = patient_id)
fit_competing_hcq

hypothesis_res_list[["competing_hcq"]] <-
  rbind(
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_death, fit_competing_hcq,
                                        coefficient_name = "took_hcqTRUE", transition = "1:3", adjusted = "age, sex"),
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_hospital, fit_competing_hcq,
                                        coefficient_name = "took_hcqTRUE", transition = "1:2", adjusted = "age, sex", multiplier = -1)
  )

#hypothesis_res_list[["competing_hcq"]]

Age + sex + hydroxychloroquine + azithromycin

The formula becomes:

f_competing_hcq_az <- update.formula(f_competing_hcq, . ~ . + took_az)
f_competing_hcq_az

The fitted coefficients are:

fit_competing_hcq_az <- coxph(f_competing_hcq_az, data = data_surv, id = patient_id)

fit_competing_hcq_az

hypothesis_res_list[["competing_hcq_az"]] <-
  rbind(
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_death, fit_competing_hcq_az,
                                        coefficient_name = "took_hcqTRUE", transition = "1:3", adjusted = "age, sex, az"),
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_hospital, fit_competing_hcq_az,
                                        coefficient_name = "took_hcqTRUE", transition = "1:2", adjusted = "age, sex, az", multiplier = -1),
  frequentist_hypothesis_res_from_coxph(hypotheses$az_death, fit_competing_hcq_az,
                                        coefficient_name = "took_azTRUE", transition = "1:3", adjusted = "age, sex, hcq"),
  frequentist_hypothesis_res_from_coxph(hypotheses$az_hospital, fit_competing_hcq_az,
                                        coefficient_name = "took_azTRUE", transition = "1:2", adjusted = "age, sex, hcq", multiplier = -1)
  )
#hypothesis_res_list[["competing_hcq_az"]]

Age + sex + hydroxychloroquine + azithromycin + favipiravir

The formula is:

f_competing_hcq_az_fpv <- update.formula(f_competing_hcq_az, . ~ . + took_favipiravir)
f_competing_hcq_az_fpv

And the fitted model is:

fit_competing_hcq_az_fpv <- coxph(f_competing_hcq_az_fpv, data = data_surv, id = patient_id)

fit_competing_hcq_az_fpv

hypothesis_res_list[["competing_hcq_az_fpv"]] <-
  rbind(
  frequentist_hypothesis_res_from_coxph(hypotheses$favipiravir_death, fit_competing_hcq_az_fpv,
                                        coefficient_name = "took_favipiravirTRUE", transition = "1:3", adjusted = "age, sex, hcq, az", model_check = "Problematic"),
  frequentist_hypothesis_res_from_coxph(hypotheses$favipiravir_hospital, fit_competing_hcq_az_fpv,
                                        coefficient_name = "took_favipiravirTRUE", transition = "1:2", adjusted = "age, sex, hcq, az", model_check = "Problematic", multiplier = -1)
  )
#hypothesis_res_list[["competing_hcq_az_fpv"]]

Age + sex + hydroxychloroquine + azithromycin + admitted

The formula is:

f_competing_hcq_az_admitted <- update.formula(f_competing_hcq_az, . ~ . + admitted_for_covid)
f_competing_hcq_az_admitted

The fitted coefficients are:

fit_competing_hcq_az_admitted <- coxph(f_competing_hcq_az_admitted, data = data_surv, id = patient_id)

fit_competing_hcq_az_admitted

hypothesis_res_list[["competing_hcq_az_admitted"]] <-
  rbind(
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_death, fit_competing_hcq_az_admitted,
                                        coefficient_name = "took_hcqTRUE", transition = "1:3", adjusted = "age, sex, admitted, az"),
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_hospital, fit_competing_hcq_az_admitted,
                                        coefficient_name = "took_hcqTRUE", transition = "1:2", adjusted = "age, sex, admitted, az", multiplier = -1),
  frequentist_hypothesis_res_from_coxph(hypotheses$az_death, fit_competing_hcq_az_admitted,
                                        coefficient_name = "took_azTRUE", transition = "1:3", adjusted = "age, sex, admitted, hcq"),
  frequentist_hypothesis_res_from_coxph(hypotheses$az_hospital, fit_competing_hcq_az_admitted,
                                        coefficient_name = "took_azTRUE", transition = "1:2", adjusted = "age, sex, admitted, hcq", multiplier = -1)
  )
#hypothesis_res_list[["competing_hcq_az_admitted"]]

Age + sex + hydroxychloroquine + azithromycin + hospital

The new formula is:

f_competing_hcq_az_hospital <- update.formula(f_competing_hcq_az, . ~ . + strata(hospital_id))
f_competing_hcq_az_hospital

The resulting coefficients are:

data_surv_hospital <- data_surv %>% 
  mutate(hospital_id = if_else(hospital_id %in% c("DJkhZ", "lRBHa", "nDINR"), "small", hospital_id))
fit_competing_hcq_az_hospital <- coxph(f_competing_hcq_az_hospital, data = data_surv_hospital, id = patient_id)
fit_competing_hcq_az_hospital

hypothesis_res_list[["competing_hcq_az_hospital"]] <-
  rbind(
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_death, fit_competing_hcq_az_hospital,
                                        coefficient_name = "took_hcqTRUE", transition = "1:3", adjusted = "age, sex, az, hospital"),
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_hospital, fit_competing_hcq_az_hospital,
                                        coefficient_name = "took_hcqTRUE", transition = "1:2", adjusted = "age, sex, az, hospital", multiplier = -1),
  frequentist_hypothesis_res_from_coxph(hypotheses$az_death, fit_competing_hcq_az_hospital,
                                        coefficient_name = "took_azTRUE", transition = "1:3", adjusted = "age, sex, hcq, hospital"),
  frequentist_hypothesis_res_from_coxph(hypotheses$az_hospital, fit_competing_hcq_az_hospital,
                                        coefficient_name = "took_azTRUE", transition = "1:2", adjusted = "age, sex, hcq, hospital", multiplier = -1)
  )
#hypothesis_res_list[["competing_hcq_az_hospital"]]

Age + sex + hydroxychloroquine + azithromycin + comorbidities

The formula becomes:

f_competing_hcq_az_comorbidities_sum <- update.formula(f_competing_hcq_az, . ~ . + comorbidities_sum)
f_competing_hcq_az_comorbidities_sum

Model fit below:

fit_competing_hcq_az_comorbidities_sum <- coxph(f_competing_hcq_az_comorbidities_sum, data = data_surv, id = patient_id)
fit_competing_hcq_az_comorbidities_sum

hypothesis_res_list[["competing_hcq_az_comorbidities_sum"]] <-
  rbind(
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_death, fit_competing_hcq_az_comorbidities_sum,
                                        coefficient_name = "took_hcqTRUE", transition = "1:3", adjusted = "age, sex, az, comorbidities (sum)"),
  frequentist_hypothesis_res_from_coxph(hypotheses$hcq_hospital, fit_competing_hcq_az_comorbidities_sum,
                                        coefficient_name = "took_hcqTRUE", transition = "1:2", adjusted = "age, sex, az, comorbidities (sum)", multiplier = -1)
  )
#hypothesis_res_list[["competing_hcq_az_comorbidities_sum"]]

Marker effects

Age + sex + IL-6

The formula is:

f_competing_marker_base <- Surv(daym1, day,status) ~ age + sex
f_competing_il_6 <- update.formula(f_competing_marker_base, . ~ . + log_peak_IL_6)
f_competing_il_6

where log_peak_IL_6 is the logarithm of the peak IL-6 value seen up to the given time. The fitted model is:

data_IL_6 <- data_surv %>%
  group_by(patient_id) %>%
  filter(any(!is.na(IL_6))) %>%
  ungroup() %>%
  compute_marker_peak(column_name = "IL_6", new_column_name = "peak_IL_6", initial_value = 1) %>%
  mutate(log_peak_IL_6 = log(peak_IL_6))


fit_competing_il_6 <- coxph(f_competing_il_6, data = data_IL_6, id = patient_id, control = coxph.control(iter.max = 200))
fit_competing_il_6

hypothesis_res_list[["competing_il_6"]] <-
  frequentist_hypothesis_res_from_coxph(hypotheses$IL_6_death, fit_competing_il_6,
                                        coefficient_name = "log_peak_IL_6", transition = "1:3", adjusted = "age, sex", model_check = "Suspicious")
#hypothesis_res_list[["competing_il_6"]]

Note the warning on convergence and the huge se for the log_peak_IL_6 coefficient for the 1:2 transition. We mark this model as "Suspicious".

Age + sex + D-dimer

The formula is:

f_competing_d_dimer <- update.formula(f_competing_marker_base, . ~ . + log_peak_d_dimer)
f_competing_d_dimer

where log_peak_d_dimer is the logarithm of the peak D-dimer value seen up to the given time. The fitted model is:

data_d_dimer <- data_surv %>%
  group_by(patient_id) %>%
  filter(any(!is.na(d_dimer))) %>%
  ungroup() %>%
  compute_marker_peak(column_name = "d_dimer", new_column_name = "peak_d_dimer", initial_value = 0) %>%
  mutate(log_peak_d_dimer = log(peak_d_dimer + 1) - 4)


fit_competing_d_dimer <- coxph(f_competing_d_dimer, data = data_d_dimer, id = patient_id, control = coxph.control(iter.max = 2000))
fit_competing_d_dimer

hypothesis_res_list[["competing_d_dimer"]] <-
  frequentist_hypothesis_res_from_coxph(hypotheses$d_dimer_death, fit_competing_d_dimer,
                                        coefficient_name = "log_peak_d_dimer", transition = "1:3", adjusted = "age, sex", model_check = "Suspicious")
#hypothesis_res_list[["competing_d_dimer"]]

We see similar, but even more extreme convergence issues as before, we mark the model "Suspicious".

hypothesis_res_all <- do.call(rbind, hypothesis_res_list)
write_csv(hypothesis_res_all, path = here::here("manuscript", "/competing_risks_res.csv"))

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