Repeated measures logistic regression in brms producing unusual estimates

Modeling
1

I am trying to conduct a repeated measures bernoulli regression on the following data

any_11toPCUD.RData (10.2 KB)

Two measurement occasions - all11 followed by pCUD of the same Yes/No outcome meetsCrit_any for each person. These are the proportions answering Yes/No at each timepoint.

  time  meetsCrit_any count   tot  perc
  <fct> <fct>         <int> <int> <dbl>
1 all11 No             1306  2054  0.64
2 all11 Yes             748  2054  0.36
3 pCUD  No             1414  2054  0.69
4 pCUD  Yes             640  2054  0.31

Odds at timepoint all11 is 0.36/(1-0.36) = 0.56 and at timepoint pCUD is 0.31/(1-0.31) = 0.45, resulting in a manual odds-ratio of pCUD/all11 of 0.80. Admittedly this manual assumes independence of observations, which is not true. So with these hand calculations I proceed to the repeated measures logistic regression

brm(formula = meetsCrit_any ~ time + (1|record_id),
    data = any_11toPCUD,
    family = bernoulli) -> mod_bayes

mod_bayes

# output
# Regression Coefficients:
#           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept   -12.01      3.60   -21.18    -7.19 1.02      288      463
# timepCUD     -5.27      1.55    -9.33    -3.24 1.01      411      526

Now clearly this log-odds coefficient is ridiculous given the relatively small difference in proportion, and the Rhat and ESS are very poor given 4000 samples. I am perplexed: thought brms would handle such a simple model effortlessly. Perhaps not so simple.

Can anyone give me some advice on what I am doing wrong and why such a simple model should fail so spectacularly?

Iā€™d really appreciate the insights of the forum here. After finding a spookily similar previous post of mine evidently from a time when I could (perfunctorily) manually code in stan, I followed of the responders advice and did some prior predictive checks. I tried a normal(0, 0.5) on the b and a normal(0,1) on the intercept and the diagnostics were all much better. However still a very large effect of log-odds = -1.9, which is yields an odds ratio of 0.15, nowhere near the 0.8 I calculated manually. Can anyone give me some advice?

Postscript

For reference an equivalent nhst analysis in the GLMMadaptive package using marginal coefficients yields much more sensible results. I include code here for completeness

library(GLMMadaptive)
mixed_model(fixed = meetsCrit_any ~ time,
            random = ~1|record_id,
            data = any_11toPCUD,
            family = binomial) -> mod_glmm

marginal_coefs(object = mod_glmm,
               std_errors = T,
               sandwich = T) -> margModTime

margModTime

#output
#             Estimate Std.Err  z-value p-value
# (Intercept)  -0.7327  0.0189 -38.6803 < 1e-04
# timepCUD     -0.2330  0.0299  -7.8025 < 1e-04

When exponentiated, -0.2330 yields an odds ratio of 0.79, very close to what I calculated manually above.

2

Two bernoulli observations per random effect level contain virtually no information about the value of the parameter. With an underlying bernoulli probability of 0.5, the maximum likelihood estimate would be infinite half the time. My guess is that the rest of your problems all stem from this

2 Likes
3

Sounds reasonable when you put it like that.

4

I just wanted to note that your main interest seems to be in the marginal or population-averaged odds-ratio. The output from the brms hierarchical model gives you the conditional estimates. Although these conditional estimates may look off (mainly due to the reasons @jsokolarjsokolar mentioned), we can still obtain the marginal estimates using the brmsmargins
package (continuing your example):

library(brmsmargins)

ames0 ā† brmsmargins(
object = mod_bayes,
at = data.frame(time = c(ā€œall11ā€, ā€œpCUDā€)),
contrasts = cbind(ā€œAME timeā€ = c(-1, 1)),
effects = ā€œintegrateoutREā€,
backtrans = ā€œresponseā€,
k = 100L,
seed = 1234,
CI = 0.95
)
str(ames0, 1)
## List of 4
## $ Posterior : num [1:4000, 1:2] 0.229 0.382 0.341 0.454 0.457 ...
## $ Summary :Classes ā€™data.tableā€™ and ā€™data.frameā€™: 2 obs. of 10 variables:
## ..- attr(*, ".internal.selfref")=<externalptr>
## $ Contrasts : num [1:4000, 1] -0.0288 -0.0755 -0.0298 -0.0281 -0.0889 ...
## ..- attr(*, "dimnames")=List of 2
## $ ContrastSummary:Classes ā€™data.tableā€™ and ā€™data.frameā€™: 1 obs. of 11 variables:
## ..- attr(*, ".internal.selfref")=<externalptr>

ames0$Summary # returns the marginal probabilities
## M Mdn LL UL PercentROPE PercentMID CI CIType
## <num> <num> <num> <num> <num> <num> <num> <char>
## 1: 0.3685718 0.3691495 0.2700247 0.4595115 NA NA 0.95 HDI
## 2: 0.3141489 0.3142302 0.2250771 0.4041698 NA NA 0.95 HDI
## ROPE MID
## <char> <char>
## 1: <NA> <NA>
## 2: <NA> <NA>

ames0$ContrastSummary # difference between marginal probabilities
## M Mdn LL UL PercentROPE PercentMID CI
## <num> <num> <num> <num> <num> <num> <num>
## 1: -0.05442294 -0.05343183 -0.09096544 -0.02001107 NA NA 0.95
## CIType ROPE MID Label
## <char> <char> <char> <char>
## 1: HDI <NA> <NA> AME time

# calculate marginal (population-averaged) odds-ratio

odds ā† function(x) x / (1 - x)
ames0$Posterior |>
  posterior::as_draws() |>
  posterior::mutate_variables(
    odds_1 = odds(ā€¦1),
    odds_2 = odds(ā€¦2),
    odds_ratio = odds_2 / odds_1
  ) |>
  posterior::summarise_draws(posterior::default_summary_measures())
## # A tibble: 5 x 7
## variable mean median sd mad q5 q95
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ...1 0.369 0.369 0.0483 0.0474 0.290 0.446
## 2 ...2 0.314 0.314 0.0468 0.0470 0.238 0.392
## 3 odds_1 0.592 0.584 0.124 0.121 0.404 0.802
## 4 odds_2 0.464 0.458 0.102 0.0999 0.309 0.639
## 5 odds_ratio 0.784 0.785 0.0652 0.0649 0.676 0.888

So the (population-averaged, marginal) odds-ratio calculated from your brms model is 0.784 [0.676 - 0.888], which is in line with your back of the envelope calculation.

The conditional and marginal values differ a lot because the sd for your random intercept is very large (~36) due to most persons rarely changing status: always No or always Yes). The high sd pushes these personsā€™ probabilities to either 0 or 1.

# A tibble: 3 Ɨ 3
  all11 pCUD      n
  <fct> <fct> <int>
1 No    No     1306
2 Yes   No      108
3 Yes   Yes     640
2 Likes
5

Amazing @hansvancalster. I was aware of individual vs population-level marginal coefficients in the world of classical statistics (see my question on cross validated from many years ago), but thought, incorrectly it seems, that Bayesian analyses have an inherently population-averaged interpretation. This is great. Much appreciated.

6

I edited my first post because it contained a mistake in the odds function. The resulting marginal odds-ratio is now fully in line with your calculations.

7

Yes I noticed the + where a - should be. Thank you again this was very useful.