Hi, I am fitting a Proportional Odds model in brms (`family = cumulative("logit")`

) with random-intercepts.

@joels found greater uncertainties around predictions with `posterior_predict(..., re_formula = NA)`

than with `re_formula = NULL`

, which we wouldn’t expect.

As per brms documention:

`re_formula`

formula containing group-level effects to be considered in the prediction. If`NULL`

(default), include all group-level effects; if`NA`

, include no group-level effects.

What is the explanation for this phenomenon?

```
pacman::p_load(brms,
patchwork)
# Data
a <- c(rep(0,0), rep(1,3), rep(2,5), rep(3,5), rep(4,25), rep(5,40), rep(6,24))
b <- c(rep(0,2), rep(1,3), rep(2,9), rep(3,17), rep(4,33), rep(5,18), rep(6,18))
x <- c(rep('medical', length(a)), rep('endovascular', length(b)))
y <- c(a, b)
d1 = data.frame(x, y, study = "RESCUE-Japan LIMIT")
a <- c(rep(0,0), rep(1,3), rep(2,9), rep(3,20), rep(4,36), rep(5,32), rep(6,71))
b <- c(rep(0,2), rep(1,9), rep(2,25), rep(3,31), rep(4,27), rep(5,15), rep(6,68))
x <- c(rep('medical', length(a)), rep('endovascular', length(b)))
y <- c(a, b)
d2 = data.frame(x, y, study = "SELECT2")
a <- c(rep(0,0), rep(1,8), rep(2,18), rep(3,49), rep(4,60), rep(5,45), rep(6,45))
b <- c(rep(0,9), rep(1,19), rep(2,41), rep(3,39), rep(4,45), rep(5,27), rep(6,50))
x <- c(rep('medical', length(a)), rep('endovascular', length(b)))
y <- c(a, b)
d3 = data.frame(x, y, study = "ANGEL-ASPECT")
d = rbind(d1, d2, d3)
d$y <-factor(d$y, ordered = T)
d$x<-factor(d$x, levels = c("endovascular", "medical"))
PO_brms <- brms::brm(
family = cumulative("logit"), # PO model
formula = y ~ x + (1 | study),
prior =
prior(normal(0, 1.5), class = "b") +
prior(student_t(4, 0, 2.5), class = "sd"),
iter = 2000,
chain = 4,
cores = parallel::detectCores(),
control = list(adapt_delta = .99),
backend = "cmdstanr",
seed = 123,
data = d
)
p1a = pp_check(PO_brms, group = "x", type = "bars_grouped", ndraws = 1000,
re_formula=NULL) +
labs(title = "re_formula=NULL") + coord_cartesian(ylim = c(0, 220))
p1b = pp_check(PO_brms, group = "x", type = "bars_grouped", ndraws = 1000,
re_formula=NA) +
labs(title = "re_formula=NA") + coord_cartesian(ylim = c(0, 220))
p1a + p1b + patchwork::plot_layout(guides = 'collect')
```