2
$\begingroup$

I am analysing 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 over the course of the trial 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.

df <- data.frame(id = factor(rep(1:128,each=5)),
                 time = factor(rep(c(0,4,8,12,24),times=128)),
                 gruop = 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))

The proportions of the observed cases who were abstinent over time look like this

enter image description here

Which, on face value, looks like a decent treatment effect, especially at week 24.

Here is the model in brms, a package that runs stan models using lme4-based syntax. The priors for the fixed effects of group and time on log-odds of abstinence are wide, normal distributions centred on zero.

m <- brm(abs ~ treatment*time + (1|id),
         data = df,
         family = 'bernoulli',
         prior = set_prior("normal(0,5)"),
         warmup = 1e3,
         iter = 1.1e4,
         chains = 4, 
         cores = 4,
         seed = 1234)

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.12      0.84     2.73     6.02 1.00    12899    21612
# 
# Population-Level Effects: 
#                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept               -10.60      1.90   -14.67    -7.24 1.00    16235    21747
# grouptreatment           -0.59      1.94    -4.49     3.13 1.00    19024    24606
# time4                     4.93      1.53     2.17     8.15 1.00    21260    25054
# time8                     6.60      1.57     3.76     9.90 1.00    20332    23917
# time12                    6.35      1.59     3.46     9.69 1.00    21334    25006
# time24                    5.71      1.56     2.85     8.98 1.00    22546    25676
# grouptreatment:time4      1.62      1.95    -2.10     5.53 1.00    21214    27577
# grouptreatment:time8      1.85      1.95    -1.90     5.77 1.00    20731    26607
# grouptreatment:time12     1.88      1.97    -1.89     5.87 1.00    20978    28585
# grouptreatment:time24     4.26      2.04     0.37     8.37 1.00    21738    28906
# 
# Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
# is a crude measure of effective sample size, and Rhat is the potential 
# scale reduction factor on split chains (at convergence, Rhat = 1).

Now this same model fails to converge in packages that use frequentist methods, such as geepack, lme4, GLMMadaptive and glmmTMB packages. What I want to know is why does the Bayesian model converge when the frequentist versions fail and return error messages?

Is it because the regularising priors help guide the algorithm which values to spend its time exploring?

Despite its converging, the result of this model trouble me. Especially concerning is the size of the log-odds coefficient for the grouptreatment:time24effect. Exponentiating 4.26 yields an odds ratio of 70.8! This does not seem like a reliable estimate to me, given that the proportion abstinent at this time point in the placebo group is 6/29 = 20.7% and the proportion in the treatment group is 14/26 = 53.8%

So my second question is why are these log-odds coefficients so high?

The output calls them 'population-level' effects but is it possible they in fact subject-specific effects of the sort usually returned by frequentist version of these models (see here?)

Another possibility that occurred to me is that the estimates at week 24 are unstable because there are so few data at this time point (n=55) compared to week 0 (n=128). I have run frequentist repeated-measures logistic regressions before and not encountered convergence problems. So I guess another supplementary question I have is why do the frequentist models fail to converge in the first place?

I guess my problem is I don't run a whole lot of panel models with binary outcomes. Maybe the estimates of the week-24 group effect are fine and the odds of abstinence at week 24 in the treatment group controlling for all other time points are actually seventy times higher than in the placebo group, but my gut tells me they're not. I just need some advice from people with more experience with these sorts of models, frequentist, Bayesian or otherwise.

$\endgroup$
  • $\begingroup$ it's probably the 0 cases at time zero - note the -10 intercept coefficient. you might need to add a prior to the intercept, or add pseudo data (similar to adding 0.5 success and 0.5 failures under Jeffrey prior). $\endgroup$ – probabilityislogic May 25 at 3:39
  • $\begingroup$ Thank you so much @probabilityislogic. What about just ommitting the baseline data from the dataset? In any case you've given me a much-needed direction to try. $\endgroup$ – llewmills May 26 at 2:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.