Bayesian analysis of MADRS primary outcome data for the DRESS trial
Suresh Muthukumaraswamy
20:05:2021
We specify a Bayesian mixed effects model equivalent to that in lme4 as:
\[\begin{align*} y & \sim \operatorname{Normal}(\mu, \sigma_{t}) \\ \mu & = \beta_{0} + s_j + (\overrightarrow{\beta_{Time}} \times \overrightarrow{X_{Time}}) + (\overrightarrow{\beta_{Group}} \times \overrightarrow{X_{Group}}) + (\overrightarrow{\beta_{Int}} \times \overrightarrow{X_{Int}})\\ \overrightarrow{\beta_{Age}},\overrightarrow{\beta_{Time}},\overrightarrow{\beta_{Int}} & \sim \operatorname{Normal}(0, 10) \\ \beta_{0} & \sim \operatorname{Normal}(30, 5) \\ s_{j} & \sim \operatorname{Normal}(0, \sigma_{s}) \\ \sigma_{s},\sigma_{t} & \sim \operatorname{Half Student}(3,0,mad(y)) \\ \end{align*}\]
Overarrows are used to indicate vectors with dummy coding of categorical fixed effect variables and \(s_j\) corresponds to the random efffects of participant j. Relatively non-committal priors are used. For the overall intercept a normal prior with mean of 30 with sd of 5 gives 95% coverage of the range 20-40. 20 is the minimum entry criteria for the study while a patient with severe MADRS score of > 40 would typically be excluded. For fixed effects, normal priors are used with an sd of 10 which allows for a range of effect sizes of size comparable to that seen in prior studies of scopolamine (e.g. Drevets and Furey (2010)) and a relatively similar local population for which we tested ketamine as an antidepressant (Sumner et al. (2020)). sd priors are specified as default half student t priors as suggested by Gelman (2006) using the median absolute deviation (mad) of y as the scale factor (Bürkner (2017)). For MCMC modelling Stan is acccessed through the brms interface (Bürkner (2017)). For computing Bayes Factors, a ROPE of range |sd/10| is used (Kruschke and Liddell (2018)).
#*******************************************************************************
# Load libraries, data and set some variables ---------------------------------
library(ggplot2)
library(tidyverse)
library(brms)
library(bayesplot)
library(bayestestR)
library(see)
library(cowplot)
library(ggdist)
library(here)
dresswide <- read_csv(here("DRESS_alldata_anon_publish.csv"))
dressMADRS <- pivot_longer(dresswide, names_to = "Time", cols = starts_with("madrs"), values_to = "MADRS") %>%
select(anon_id, Time, drug, MADRS) %>% rename(ID = anon_id) %>% rename(Group = drug) %>%
filter(Time != "madrs_0_total")
dressMADRS$ID <- factor(dressMADRS$ID) #Participant ID
dressMADRS$Group <- factor(dressMADRS$Group) #Scopolamine or Placebo
dressMADRS$Time <- factor(dressMADRS$Time,
labels = c("Base", "3 Hour", "1 Day", "3 Day", "1 Week", "2 Week", "4 Week", "6 Week"))#*******************************************************************************
# Plot the model outputs - examine chains and overlays and print summary /BF ---
plot(model1)
mcmc_dens_overlay(postModel1, regex_pars = "b_") 
mcmc_dens_overlay(postModel1, pars = c("sigma", "sd_ID__Intercept"))
summary(model1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: MADRS ~ 0 + Intercept + Group + Time + Group:Time + (1 | ID)
## Data: dressMADRS (Number of observations: 320)
## Samples: 3 chains, each with iter = 10000; warmup = 500; thin = 1;
## total post-warmup samples = 28500
##
## Group-Level Effects:
## ~ID (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 6.04 0.80 4.68 7.80 1.00 7216 12236
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 27.14 1.82 23.56 30.75 1.00 6105 11297
## Groupscop 1.00 2.37 -3.69 5.63 1.00 6514 11515
## Time3Hour -7.26 1.87 -10.92 -3.62 1.00 13542 18803
## Time1Day -10.50 1.87 -14.17 -6.84 1.00 14237 19526
## Time3Day -10.22 1.86 -13.92 -6.58 1.00 13154 18227
## Time1Week -9.07 1.85 -12.68 -5.41 1.00 13005 19315
## Time2Week -6.75 1.86 -10.43 -3.09 1.00 13168 18837
## Time4Week -6.09 1.87 -9.73 -2.46 1.00 13598 19716
## Time6Week -6.96 1.88 -10.64 -3.27 1.00 14109 19526
## Groupscop:Time3Hour 1.92 2.38 -2.72 6.57 1.00 14019 19330
## Groupscop:Time1Day -3.46 2.39 -8.20 1.16 1.00 14428 19224
## Groupscop:Time3Day -2.22 2.36 -6.85 2.39 1.00 13689 18599
## Groupscop:Time1Week 0.21 2.37 -4.52 4.84 1.00 13591 19199
## Groupscop:Time2Week -0.57 2.40 -5.31 4.08 1.00 13475 19311
## Groupscop:Time4Week -1.05 2.38 -5.73 3.60 1.00 13946 19144
## Groupscop:Time6Week -1.29 2.39 -5.99 3.40 1.00 14382 19196
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 5.81 0.25 5.35 6.34 1.00 37001 21957
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).sd_y = sd(dressMADRS$MADRS, na.rm = TRUE)
bf1 <- bayesfactor_parameters(model1, prior = NULL, null = c(-sd_y/10, sd_y/10))
bf1 %>% mutate(int = effectsize::interpret_bf((1/bf1$BF), include_value = TRUE))## # Bayes Factor (Null-Interval)
##
## Parameter | BF | int
## -----------------------------------------------------------------------------
## Intercept | 3.190e+09 | extreme evidence (BF = 1/3.19e+09) against
## Groupscop | 0.206 | moderate evidence (BF = 4.86) in favour of
## Time3Hour | 135.933 | extreme evidence (BF = 1/135.93) against
## Time1Day | 11409.734 | extreme evidence (BF = 1/1.14e+04) against
## Time3Day | 14370.537 | extreme evidence (BF = 1/1.44e+04) against
## Time1Week | 2701.094 | extreme evidence (BF = 1/2.70e+03) against
## Time2Week | 70.516 | very strong evidence (BF = 1/70.52) against
## Time4Week | 32.351 | very strong evidence (BF = 1/32.35) against
## Time6Week | 97.653 | very strong evidence (BF = 1/97.65) against
## Groupscop.Time3Hour | 0.278 | moderate evidence (BF = 3.60) in favour of
## Groupscop.Time1Day | 0.630 | anecdotal evidence (BF = 1.59) in favour of
## Groupscop.Time3Day | 0.319 | moderate evidence (BF = 3.14) in favour of
## Groupscop.Time1Week | 0.185 | moderate evidence (BF = 5.40) in favour of
## Groupscop.Time2Week | 0.191 | moderate evidence (BF = 5.24) in favour of
## Groupscop.Time4Week | 0.210 | moderate evidence (BF = 4.76) in favour of
## Groupscop.Time6Week | 0.227 | moderate evidence (BF = 4.40) in favour of
##
## * Evidence Against The Null: [-0.892, 0.892]#*******************************************************************************
# Plot posterior densities and HDIs -------------------------------------------
cnames <- c("Group","Group x 3 Hour", "Group x 1 Day", "Group x 3 Day", "Group x 1 x Week", "Group x 2 Week",
"Group x 4 Week", "Group x 6 Week", "3 Hour", "1 Day", "3 Day", "1 Week", "2 Week", "4 Week",
"6 Week")
ccc <- postModel1 %>% select(starts_with(c("b_G", "b_T"))) %>%
tibble() %>%
rename_all(function(.){cnames}) %>%
pivot_longer(names_to = "cat", values_to = "val", cols = everything()) %>%
mutate(cat = factor(cat))
ccc %>% ggplot(aes(x=val, y = 0, group = cat)) +
stat_halfeye(point_interval = mode_hdi, .width = 0.95, fill = "skyblue") +
geom_vline(xintercept = 0, linetype = 1, size = 1, colour = "red")+
ylab("Probability Density") +
xlab(expression(beta)) +
labs(title =paste("Posterior densities of", expression(beta), "estimates with 95% HDI")) +
theme_cowplot() +
scale_y_continuous(expand = expansion(mult = c(0,0.05))) +
theme(strip.background = element_blank()) +
facet_wrap(~cat, scales = "free") 
# Plot prior and posterior densities with ROPE
bf1$Parameter[2:16] <- cnames
pl <- plot(bf1, stack = TRUE)
pl + xlim(c(-15,15)) +
labs(title = 'Bayes factor analysis - priors, posteriors and ROPE') +
theme_cowplot()
#*******************************************************************************
# Make a summary table --------------------------------------------------------
bnames <- cbind(c("Intercept","Group", "3 Hour", "1 Day", "3 Day", "1 Week", "2 Week", "4 Week", "6 Week",
"Group x 3 Hour", "Group x 1 Day", "Group x 3 Day", "Group x 1 Week", "Group x 2 Week",
"Group x 4 Week", "Group x 6 Week"))
s <- summary(model1)
s <- data.frame(s$fixed) %>%
mutate(BF = bf1$BF) %>%
mutate(`Interpretation for null` = effectsize::interpret_bf((1/bf1$BF), include_value = FALSE)) %>%
mutate(Bulk_ESS = NULL) %>%
mutate(Tail_ESS = NULL) %>%
mutate(Est.Error = NULL) %>%
mutate(Rhat = NULL) %>%
mutate(Effect = bnames) %>%
select(Effect,Estimate, l.95..CI, u.95..CI, BF, `Interpretation for null`) %>%
slice(2:16) %>%
print(digits = 5, row.names = FALSE)## Effect Estimate l.95..CI u.95..CI BF
## Group 0.99595 -3.6908 5.6345 2.0576e-01
## 3 Hour -7.26421 -10.9192 -3.6165 1.3593e+02
## 1 Day -10.49908 -14.1664 -6.8351 1.1410e+04
## 3 Day -10.21882 -13.9212 -6.5795 1.4371e+04
## 1 Week -9.06840 -12.6804 -5.4141 2.7011e+03
## 2 Week -6.75261 -10.4306 -3.0864 7.0516e+01
## 4 Week -6.08762 -9.7292 -2.4554 3.2351e+01
## 6 Week -6.96040 -10.6350 -3.2728 9.7653e+01
## Group x 3 Hour 1.92228 -2.7164 6.5664 2.7815e-01
## Group x 1 Day -3.45865 -8.1963 1.1647 6.3018e-01
## Group x 3 Day -2.22364 -6.8506 2.3915 3.1895e-01
## Group x 1 Week 0.20670 -4.5229 4.8378 1.8530e-01
## Group x 2 Week -0.57386 -5.3054 4.0811 1.9068e-01
## Group x 4 Week -1.05211 -5.7269 3.6005 2.1028e-01
## Group x 6 Week -1.28786 -5.9921 3.3972 2.2724e-01
## Interpretation for null
## moderate evidence in favour of
## extreme evidence against
## extreme evidence against
## extreme evidence against
## extreme evidence against
## very strong evidence against
## very strong evidence against
## very strong evidence against
## moderate evidence in favour of
## anecdotal evidence in favour of
## moderate evidence in favour of
## moderate evidence in favour of
## moderate evidence in favour of
## moderate evidence in favour of
## moderate evidence in favour ofwrite.table(s, "Dress_Bayes.txt") References
Bürkner, Paul-Christian. 2017. “Brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal Article. 2017 80 (1): 28. https://doi.org/10.18637/jss.v080.i01.
Drevets, W. C., and M. L. Furey. 2010. “Replication of Scopolamine’s Antidepressant Efficacy in Major Depressive Disorder: A Randomized, Placebo-Controlled Clinical Trial.” Journal Article. Biol Psychiatry 67 (5): 432–8. https://doi.org/10.1016/j.biopsych.2009.11.021.
Gelman, A. 2006. “Prior Distributions for Variance Parameters in Hierarchical Models (Comment on Article by Browne and Draper).” Journal Article. Bayesian Analysis 1 (3): 515–34. https://doi.org/10.1214/06-BA117A.
Kruschke, John K., and Torrin M. Liddell. 2018. “The Bayesian New Statistics: Hypothesis Testing, Estimation, Meta-Analysis, and Power Analysis from a Bayesian Perspective.” Journal Article. Psychonomic Bulletin & Review 25 (1): 178–206. https://doi.org/10.3758/s13423-016-1221-4.
Sumner, R. L., R. McMillan, M. J. Spriggs, D. Campbell, G. Malpas, E. Maxwell, C. Deng, et al. 2020. “Ketamine Enhances Visual Sensory Evoked Potential Long-Term Potentiation in Patients with Major Depressive Disorder.” Journal Article. Biol Psychiatry Cogn Neurosci Neuroimaging 5 (1): 45–55. https://doi.org/10.1016/j.bpsc.2019.07.002.