In cas-bioinf/covid19retrospective: What the Package Does (One Line, Title Case)
Bayesian models with brms
knitr::opts_chunk$set(cache = TRUE, echo=FALSE, fig.width = 10)
devtools::load_all()
library(tidyverse)
library(tidybayes)
library(brms)
options(width = 160)
options(brms.backend = "cmdstanr")
options(mc.cores = 4)
#theme_set(bayesplot::theme_default())
cache_dir <- here::here("local_temp_data", "analysis_brms")
if (!dir.exists(cache_dir)) {
dir.create(cache_dir, recursive = TRUE)
}
scale2 <- function(x, na.rm = TRUE) {
(x - mean(x, na.rm = na.rm)) / sd(x, na.rm = na.rm)
}
center <- function(x, na.rm = TRUE) {
x - mean(x, na.rm = na.rm)
}
SW <- suppressWarnings
sum_coding <- function(x, lvls = levels(x)) {
# codes the first category with -1
nlvls <- length(lvls)
stopifnot(nlvls > 1)
cont <- diag(nlvls)[, -nlvls, drop = FALSE]
cont[nlvls, ] <- -1
cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE]
colnames(cont) <- lvls[-1]
x <- factor(x, levels = lvls)
contrasts(x) <- cont
x
}
dummy_coding <- function(x, lvls = levels(x)) {
x <- factor(x, levels = lvls)
contrasts(x) <- contr.treatment(levels(x))
x
}
default_control <- list(adapt_delta = 0.95)
Here we use a set of Bayesian linear models using the brms package. There are several classes of model:
- Binary survival models - logistic regression with survival as the outcome.
- Survival time models - Bayesian version of the single-outcome Cox proportional hazards model for the "Death" event.
- Binary discharged models - logistic regression with discharged as the outcome (differs from binary survival in the handling of censored patients)
- Discharged time models - Bayesian version of the single-outcome Cox proportional hazards model for the "Discharged" event
- Categorical models - categorical (multinomial) regression of the "Discharged"/"Hospitalized"/"Death" outcome.
For the final outcome models (binary and categorical) we use three variants of dealing with censoring:
- All: Remove censored patients and keep only those that have data until final outcome
- 7 days: Remove patients with less than 7 days of data and treat the state at day 7 post admission as the outcome
- 28 days: Same as 7 days, but with a longer time frame
The 7 and 28 days models should be somewhat less biased as censored patients are more likely those that spent longer time in hospital and thus have likely more severe disease. We marked the "All" models as "Suspicious" for the purpose of the multiverse analysis.
For each model we show the brms formula used, summary of fitted coefficients. See the brms manual for detailed explanation of each model. Additionally we show some posterior predictive checks [@http://zotero.org/users/5567156/items/42DSAQX9; @http://zotero.org/users/5567156/items/ZKWU8INU] to verify that models represent the data well. Consult the bayesplot package manual for details on the checks used. As another way to check for potential overfitting we use the loo package [@http://zotero.org/users/5567156/items/88ANHBVF] and its Pareto-k diagnostic.
All hypotheses are evaluated via model coefficients of predictors representing treatments, except for sign corrections (e.g. in binary survival models positive effects represent increased survival while for the main analysis we take positive effect as increased mortality, so we flip the sign).
Some derived predictors are used age_norm (age, normalized to have roughly mean 0 and standard deviation of 1), ever_az, ever_hcq, ever_favipiravir (whether the patient took the respective medications at any time) and comorbidities_sum - the sum of comorbidities as described in Section \@ref(comorbiditiessum).
data <- read_data_for_analysis()
data$patient_data <- data$patient_data %>%
mutate(
survival = outcome != "Death",
discharged = outcome == "Discharged",
outcome3 = fct_collapse(outcome, Discharged = c("Discharged", "Transferred")) %>%
fct_relevel("Hospitalized", "Discharged", "Death")
)
data$surv_data <- data$marker_data_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"),
censored = case_when(
status == "Death" ~ "none",
status %in% c("Hospitalized", "Discharged") ~ "right"
),
censored_discharged = case_when(
status == "Discharged" ~ "none",
# approximation since patients who died cannot be
# discharged in the future
status %in% c("Hospitalized", "Death") ~ "right"
)
) %>%
inner_join(data$patient_data, by = c("hospital_id", "patient_id"))
data$surv_data_last_day <- data$surv_data %>%
group_by(patient_id) %>%
filter(day == max(day)) %>%
ungroup()
# we need to have "state at day 7", i.e. patients that died/were discharged
# after day 7 are still treated as "Hospitalized" for the analysis. Other
# studies (e.g. RECOVERY) work this way. Code provided by Martin.
surv_data_at_day <- function(patient_data, day) {
res <- patient_data %>%
filter(!(last_record < day & outcome3 == "Hospitalized"))
# Non tidy to avoid hassle with factor levels
res$outcome3[res$last_record > day] <- "Hospitalized"
res
}
data$surv_data_7 <- surv_data_at_day(data$surv_data, 7)
data$surv_data_28 <- surv_data_at_day(data$surv_data, 28)
hypothesis_res_list <- list()
Binary survival models
All variants of the model use the same brms formula, with family bernoulli("logit").
form_survival_base <- bf(
survival ~ sex + age_norm +
smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir +
(1 | hospital_id),
family = bernoulli("logit")
)
print_formula(form_survival_base)
Basic binary model
prior_survival_base <-
prior(normal(0, 5), class = "b") +
prior(normal(0, 5), class = "sd")
fit_survival_base <- brm(
formula = form_survival_base,
data = data$patient_data,
prior = prior_survival_base,
control = default_control,
file = paste0(cache_dir, "/fit_survival_base")
)
fit_survival_base <- add_criterion(fit_survival_base, "loo")
The fitted coefficients:
summary(fit_survival_base)
invisible(loo(fit_survival_base)) #loo used primarily for checks not for the actual values
The error_binned posterior predictive check - showing the distribution of errors per predicted probability bin.
pp_check(fit_survival_base, "error_binned", nsamples = 20)
#Marking the results as suspicious because of the potential bias
draws <- tidy_draws(fit_survival_base)
hypothesis_res_list[["hcq_az_survival_base"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_hcqTRUE,
model = "brms_binary",
estimand = "log(OR)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_azTRUE,
model = "brms_binary",
estimand = "log(OR)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_favipiravirTRUE,
model = "brms_binary",
estimand = "log(OR)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
)
)
Basic binary model (7 days)
fit_survival_base_7 <- brm(
formula = form_survival_base,
data = data$surv_data_7,
prior = prior_survival_base,
control = default_control,
file = paste0(cache_dir, "/fit_survival_base_7")
)
fit_survival_base_7 <- add_criterion(fit_survival_base_7, "loo")
The fitted coefficients:
summary(fit_survival_base_7)
# ce <- conditional_effects(fit_survival_base_7)
# plot(ce, ask = FALSE)
The error_binned check:
pp_check(fit_survival_base_7, "error_binned", nsamples = 20)
invisible(loo(fit_survival_base_7))
draws <- tidy_draws(fit_survival_base_7)
hypothesis_res_list[["hcq_az_survival_base_7"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_hcqTRUE,
model = "brms_binary_7",
estimand = "log(OR)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_azTRUE,
model = "brms_binary_7",
estimand = "log(OR)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_favipiravirTRUE,
model = "brms_binary_7",
estimand = "log(OR)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
)
)
Basic binary model (28 days)
fit_survival_base_28 <- brm(
formula = form_survival_base,
data = data$surv_data_28,
prior = prior_survival_base,
control = default_control,
file = paste0(cache_dir, "/fit_survival_base_28")
)
fit_survival_base_28 <- add_criterion(fit_survival_base_28, "loo")
The fitted coefficients:
summary(fit_survival_base_28)
# ce <- conditional_effects(fit_survival_base_28)
# plot(ce, ask = FALSE)
Note that there is a divergent transition indicating potential computational problems, we will mark the model as "Suspicious".
The error_binned check:
pp_check(fit_survival_base_28, "error_binned", nsamples = 20)
invisible(loo(fit_survival_base_28))
draws <- tidy_draws(fit_survival_base_28)
hypothesis_res_list[["hcq_az_survival_base_28"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_hcqTRUE,
model = "brms_binary_28",
estimand = "log(OR)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_azTRUE,
model = "brms_binary_28",
estimand = "log(OR)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_favipiravirTRUE,
model = "brms_binary_28",
estimand = "log(OR)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
)
)
Survival time models
prior_survival_time <-
prior(normal(0, 5), class = "b") +
prior(normal(0, 5), class = "sd")
Fewer adjustments
Those models are described by the brms formula with family brmsfamily("cox", "log"):
form_survival_time_small <- bf(
day | cens(censored) ~ sex + age_norm + ever_hcq + ever_az + ever_favipiravir +
(1 | hospital_id),
family = brmsfamily("cox", "log")
)
print_formula(form_survival_time_small)
fit_survival_time_small <- brm(
formula = form_survival_time_small,
data = data$surv_data_last_day,
prior = prior_survival_time,
inits = 0,
control = default_control,
file = paste0(cache_dir, "/fit_survival_time_small")
)
fit_survival_time_small <- add_criterion(fit_survival_time_small, "loo")
Here are the fitted coefficients:
summary(fit_survival_time_small)
invisible(loo(fit_survival_time_small))
draws <- tidy_draws(fit_survival_time_small)
hypothesis_res_list[["hcq_az_survival_time_small"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_hcqTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm",
model_check = "some influential observations"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_azTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm",
model_check = "some influential observations"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_favipiravirTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm",
model_check = "some influential observations"
)
)
Fewer adjustments - First wave only
The same model as above, but using data from the first wave only. Here is the summary of the fit:
fit_survival_time_small_fw <- brm(
formula = form_survival_time_small,
data = data$surv_data_last_day %>% filter(first_wave),
prior = prior_survival_time,
inits = 0,
control = default_control,
file = paste0(cache_dir, "/fit_survival_time_small_fw")
)
fit_survival_time_small_fw <- add_criterion(fit_survival_time_small_fw, "loo")
summary(fit_survival_time_small_fw)
loo(fit_survival_time_small_fw)
draws <- tidy_draws(fit_survival_time_small_fw)
hypothesis_res_list[["hcq_az_survival_time_small_fw"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_hcqTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, first_wave",
model_check = "some influential observations"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_azTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, first_wave",
model_check = "some influential observations"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_favipiravirTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, first_wave",
model_check = "some influential observations"
)
)
More adjustments
For those models the formula changes to:
form_survival_time <- bf(
day | cens(censored) ~ sex + age_norm +
smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir +
(1 | hospital_id),
family = brmsfamily("cox", "log")
)
print_formula(form_survival_time)
fit_survival_time <- brm(
formula = form_survival_time,
data = data$surv_data_last_day,
prior = prior_survival_time,
inits = 0,
control = default_control,
file = paste0(cache_dir, "/fit_survival_time")
)
fit_survival_time <- add_criterion(fit_survival_time, "loo")
Here are the fitted coefficients - there is a divergent transition so we are possibly fitting an overparemtrized model and we will mark it as "Suspicious":
summary(fit_survival_time)
In addition, the LOO diagnostics warn us that there were some very influential observations, indicating possible overfitting.
loo(fit_survival_time)
draws <- tidy_draws(fit_survival_time)
hypothesis_res_list[["hcq_az_survival_time"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_hcqTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_azTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_favipiravirTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
)
)
More adjustments - First wave only
The same model as above, but fitting only to patients from first wave. Here are the fitted model coefficients:
fit_survival_time_fw <- brm(
formula = form_survival_time,
data = data$surv_data_last_day %>% filter(first_wave),
prior = prior_survival_time,
inits = 0,
control = default_control,
file = paste0(cache_dir, "/fit_survival_time_fw")
)
fit_survival_time_fw <- add_criterion(fit_survival_time_fw, "loo")
summary(fit_survival_time_fw)
We do not see divergent transitions, but some overly influential observations (according to loo) remain, the model will be marked as "Suspicious".
loo(fit_survival_time_fw)
draws <- tidy_draws(fit_survival_time_fw)
hypothesis_res_list[["hcq_az_survival_time_fw"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_hcqTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_azTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_ever_favipiravirTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave",
model_check = "Suspicious"
)
)
Binary discharged models
All the models in this section use the bernoulli("logit") family (i.e. logistic regression) and follow the brms formula
form_discharged_base <- bf(
discharged ~ sex + age_norm +
smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir +
(1 | hospital_id),
family = bernoulli("logit")
)
print_formula(form_discharged_base)
Basic binary model
prior_discharged_base <-
prior(normal(0, 5), class = "b") +
prior(normal(0, 5), class = "sd")
fit_discharged_base <- brm(
formula = form_discharged_base,
data = data$patient_data,
prior = prior_discharged_base,
control = default_control,
file = paste0(cache_dir, "/fit_discharged_base")
)
fit_discharged_base <- add_criterion(fit_discharged_base, "loo")
Here is the summary of the fitted coefficients
summary(fit_discharged_base)
invisible(loo(fit_discharged_base))
# ce <- conditional_effects(fit_discharged_base)
# plot(ce, ask = FALSE)
And the error_binned check:
pp_check(fit_discharged_base, "error_binned", nsamples = 20)
draws <- tidy_draws(fit_discharged_base)
hypothesis_res_list[["hcq_az_discharged_base"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_hcqTRUE,
model = "brms_binary",
estimand = "log(OR)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_azTRUE,
model = "brms_binary",
estimand = "log(OR)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_favipiravirTRUE,
model = "brms_binary",
estimand = "log(OR)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
)
)
Basic binary model (7 Days)
fit_discharged_base_7 <- brm(
formula = form_discharged_base,
data = data$surv_data_7,
prior = prior_discharged_base,
control = default_control,
file = paste0(cache_dir, "/fit_discharged_base_7")
)
fit_discharged_base_7 <- add_criterion(fit_discharged_base_7, "loo")
Here is the summary of the fitted coefficients:
summary(fit_discharged_base_7)
invisible(loo(fit_discharged_base_7))
# ce <- conditional_effects(fit_discharged_base_7)
# plot(ce, ask = FALSE)
And the error_binned check.
pp_check(fit_discharged_base_7, "error_binned", nsamples = 20)
draws <- tidy_draws(fit_discharged_base_7)
hypothesis_res_list[["hcq_az_discharged_base_7"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_hcqTRUE,
model = "brms_binary_7",
estimand = "log(OR)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_azTRUE,
model = "brms_binary_7",
estimand = "log(OR)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_favipiravirTRUE,
model = "brms_binary_7",
estimand = "log(OR)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
)
)
Basic binary model (28 Days)
fit_discharged_base_28 <- brm(
formula = form_discharged_base,
data = data$surv_data_28,
prior = prior_discharged_base,
control = default_control,
file = paste0(cache_dir, "/fit_discharged_base_28")
)
fit_discharged_base_28 <- add_criterion(fit_discharged_base_28, "loo")
Here are the fitted coefficients:
summary(fit_discharged_base_28)
invisible(loo(fit_discharged_base_28))
# ce <- conditional_effects(fit_discharged_base_28)
# plot(ce, ask = FALSE)
And the error_binned check:
pp_check(fit_discharged_base_28, "error_binned")
draws <- tidy_draws(fit_discharged_base_28)
hypothesis_res_list[["hcq_az_discharged_base_28"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_hcqTRUE,
model = "brms_binary_28",
estimand = "log(OR)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_azTRUE,
model = "brms_binary_28",
estimand = "log(OR)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_favipiravirTRUE,
model = "brms_binary_28",
estimand = "log(OR)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
)
)
Discharged time model
The model uses the family brmsfamily("cox", "log") and the formula:
form_discharged_time <- bf(
day | cens(censored_discharged) ~ sex + age_norm +
smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir +
(1 | hospital_id),
family = brmsfamily("cox", "log")
)
print_formula(form_discharged_time)
prior_discharged_time <-
prior(normal(0, 5), class = "b") +
prior(normal(0, 5), class = "sd")
fit_discharged_time <- brm(
formula = form_discharged_time,
data = data$surv_data_last_day,
prior = prior_discharged_time,
inits = 0,
control = list(adapt_delta = 0.98),
file = paste0(cache_dir, "/fit_discharged_time")
)
fit_discharged_time <- add_criterion(fit_discharged_time, "loo")
The fitted coefficients are:
summary(fit_discharged_time)
The loo diagnostics warn us that there were some very influential observations.
loo(fit_discharged_time)
draws <- tidy_draws(fit_discharged_time)
hypothesis_res_list[["hcq_az_discharged_time"]] <- rbind(
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_hcqTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "some influential observations"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_azTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "some influential observations"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_ever_favipiravirTRUE,
model = "brms_cox",
estimand = "log(HR)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "some influential observations"
)
)
Categorical models
All the categorical models use the categorical("logit") family and the following formula:
form_outcome_base <- bf(
outcome3 ~ sex + age_norm +
smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir +
(1 | hospital_id),
family = categorical("logit")
)
print_formula(form_outcome_base)
Basic Categorical model
prior_outcome_base <-
prior(normal(0, 5), class = "b", dpar = "muDeath") +
prior(normal(0, 5), class = "sd", dpar = "muDeath") +
prior(normal(0, 5), class = "b", dpar = "muDischarged") +
prior(normal(0, 5), class = "sd", dpar = "muDischarged")
fit_outcome_base <- brm(
formula = form_outcome_base,
data = data$patient_data,
prior = prior_outcome_base,
control = default_control,
file = paste0(cache_dir, "/fit_outcome_base")
)
fit_outcome_base <- add_criterion(fit_outcome_base, "loo")
Fitted coefficients below:
summary(fit_outcome_base)
invisible(loo(fit_outcome_base))
# ce <- conditional_effects(fit_outcome_base)
# plot(ce, ask = FALSE)
draws <- tidy_draws(fit_outcome_base)
hypothesis_res_list[["hcq_az_outcome_base"]] <- rbind(
# Death
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_hcqTRUE,
model = "brms_categorical",
estimand = "logit(death)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_azTRUE,
model = "brms_categorical",
estimand = "logit(death)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_favipiravirTRUE,
model = "brms_categorical",
estimand = "logit(death)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
# Discharged
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_hcqTRUE,
model = "brms_categorical",
estimand = "logit(discharged)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_azTRUE,
model = "brms_categorical",
estimand = "logit(discharged)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_favipiravirTRUE,
model = "brms_categorical",
estimand = "logit(discharged)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "Suspicious"
)
)
Basic Categorical model (7 Days)
fit_outcome_base_7 <- brm(
formula = form_outcome_base,
data = data$surv_data_7,
prior = prior_outcome_base,
control = list(adapt_delta = 0.98),
file = paste0(cache_dir, "/fit_outcome_base_7")
)
fit_outcome_base_7 <- add_criterion(fit_outcome_base_7, "loo")
Summary of model fit:
summary(fit_outcome_base_7)
invisible(loo(fit_outcome_base_7))
# ce <- conditional_effects(fit_outcome_base_7)
# plot(ce, ask = FALSE)
#pp_check(fit_outcome_base_7, type = "bars", nsamples = 200)
draws <- tidy_draws(fit_outcome_base_7)
hypothesis_res_list[["hcq_az_outcome_base_7"]] <- rbind(
# Death
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_hcqTRUE,
model = "brms_categorical_7",
estimand = "logit(death)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_azTRUE,
model = "brms_categorical_7",
estimand = "logit(death)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_favipiravirTRUE,
model = "brms_categorical_7",
estimand = "logit(death)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
# Discharged
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_hcqTRUE,
model = "brms_categorical_7",
estimand = "logit(discharged)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_azTRUE,
model = "brms_categorical_7",
estimand = "logit(discharged)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_favipiravirTRUE,
model = "brms_categorical_7",
estimand = "logit(discharged)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
)
)
Basic Categorical model (7 Days), first wave only
The same model as above, using only data from the first wave. Summary below:
fit_outcome_base_7_fw <- brm(
formula = form_outcome_base,
data = data$surv_data_7 %>% filter(first_wave),
prior = prior_outcome_base,
control = default_control,
file = paste0(cache_dir, "/fit_outcome_base_7_fw")
)
fit_outcome_base_7_fw <- add_criterion(fit_outcome_base_7_fw, "loo")
summary(fit_outcome_base_7_fw)
invisible(loo(fit_outcome_base_7_fw))
# ce <- conditional_effects(fit_outcome_base_7_fw)
# plot(ce, ask = FALSE)
#pp_check(fit_outcome_base_7_fw, type = "bars", nsamples = 200)
draws <- tidy_draws(fit_outcome_base_7_fw)
hypothesis_res_list[["hcq_az_outcome_base_7_fw"]] <- rbind(
# Death
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_hcqTRUE,
model = "brms_categorical_7",
estimand = "logit(death)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_azTRUE,
model = "brms_categorical_7",
estimand = "logit(death)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_favipiravirTRUE,
model = "brms_categorical_7",
estimand = "logit(death)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
# Discharged
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_hcqTRUE,
model = "brms_categorical_7",
estimand = "logit(discharged)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_azTRUE,
model = "brms_categorical_7",
estimand = "logit(discharged)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_favipiravirTRUE,
model = "brms_categorical_7",
estimand = "logit(discharged)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
)
)
Basic Categorical model (28 Days)
fit_outcome_base_28 <- brm(
formula = form_outcome_base,
data = data$surv_data_28,
prior = prior_outcome_base,
control = default_control,
file = paste0(cache_dir, "/fit_outcome_base_28")
)
fit_outcome_bas_28 <- add_criterion(fit_outcome_base_28, "loo")
Summary of the model is:
summary(fit_outcome_base_28)
#pp_check(fit_outcome_base_28)
invisible(loo(fit_outcome_base_28))
draws <- tidy_draws(fit_outcome_base_28)
hypothesis_res_list[["hcq_az_outcome_base_28"]] <- rbind(
# Death
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_hcqTRUE,
model = "brms_categorical_28",
estimand = "logit(death)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_azTRUE,
model = "brms_categorical_28",
estimand = "logit(death)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_favipiravirTRUE,
model = "brms_categorical_28",
estimand = "logit(death)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
# Discharged
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_hcqTRUE,
model = "brms_categorical_28",
estimand = "logit(discharged)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_azTRUE,
model = "brms_categorical_28",
estimand = "logit(discharged)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_favipiravirTRUE,
model = "brms_categorical_28",
estimand = "logit(discharged)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity",
model_check = "probably overconfident"
)
)
Basic Categorical model (28 Days), first wave only
The same model, but only for patients from the first wave, summary of the fit:
fit_outcome_base_28_fw <- brm(
formula = form_outcome_base,
data = data$surv_data_28,
prior = prior_outcome_base,
control = default_control,
file = paste0(cache_dir, "/fit_outcome_base_28_fw")
)
fit_outcome_bas_28 <- add_criterion(fit_outcome_base_28_fw, "loo")
summary(fit_outcome_base_28_fw)
#pp_check(fit_outcome_base_28_fw)
invisible(loo(fit_outcome_base_28_fw))
draws <- tidy_draws(fit_outcome_base_28_fw)
hypothesis_res_list[["hcq_az_outcome_base_28_fw"]] <- rbind(
# Death
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_hcqTRUE,
model = "brms_categorical_28",
estimand = "logit(death)",
hypothesis = hypotheses$hcq_death,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_azTRUE,
model = "brms_categorical_28",
estimand = "logit(death)",
hypothesis = hypotheses$az_death,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = draws$b_muDeath_ever_favipiravirTRUE,
model = "brms_categorical_28",
estimand = "logit(death)",
hypothesis = hypotheses$favipiravir_death,
adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
# Discharged
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_hcqTRUE,
model = "brms_categorical_28",
estimand = "logit(discharged)",
hypothesis = hypotheses$hcq_hospital,
adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_azTRUE,
model = "brms_categorical_28",
estimand = "logit(discharged)",
hypothesis = hypotheses$az_hospital,
adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
),
bayesian_hypothesis_res_from_draws(
draws = -draws$b_muDischarged_ever_favipiravirTRUE,
model = "brms_categorical_28",
estimand = "logit(discharged)",
hypothesis = hypotheses$favipiravir_hospital,
adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave",
model_check = "probably overconfident"
)
)
Worst Breathing {.tabset}
Not reported in the main manuscript but we also experimented with ordinal (stopping ratio) models for the worst breathing level experienced by a patient.
form_worst_breathing_base <- bf(
worst_breathing ~ sex + age_norm +
smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir +
(1 | hospital_id),
family = sratio("logit")
)
prior_worst_breathing_base <-
prior(normal(0, 5), class = "b") +
prior(normal(0, 5), class = "sd")
fit_worst_breathing_base <- brm(
formula = form_worst_breathing_base,
data = data$patient_data,
prior = prior_worst_breathing_base,
control = default_control,
file = paste0(cache_dir, "/fit_worst_breathing_base")
)
fit_worst_breathing_base <- add_criterion(fit_worst_breathing_base, "loo")
summary(fit_worst_breathing_base)
invisible(loo(fit_worst_breathing_base))
# ce <- conditional_effects(fit_worst_breathing_base, categorical = TRUE)
# plot(ce, ask = FALSE)
hypothesis_res_all <- do.call(rbind, hypothesis_res_list)
write_csv(hypothesis_res_all, path = here::here("manuscript/analysis_brms.csv"))
cas-bioinf/covid19retrospective documentation built on Sept. 7, 2021, 6:19 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Bayesian models with brms
knitr::opts_chunk$set(cache = TRUE, echo=FALSE, fig.width = 10)
devtools::load_all() library(tidyverse) library(tidybayes) library(brms) options(width = 160) options(brms.backend = "cmdstanr") options(mc.cores = 4) #theme_set(bayesplot::theme_default()) cache_dir <- here::here("local_temp_data", "analysis_brms") if (!dir.exists(cache_dir)) { dir.create(cache_dir, recursive = TRUE) }
scale2 <- function(x, na.rm = TRUE) { (x - mean(x, na.rm = na.rm)) / sd(x, na.rm = na.rm) } center <- function(x, na.rm = TRUE) { x - mean(x, na.rm = na.rm) } SW <- suppressWarnings sum_coding <- function(x, lvls = levels(x)) { # codes the first category with -1 nlvls <- length(lvls) stopifnot(nlvls > 1) cont <- diag(nlvls)[, -nlvls, drop = FALSE] cont[nlvls, ] <- -1 cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE] colnames(cont) <- lvls[-1] x <- factor(x, levels = lvls) contrasts(x) <- cont x } dummy_coding <- function(x, lvls = levels(x)) { x <- factor(x, levels = lvls) contrasts(x) <- contr.treatment(levels(x)) x } default_control <- list(adapt_delta = 0.95)
Here we use a set of Bayesian linear models using the brms package. There are several classes of model:
- Binary survival models - logistic regression with survival as the outcome.
- Survival time models - Bayesian version of the single-outcome Cox proportional hazards model for the "Death" event.
- Binary discharged models - logistic regression with discharged as the outcome (differs from binary survival in the handling of censored patients)
- Discharged time models - Bayesian version of the single-outcome Cox proportional hazards model for the "Discharged" event
- Categorical models - categorical (multinomial) regression of the "Discharged"/"Hospitalized"/"Death" outcome.
For the final outcome models (binary and categorical) we use three variants of dealing with censoring:
- All: Remove censored patients and keep only those that have data until final outcome
- 7 days: Remove patients with less than 7 days of data and treat the state at day 7 post admission as the outcome
- 28 days: Same as 7 days, but with a longer time frame
The 7 and 28 days models should be somewhat less biased as censored patients are more likely those that spent longer time in hospital and thus have likely more severe disease. We marked the "All" models as "Suspicious" for the purpose of the multiverse analysis.
For each model we show the brms formula used, summary of fitted coefficients. See the brms manual for detailed explanation of each model. Additionally we show some posterior predictive checks [@http://zotero.org/users/5567156/items/42DSAQX9; @http://zotero.org/users/5567156/items/ZKWU8INU] to verify that models represent the data well. Consult the bayesplot package manual for details on the checks used. As another way to check for potential overfitting we use the loo package [@http://zotero.org/users/5567156/items/88ANHBVF] and its Pareto-k diagnostic.
All hypotheses are evaluated via model coefficients of predictors representing treatments, except for sign corrections (e.g. in binary survival models positive effects represent increased survival while for the main analysis we take positive effect as increased mortality, so we flip the sign).
Some derived predictors are used age_norm (age, normalized to have roughly mean 0 and standard deviation of 1), ever_az, ever_hcq, ever_favipiravir (whether the patient took the respective medications at any time) and comorbidities_sum - the sum of comorbidities as described in Section \@ref(comorbiditiessum).
data <- read_data_for_analysis()
data$patient_data <- data$patient_data %>% mutate( survival = outcome != "Death", discharged = outcome == "Discharged", outcome3 = fct_collapse(outcome, Discharged = c("Discharged", "Transferred")) %>% fct_relevel("Hospitalized", "Discharged", "Death") )
data$surv_data <- data$marker_data_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"), censored = case_when( status == "Death" ~ "none", status %in% c("Hospitalized", "Discharged") ~ "right" ), censored_discharged = case_when( status == "Discharged" ~ "none", # approximation since patients who died cannot be # discharged in the future status %in% c("Hospitalized", "Death") ~ "right" ) ) %>% inner_join(data$patient_data, by = c("hospital_id", "patient_id")) data$surv_data_last_day <- data$surv_data %>% group_by(patient_id) %>% filter(day == max(day)) %>% ungroup()
# we need to have "state at day 7", i.e. patients that died/were discharged # after day 7 are still treated as "Hospitalized" for the analysis. Other # studies (e.g. RECOVERY) work this way. Code provided by Martin. surv_data_at_day <- function(patient_data, day) { res <- patient_data %>% filter(!(last_record < day & outcome3 == "Hospitalized")) # Non tidy to avoid hassle with factor levels res$outcome3[res$last_record > day] <- "Hospitalized" res } data$surv_data_7 <- surv_data_at_day(data$surv_data, 7) data$surv_data_28 <- surv_data_at_day(data$surv_data, 28)
hypothesis_res_list <- list()
Binary survival models
All variants of the model use the same brms formula, with family bernoulli("logit").
form_survival_base <- bf( survival ~ sex + age_norm + smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir + (1 | hospital_id), family = bernoulli("logit") ) print_formula(form_survival_base)
Basic binary model
prior_survival_base <- prior(normal(0, 5), class = "b") + prior(normal(0, 5), class = "sd") fit_survival_base <- brm( formula = form_survival_base, data = data$patient_data, prior = prior_survival_base, control = default_control, file = paste0(cache_dir, "/fit_survival_base") ) fit_survival_base <- add_criterion(fit_survival_base, "loo")
The fitted coefficients:
summary(fit_survival_base) invisible(loo(fit_survival_base)) #loo used primarily for checks not for the actual values
The error_binned posterior predictive check - showing the distribution of errors per predicted probability bin.
pp_check(fit_survival_base, "error_binned", nsamples = 20)
#Marking the results as suspicious because of the potential bias draws <- tidy_draws(fit_survival_base) hypothesis_res_list[["hcq_az_survival_base"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_hcqTRUE, model = "brms_binary", estimand = "log(OR)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_azTRUE, model = "brms_binary", estimand = "log(OR)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_favipiravirTRUE, model = "brms_binary", estimand = "log(OR)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ) )
Basic binary model (7 days)
fit_survival_base_7 <- brm( formula = form_survival_base, data = data$surv_data_7, prior = prior_survival_base, control = default_control, file = paste0(cache_dir, "/fit_survival_base_7") ) fit_survival_base_7 <- add_criterion(fit_survival_base_7, "loo")
The fitted coefficients:
summary(fit_survival_base_7) # ce <- conditional_effects(fit_survival_base_7) # plot(ce, ask = FALSE)
The error_binned check:
pp_check(fit_survival_base_7, "error_binned", nsamples = 20) invisible(loo(fit_survival_base_7))
draws <- tidy_draws(fit_survival_base_7) hypothesis_res_list[["hcq_az_survival_base_7"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_hcqTRUE, model = "brms_binary_7", estimand = "log(OR)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_azTRUE, model = "brms_binary_7", estimand = "log(OR)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_favipiravirTRUE, model = "brms_binary_7", estimand = "log(OR)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ) )
Basic binary model (28 days)
fit_survival_base_28 <- brm( formula = form_survival_base, data = data$surv_data_28, prior = prior_survival_base, control = default_control, file = paste0(cache_dir, "/fit_survival_base_28") ) fit_survival_base_28 <- add_criterion(fit_survival_base_28, "loo")
The fitted coefficients:
summary(fit_survival_base_28) # ce <- conditional_effects(fit_survival_base_28) # plot(ce, ask = FALSE)
Note that there is a divergent transition indicating potential computational problems, we will mark the model as "Suspicious".
The error_binned check:
pp_check(fit_survival_base_28, "error_binned", nsamples = 20) invisible(loo(fit_survival_base_28))
draws <- tidy_draws(fit_survival_base_28) hypothesis_res_list[["hcq_az_survival_base_28"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_hcqTRUE, model = "brms_binary_28", estimand = "log(OR)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_azTRUE, model = "brms_binary_28", estimand = "log(OR)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_favipiravirTRUE, model = "brms_binary_28", estimand = "log(OR)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ) )
Survival time models
prior_survival_time <- prior(normal(0, 5), class = "b") + prior(normal(0, 5), class = "sd")
Fewer adjustments
Those models are described by the brms formula with family brmsfamily("cox", "log"):
form_survival_time_small <- bf( day | cens(censored) ~ sex + age_norm + ever_hcq + ever_az + ever_favipiravir + (1 | hospital_id), family = brmsfamily("cox", "log") ) print_formula(form_survival_time_small)
fit_survival_time_small <- brm( formula = form_survival_time_small, data = data$surv_data_last_day, prior = prior_survival_time, inits = 0, control = default_control, file = paste0(cache_dir, "/fit_survival_time_small") ) fit_survival_time_small <- add_criterion(fit_survival_time_small, "loo")
Here are the fitted coefficients:
summary(fit_survival_time_small) invisible(loo(fit_survival_time_small))
draws <- tidy_draws(fit_survival_time_small) hypothesis_res_list[["hcq_az_survival_time_small"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = draws$b_ever_hcqTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm", model_check = "some influential observations" ), bayesian_hypothesis_res_from_draws( draws = draws$b_ever_azTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm", model_check = "some influential observations" ), bayesian_hypothesis_res_from_draws( draws = draws$b_ever_favipiravirTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm", model_check = "some influential observations" ) )
Fewer adjustments - First wave only
The same model as above, but using data from the first wave only. Here is the summary of the fit:
fit_survival_time_small_fw <- brm( formula = form_survival_time_small, data = data$surv_data_last_day %>% filter(first_wave), prior = prior_survival_time, inits = 0, control = default_control, file = paste0(cache_dir, "/fit_survival_time_small_fw") ) fit_survival_time_small_fw <- add_criterion(fit_survival_time_small_fw, "loo")
summary(fit_survival_time_small_fw) loo(fit_survival_time_small_fw)
draws <- tidy_draws(fit_survival_time_small_fw) hypothesis_res_list[["hcq_az_survival_time_small_fw"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = draws$b_ever_hcqTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, first_wave", model_check = "some influential observations" ), bayesian_hypothesis_res_from_draws( draws = draws$b_ever_azTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, first_wave", model_check = "some influential observations" ), bayesian_hypothesis_res_from_draws( draws = draws$b_ever_favipiravirTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, first_wave", model_check = "some influential observations" ) )
More adjustments
For those models the formula changes to:
form_survival_time <- bf( day | cens(censored) ~ sex + age_norm + smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir + (1 | hospital_id), family = brmsfamily("cox", "log") ) print_formula(form_survival_time)
fit_survival_time <- brm( formula = form_survival_time, data = data$surv_data_last_day, prior = prior_survival_time, inits = 0, control = default_control, file = paste0(cache_dir, "/fit_survival_time") ) fit_survival_time <- add_criterion(fit_survival_time, "loo")
Here are the fitted coefficients - there is a divergent transition so we are possibly fitting an overparemtrized model and we will mark it as "Suspicious":
summary(fit_survival_time)
In addition, the LOO diagnostics warn us that there were some very influential observations, indicating possible overfitting.
loo(fit_survival_time)
draws <- tidy_draws(fit_survival_time) hypothesis_res_list[["hcq_az_survival_time"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = draws$b_ever_hcqTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = draws$b_ever_azTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = draws$b_ever_favipiravirTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ) )
More adjustments - First wave only
The same model as above, but fitting only to patients from first wave. Here are the fitted model coefficients:
fit_survival_time_fw <- brm( formula = form_survival_time, data = data$surv_data_last_day %>% filter(first_wave), prior = prior_survival_time, inits = 0, control = default_control, file = paste0(cache_dir, "/fit_survival_time_fw") ) fit_survival_time_fw <- add_criterion(fit_survival_time_fw, "loo")
summary(fit_survival_time_fw)
We do not see divergent transitions, but some overly influential observations (according to loo) remain, the model will be marked as "Suspicious".
loo(fit_survival_time_fw)
draws <- tidy_draws(fit_survival_time_fw) hypothesis_res_list[["hcq_az_survival_time_fw"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = draws$b_ever_hcqTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = draws$b_ever_azTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = draws$b_ever_favipiravirTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave", model_check = "Suspicious" ) )
Binary discharged models
All the models in this section use the bernoulli("logit") family (i.e. logistic regression) and follow the brms formula
form_discharged_base <- bf( discharged ~ sex + age_norm + smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir + (1 | hospital_id), family = bernoulli("logit") ) print_formula(form_discharged_base)
Basic binary model
prior_discharged_base <- prior(normal(0, 5), class = "b") + prior(normal(0, 5), class = "sd") fit_discharged_base <- brm( formula = form_discharged_base, data = data$patient_data, prior = prior_discharged_base, control = default_control, file = paste0(cache_dir, "/fit_discharged_base") ) fit_discharged_base <- add_criterion(fit_discharged_base, "loo")
Here is the summary of the fitted coefficients
summary(fit_discharged_base) invisible(loo(fit_discharged_base)) # ce <- conditional_effects(fit_discharged_base) # plot(ce, ask = FALSE)
And the error_binned check:
pp_check(fit_discharged_base, "error_binned", nsamples = 20)
draws <- tidy_draws(fit_discharged_base) hypothesis_res_list[["hcq_az_discharged_base"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_hcqTRUE, model = "brms_binary", estimand = "log(OR)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_azTRUE, model = "brms_binary", estimand = "log(OR)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_favipiravirTRUE, model = "brms_binary", estimand = "log(OR)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ) )
Basic binary model (7 Days)
fit_discharged_base_7 <- brm( formula = form_discharged_base, data = data$surv_data_7, prior = prior_discharged_base, control = default_control, file = paste0(cache_dir, "/fit_discharged_base_7") ) fit_discharged_base_7 <- add_criterion(fit_discharged_base_7, "loo")
Here is the summary of the fitted coefficients:
summary(fit_discharged_base_7) invisible(loo(fit_discharged_base_7)) # ce <- conditional_effects(fit_discharged_base_7) # plot(ce, ask = FALSE)
And the error_binned check.
pp_check(fit_discharged_base_7, "error_binned", nsamples = 20)
draws <- tidy_draws(fit_discharged_base_7) hypothesis_res_list[["hcq_az_discharged_base_7"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_hcqTRUE, model = "brms_binary_7", estimand = "log(OR)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_azTRUE, model = "brms_binary_7", estimand = "log(OR)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_favipiravirTRUE, model = "brms_binary_7", estimand = "log(OR)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ) )
Basic binary model (28 Days)
fit_discharged_base_28 <- brm( formula = form_discharged_base, data = data$surv_data_28, prior = prior_discharged_base, control = default_control, file = paste0(cache_dir, "/fit_discharged_base_28") ) fit_discharged_base_28 <- add_criterion(fit_discharged_base_28, "loo")
Here are the fitted coefficients:
summary(fit_discharged_base_28) invisible(loo(fit_discharged_base_28)) # ce <- conditional_effects(fit_discharged_base_28) # plot(ce, ask = FALSE)
And the error_binned check:
pp_check(fit_discharged_base_28, "error_binned")
draws <- tidy_draws(fit_discharged_base_28) hypothesis_res_list[["hcq_az_discharged_base_28"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_hcqTRUE, model = "brms_binary_28", estimand = "log(OR)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_azTRUE, model = "brms_binary_28", estimand = "log(OR)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_favipiravirTRUE, model = "brms_binary_28", estimand = "log(OR)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ) )
Discharged time model
The model uses the family brmsfamily("cox", "log") and the formula:
form_discharged_time <- bf( day | cens(censored_discharged) ~ sex + age_norm + smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir + (1 | hospital_id), family = brmsfamily("cox", "log") ) print_formula(form_discharged_time)
prior_discharged_time <- prior(normal(0, 5), class = "b") + prior(normal(0, 5), class = "sd") fit_discharged_time <- brm( formula = form_discharged_time, data = data$surv_data_last_day, prior = prior_discharged_time, inits = 0, control = list(adapt_delta = 0.98), file = paste0(cache_dir, "/fit_discharged_time") ) fit_discharged_time <- add_criterion(fit_discharged_time, "loo")
The fitted coefficients are:
summary(fit_discharged_time)
The loo diagnostics warn us that there were some very influential observations.
loo(fit_discharged_time)
draws <- tidy_draws(fit_discharged_time) hypothesis_res_list[["hcq_az_discharged_time"]] <- rbind( bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_hcqTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "some influential observations" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_azTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "some influential observations" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_ever_favipiravirTRUE, model = "brms_cox", estimand = "log(HR)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "some influential observations" ) )
Categorical models
All the categorical models use the categorical("logit") family and the following formula:
form_outcome_base <- bf( outcome3 ~ sex + age_norm + smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir + (1 | hospital_id), family = categorical("logit") ) print_formula(form_outcome_base)
Basic Categorical model
prior_outcome_base <- prior(normal(0, 5), class = "b", dpar = "muDeath") + prior(normal(0, 5), class = "sd", dpar = "muDeath") + prior(normal(0, 5), class = "b", dpar = "muDischarged") + prior(normal(0, 5), class = "sd", dpar = "muDischarged") fit_outcome_base <- brm( formula = form_outcome_base, data = data$patient_data, prior = prior_outcome_base, control = default_control, file = paste0(cache_dir, "/fit_outcome_base") ) fit_outcome_base <- add_criterion(fit_outcome_base, "loo")
Fitted coefficients below:
summary(fit_outcome_base) invisible(loo(fit_outcome_base)) # ce <- conditional_effects(fit_outcome_base) # plot(ce, ask = FALSE)
draws <- tidy_draws(fit_outcome_base) hypothesis_res_list[["hcq_az_outcome_base"]] <- rbind( # Death bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_hcqTRUE, model = "brms_categorical", estimand = "logit(death)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_azTRUE, model = "brms_categorical", estimand = "logit(death)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_favipiravirTRUE, model = "brms_categorical", estimand = "logit(death)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), # Discharged bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_hcqTRUE, model = "brms_categorical", estimand = "logit(discharged)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_azTRUE, model = "brms_categorical", estimand = "logit(discharged)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "Suspicious" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_favipiravirTRUE, model = "brms_categorical", estimand = "logit(discharged)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "Suspicious" ) )
Basic Categorical model (7 Days)
fit_outcome_base_7 <- brm( formula = form_outcome_base, data = data$surv_data_7, prior = prior_outcome_base, control = list(adapt_delta = 0.98), file = paste0(cache_dir, "/fit_outcome_base_7") ) fit_outcome_base_7 <- add_criterion(fit_outcome_base_7, "loo")
Summary of model fit:
summary(fit_outcome_base_7) invisible(loo(fit_outcome_base_7)) # ce <- conditional_effects(fit_outcome_base_7) # plot(ce, ask = FALSE)
#pp_check(fit_outcome_base_7, type = "bars", nsamples = 200)
draws <- tidy_draws(fit_outcome_base_7) hypothesis_res_list[["hcq_az_outcome_base_7"]] <- rbind( # Death bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_hcqTRUE, model = "brms_categorical_7", estimand = "logit(death)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_azTRUE, model = "brms_categorical_7", estimand = "logit(death)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_favipiravirTRUE, model = "brms_categorical_7", estimand = "logit(death)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), # Discharged bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_hcqTRUE, model = "brms_categorical_7", estimand = "logit(discharged)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_azTRUE, model = "brms_categorical_7", estimand = "logit(discharged)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_favipiravirTRUE, model = "brms_categorical_7", estimand = "logit(discharged)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ) )
Basic Categorical model (7 Days), first wave only
The same model as above, using only data from the first wave. Summary below:
fit_outcome_base_7_fw <- brm( formula = form_outcome_base, data = data$surv_data_7 %>% filter(first_wave), prior = prior_outcome_base, control = default_control, file = paste0(cache_dir, "/fit_outcome_base_7_fw") ) fit_outcome_base_7_fw <- add_criterion(fit_outcome_base_7_fw, "loo")
summary(fit_outcome_base_7_fw) invisible(loo(fit_outcome_base_7_fw)) # ce <- conditional_effects(fit_outcome_base_7_fw) # plot(ce, ask = FALSE)
#pp_check(fit_outcome_base_7_fw, type = "bars", nsamples = 200)
draws <- tidy_draws(fit_outcome_base_7_fw) hypothesis_res_list[["hcq_az_outcome_base_7_fw"]] <- rbind( # Death bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_hcqTRUE, model = "brms_categorical_7", estimand = "logit(death)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_azTRUE, model = "brms_categorical_7", estimand = "logit(death)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_favipiravirTRUE, model = "brms_categorical_7", estimand = "logit(death)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), # Discharged bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_hcqTRUE, model = "brms_categorical_7", estimand = "logit(discharged)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_azTRUE, model = "brms_categorical_7", estimand = "logit(discharged)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_favipiravirTRUE, model = "brms_categorical_7", estimand = "logit(discharged)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ) )
Basic Categorical model (28 Days)
fit_outcome_base_28 <- brm( formula = form_outcome_base, data = data$surv_data_28, prior = prior_outcome_base, control = default_control, file = paste0(cache_dir, "/fit_outcome_base_28") ) fit_outcome_bas_28 <- add_criterion(fit_outcome_base_28, "loo")
Summary of the model is:
summary(fit_outcome_base_28) #pp_check(fit_outcome_base_28) invisible(loo(fit_outcome_base_28))
draws <- tidy_draws(fit_outcome_base_28) hypothesis_res_list[["hcq_az_outcome_base_28"]] <- rbind( # Death bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_hcqTRUE, model = "brms_categorical_28", estimand = "logit(death)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_azTRUE, model = "brms_categorical_28", estimand = "logit(death)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_favipiravirTRUE, model = "brms_categorical_28", estimand = "logit(death)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), # Discharged bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_hcqTRUE, model = "brms_categorical_28", estimand = "logit(discharged)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_azTRUE, model = "brms_categorical_28", estimand = "logit(discharged)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_favipiravirTRUE, model = "brms_categorical_28", estimand = "logit(discharged)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity", model_check = "probably overconfident" ) )
Basic Categorical model (28 Days), first wave only
The same model, but only for patients from the first wave, summary of the fit:
fit_outcome_base_28_fw <- brm( formula = form_outcome_base, data = data$surv_data_28, prior = prior_outcome_base, control = default_control, file = paste0(cache_dir, "/fit_outcome_base_28_fw") ) fit_outcome_bas_28 <- add_criterion(fit_outcome_base_28_fw, "loo")
summary(fit_outcome_base_28_fw) #pp_check(fit_outcome_base_28_fw) invisible(loo(fit_outcome_base_28_fw))
draws <- tidy_draws(fit_outcome_base_28_fw) hypothesis_res_list[["hcq_az_outcome_base_28_fw"]] <- rbind( # Death bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_hcqTRUE, model = "brms_categorical_28", estimand = "logit(death)", hypothesis = hypotheses$hcq_death, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_azTRUE, model = "brms_categorical_28", estimand = "logit(death)", hypothesis = hypotheses$az_death, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = draws$b_muDeath_ever_favipiravirTRUE, model = "brms_categorical_28", estimand = "logit(death)", hypothesis = hypotheses$favipiravir_death, adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), # Discharged bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_hcqTRUE, model = "brms_categorical_28", estimand = "logit(discharged)", hypothesis = hypotheses$hcq_hospital, adjusted = "az, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_azTRUE, model = "brms_categorical_28", estimand = "logit(discharged)", hypothesis = hypotheses$az_hospital, adjusted = "hcq, favipiravir, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ), bayesian_hypothesis_res_from_draws( draws = -draws$b_muDischarged_ever_favipiravirTRUE, model = "brms_categorical_28", estimand = "logit(discharged)", hypothesis = hypotheses$favipiravir_hospital, adjusted = "hcq, az, sex, age_norm, smoking, obesity, first_wave", model_check = "probably overconfident" ) )
Worst Breathing {.tabset}
Not reported in the main manuscript but we also experimented with ordinal (stopping ratio) models for the worst breathing level experienced by a patient.
form_worst_breathing_base <- bf( worst_breathing ~ sex + age_norm + smoking + mo(obesity) + ever_hcq + ever_az + ever_favipiravir + (1 | hospital_id), family = sratio("logit") ) prior_worst_breathing_base <- prior(normal(0, 5), class = "b") + prior(normal(0, 5), class = "sd") fit_worst_breathing_base <- brm( formula = form_worst_breathing_base, data = data$patient_data, prior = prior_worst_breathing_base, control = default_control, file = paste0(cache_dir, "/fit_worst_breathing_base") ) fit_worst_breathing_base <- add_criterion(fit_worst_breathing_base, "loo")
summary(fit_worst_breathing_base) invisible(loo(fit_worst_breathing_base)) # ce <- conditional_effects(fit_worst_breathing_base, categorical = TRUE) # plot(ce, ask = FALSE)
hypothesis_res_all <- do.call(rbind, hypothesis_res_list) write_csv(hypothesis_res_all, path = here::here("manuscript/analysis_brms.csv"))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.