Interpretation of conditional effects in brms

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
group

This plots the effects of time averaged across group
time

What I am confused about is the meaning of the following graph.
groupTime

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?

Hello!

Iā€™m no brms expert, but from quickly glancing over your post I think this bit is where your got something wrong:

I think conditional_effects defaults to method = "posterior_epred" (the expectations of the posterior predictive distribution), which means that those numbers are on the probability scale. I think you can use method = "posterior_linpred" to get the log-odds scale estimates.

2 Likes

Ah right thank you! That was very helpful. Makes sense given that none of the predictions in from conditional effects are below 0

1 Like