I have run an experiment looking at abstinence rates among participants in a clinical drug and alcohol trial. There were two groups, those who received the new treatment and those who received placebo. Abstinence during the previous week was measured at 0 (baseline), 4, 8, 12, and 24 weeks. There was a high rate of attrition/missing data so that of the 128 that commenced the trial, only 55 supplied data at follow-up. I want to estimate the group difference in odds of abstinence at each time point, controlling for the other time points. I am especially interested in the difference in odds of abstinence at week 24 (controlling for all the other timepoints).
Here is the data.
# analysing the abstinence data using brms
library(tidyverse)
library(brms)
df <- tibble(id = factor(rep(1:128,each=5)),
time = factor(rep(c(0,4,8,12,24),times=128)),
group = factor(c(rep("placebo",335),rep("treatment",305))),
abs = c(0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, NA, NA, NA, NA, 0, 1, 1, NA, 1, 0, 0, 0, NA, NA, 0, 0, 1, NA, 0, 0, 0, NA, NA, 1, 0, NA, NA, NA, 0, 0, 0, NA, NA, NA, 0, NA, NA, NA, 1, 0, 0, NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, NA, NA, 0, NA, NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, NA, 0, NA, NA, NA, NA, 0, 0, NA, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, NA, 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, 0, NA, NA, NA, NA, 0, 0, 0, NA, NA, 0, 1, 1, 1, NA, 0, 0, 0, 0, NA, 0, 0, 0, NA, NA, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, 0, NA, NA, NA, 0, 0, 1, 1, NA, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, NA, 0, 0, 1, 1, 1, 0, 0, 0, NA, NA, 0, 0, 0, 0, NA, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, NA, 0, 0, NA, NA, NA, 0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, NA, NA, NA, 0, NA, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, 0, 0, NA, 0, NA, NA, NA, NA, 0, 0, 0, 0, 1, 0, 1, NA, NA, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, NA, NA, NA, NA, 0, NA, NA, NA, NA, 0, NA, NA, NA, NA, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, NA, 0, NA, 0, NA, NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, NA, NA, NA, 0, 0, NA, 0, NA, 0, NA, 1, NA, 1, 0, 0, NA, NA, NA, 0, 1, 1, 1, 1, 0, 0, 0, NA, NA, 0, NA, NA, NA, NA, 0, NA, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, 0, 0, NA, 0, 0, NA, 0, NA, 0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, NA, NA, NA, NA, 0, 0, NA, NA, NA, 0, 0, 0, 0, NA, 0, 0, NA, 0, NA, 0, 0, 0, 0, NA, 0, 0, 0, 0, NA, 0, 0, 0, 0, NA, 0, NA, NA, NA, NA, 0, 1, 1, 1, 1, 0, NA, NA, NA, NA, 0, 0, NA, NA, NA, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, NA, NA, NA, 0, NA, NA, NA, NA, 0, 0, 0, 0, 1, 0, 0, NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, NA, NA, NA, 0, 1, 1, 1, 1, 0, 0, 1, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, NA, 0, NA, 0, 0, NA, 0, NA, 0, 0, 0, 0, NA, 0, 0, NA, NA, NA))
I ran a longitudinal bernoulli regression model using brms, with group and time as the fixed effects and id as the random intercept. Here is the code for the model.
m <- brm(abs ~ group*time + (1|id),
data = df,
family = 'bernoulli',
prior = set_prior("normal(0,5)"),
warmup = 1e3,
iter = 1.1e4,
chains = 4,
cores = 4,
seed = 1234)
And here is the output
print(m, pars = "b")
# output
Family: bernoulli
Links: mu = logit
Formula: abs ~ group * time + (1 | id)
Data: df (Number of observations: 440)
Samples: 4 chains, each with iter = 11000; warmup = 1000; thin = 1;
total post-warmup samples = 40000
Group-Level Effects:
~id (Number of levels: 128)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 4.11 0.84 2.72 6.02 1.00 13409 19726
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -10.59 1.91 -14.68 -7.23 1.00 15748 20406
grouptreatment -0.58 1.93 -4.45 3.12 1.00 17525 23863
time4 4.92 1.52 2.17 8.15 1.00 20424 23546
time8 6.59 1.57 3.74 9.90 1.00 19619 23639
time12 6.35 1.58 3.47 9.65 1.00 20103 23562
time24 5.71 1.55 2.91 8.97 1.00 21435 24269
grouptreatment:time4 1.60 1.93 -2.08 5.51 1.00 18437 25637
grouptreatment:time8 1.83 1.94 -1.87 5.71 1.00 17841 24707
grouptreatment:time12 1.87 1.95 -1.85 5.80 1.00 17985 26137
grouptreatment:time24 4.25 2.03 0.40 8.36 1.00 19541 27210
I passed this model into the conditional_effects() function in brms.
# condiitional effects
ce <- conditional_effects(m)
And got the following plots. This plots the effect of group averaged across time
This plots the effects of time averaged across group
What I am confused about is the meaning of the following graph.
I assume this is the log odds of abstinence in each of the ten group x time cells.
But what is being conditioned here? What is each log-odds estimate averaged across? The output from the original object calls the estimates āpopulation-averaged estimatesā yet adding up the coefficients in the original objects to calculate the estimated log-odds of abstinence in the treatment group at week 24 we get -10.59 + -0.58 + 5.71 + 4.25 = -1.21 which is totally different from the modal estimated log-odds for the same group x time cell from conditional_effects.
round(ce$`group:time`[10,8],2)
# output
[1] 0.24
So what is the difference between the population-averaged estimates and the conditional effects estimates?
Secondly does anyone know any way to extract the mcmc estimates of the conditional estimates from this conditional-effects object? Or alternatively how to calculate the conditional effects estimates manually from the original brms object?