`re_formula = NA` or `NULL`: some inconsistency between modelbased and brms

[AleAnsani]

Hello everyone!
I think I might have found an inconsistency between modelbased's (estimate_means) and brms's conditional_effects when it comes to computing the EMMs with re_formula. In particular, this regards the difference between conditional and marginal effects.

Brief intro:
I found this post by Andrew Heiss very clear when it comes to describing the marginal (re_formula = NULL, default in modelbased) vs conditional (re_formula = NA) effects.

Here's a general definition:

• Conditional effect = re_formula = NA > random offsets set to 0
• Marginal effect = re_formula = NULL: random offsets are incorporated into the estimate

This definition seems to be shared elsewhere. In several blogs, they say that re_formula = NULL should increase the uncertainty in the predictions due to the inclusion of the random effects, whereas with re_formula = NA, the predictions are based only on fixed effects and the estimated error of those parameters. Here, Paul Buerkner is saying this. In the brms manual, it says the same thing:

re_formula: A formula containing group-level effects to be considered in the conditional predictions. If NULL, include all group-level effects; if NA (default), include no group-level effects

However, in modelbased's EMMs, the Credible Intervals with NA are way larger than those with NULL. But shouldn't it be the other way around? If NULL includes all the uncertainty coming from the random part of the model, shouldn't the CIs be larger? I confess I'm a bit confused.

Another user reported a similar doubt, but the replies are still a bit confusing to me. I noticed everyone suggests using NA as default (brms too), but in modelbased NULL seems to be the default, so I wondered why.

I'm attaching here three figures derived by three models (reproducible code below):

• scale(Reaction) ~ Days + (1+Days|Subject), data=sleepstudy

• scale(height) ~ Occasion + (1+Occasion|Subject), data = Oxboys

What I noticed is that the estimates seem more similar with NA (still not identical, though), but they vary enormously with NULL. Compared to NA, modelbased's NULL decreases the variability, whereas brms' conditional_effects increases it. Brms's results seem more logical to me though; cause, again, it makes sense that if you add the variability coming from the RE, the uncertainty should be larger. However, to complicate things even more, the NULL estimates with brms give the impression that nothing is "significant", whereas in fact everything is!

• Finally, this one comes from a more complex model with one random intercept and three random slopes (1 + x + y + z | participant). Here the discrepancies seem even larger.

Here's a reproducible code for the first two figures. I hope it's helpful for you to get to the bottom of this XD
Also, I apologise if I made some mistakes in describing what I found. Maybe it's just me not getting it.

### Marginal / Conditional Effects in B-GLMMs (brms & modelbased) ###

library(pacman)
pacman::p_load(brms, modelbased, lme4)

#### Sleepstudy dataset (slope) ####
data("sleepstudy", package = "lme4")

m <- brm(scale(Reaction) ~ Days + (1+Days|Subject),
         family = gaussian(),
         data = sleepstudy,
         prior = c(
           set_prior("normal(0,1)", class = "Intercept"),
           set_prior("normal(0,1)", class = "b"),
           set_prior("exponential(1)", class = "sd"),
           set_prior("exponential(1)", class = "sigma")
         ),
         cores= 4,
         control = list(adapt_delta = 0.95),
         save_pars = save_pars(all = TRUE), 
         backend = "cmdstanr",
         seed = 1234)

## Marginal effect ##
emm_marginal <- estimate_means(m, by="Days", centrality="mean", ci=.95, ci_method="hdi", re_formula=NULL)

## Conditional effect ##
emm_conditional <- estimate_means(m, by="Days", centrality="mean", ci=.95, ci_method="hdi", re_formula=NA)

ce1 <- conditional_effects(m, effects = "Days", re_formula = NULL)
ce2 <- conditional_effects(m, effects = "Days", re_formula = NA)

p1 <- plot(ce1, plot = FALSE)[[1]]
p2 <- plot(ce2, plot = FALSE)[[1]]

# Merge all graphs
plot(emm_marginal)+labs(title="NULL | modelbased", y="Reaction")+
  plot(emm_conditional)+labs(title="NA | modelbased", y="Reaction")+
  p1 + labs(title="NULL | brms conditional_effects", y="Reaction")+
  p2 + labs(title="NA | brms conditional_effects", y="Reaction")&
  coord_cartesian(ylim=c(-2.5,3))&
  theme_bw()


#### Oxboys dataset (mean difference) ####
data("Oxboys", package = "nlme")
Oxboys <- subset(Oxboys, Occasion <= 4)
Oxboys$Occasion <- factor(Oxboys$Occasion, ordered = FALSE)

m3 <- brm(scale(height) ~ Occasion + (1+Occasion|Subject),
          family = gaussian(),
          data = Oxboys,
          prior = NULL,
          cores= 4,
          control = list(adapt_delta = 0.95),
          save_pars = save_pars(all = TRUE), 
          backend = "cmdstanr",
          seed = 1234)

## Marginal effect ##
emm_marginal3 <- estimate_means(m3, by="Occasion", centrality="mean", ci=.95, ci_method="hdi", re_formula=NULL)

## Conditional effect ##
emm_conditional3 <- estimate_means(m3, by="Occasion", centrality="mean", ci=.95, ci_method="hdi", re_formula=NA)

ce13 <- conditional_effects(m3, effects = "Occasion", re_formula = NULL)
ce23 <- conditional_effects(m3, effects = "Occasion", re_formula = NA)

p13 <- plot(ce13, plot = FALSE)[[1]]
p23 <- plot(ce23, plot = FALSE)[[1]]

# Merge all graphs
plot(emm_marginal3)+labs(title="NULL | modelbased", y="Height (z-score)")+
  plot(emm_conditional3)+labs(title="NA | modelbased", y="Height (z-score)")+
  p13 + labs(title="NULL | brms conditional_effects", y="Height (z-score)")+
  p23 + labs(title="NA | brms conditional_effects", y="Height (z-score)")&
  coord_cartesian(ylim=c(-3,3.5))&
  theme_bw()

[DominiqueMakowski]

@mattansb @strengejacke any thoughts?

[mattansb]

modelbased averages across all levels of the random grouping variable - so each prediction is conditional, the final estimate (and CI) are marginal (see "backend" argument) resulting in lower uncertainty (due to averaging).

brms on the other hand defaults to using the first level of the grouping variable when re_form = NULL.

Neither behavior seems optimal to me 🤷‍♂️

I would expect re_form = NULL to act like brms allow_new_levels = TRUE, sample_new_levels = "uncertainty" to give both the mean prediction and the overall uncertainty (which is closer to a prediction interval really).

[DominiqueMakowski]

I would expect re_form = NULL to act like brms allow_new_levels = TRUE, sample_new_levels = "uncertainty"

But that's a brms suggestion right?

[mattansb]

Right.

[mattansb]

Actually, maybe for brms and here...? What's the point for using all random parameters if you're just going to marginalize them out?

[DominiqueMakowski]

what's the alternative approach / implementation?

[AleAnsani]

Actually, maybe for brms and here...? What's the point for using all random parameters if you're just going to marginalize them out?

Thanks, @mattansb ! That's actually a question I've been asking myself a few times now. I add the random slopes to provide the model with the (un?)famous maximal random structure justified by the experimental design (Barr et al., 2013). Then, as far as I understand, given that my models often have interactions and a lot of parameters, I'd like to base my inferences on marginal contrasts. To do so, I would ideally like to have an estimate that takes into account somehow the variability coming from the random effects (or at least it shouldn't be setting all the random offsets to 0), which is what the NULL promises in brms. Or at least that's what I thought. Now that you mentioned that brms defaults to using the first level of the grouping variable when re_form = NULL, it all makes more sense. In retrospect, now I understand why Paul Buerkner said that we should use NULL just to see the estimates at specific levels of the grouping variable.

I guess then my question for you would be: What would be the best way to compute the contrasts/EMMs so that the CIs reflect the uncertainty coming from the RE but without focusing on just one level of the grouping variable? Is NA the way to go for this? Thanks in advance for any input on this :)

I'm puzzled because, from my naive point of view:

• With NA, you'd be marginalising out all the RE variance (undesirable, and indeed you guys used NULL as default in modelbased)
• With NULL (default in modelbased), you'd incur in what seems to be a paradox (again, to my very naive perspective): you're averaging across all levels of the random grouping variable; so, this averaging, even if it's technically "taking into account the RE" (by averaging across them), produces estimates whose uncertainty is even lower than the one coming from setting all random offsets to 0. In other terms, the re_form that should take into account the RE variability produces more precise estimates than the re_form that explicitly excludes that RE variability. Doesn't this sound conceptually a bit strange?

Again, I apologise if I'm saying nonsense, I'm a noob in this, and I really hope I'm not wasting your time.

Related point
As you seemed to suggest that the marginal estimates of x in y ~ x + (1 + x | id) are identical to the estimates of y ~ x + (1 | id), I gave it a try. When using re_form = NULL, the estimates of the two models are indeed very similar.

scale(Reaction) ~ Days + (1+Days|Subject) [`re_form = NULL`]
Days |  Mean |         95% CI
-----------------------------
0    | -0.84 | [-0.96, -0.71]
5    |  0.09 | [ 0.03,  0.17]
9    |  0.84 | [ 0.72,  0.96]

scale(Reaction) ~ Days + (1|Subject) [`re_form = NULL`]
Days |  Mean |         95% CI
-----------------------------
0    | -0.84 | [-0.97, -0.68]
5    |  0.09 | [ 0.01,  0.18]
9    |  0.84 | [ 0.68,  0.98]

However, when using re_form = NA, they vary:

scale(Reaction) ~ Days + (1+Days|Subject) [`re_form = NA`]
Days |  Mean |         95% CI
-----------------------------
0    | -0.84 | [-1.08, -0.58]
5    |  0.09 | [-0.27,  0.44]
9    |  0.84 | [ 0.28,  1.36]

scale(Reaction) ~ Days + (1|Subject) [`re_form = NA`]
Days |  Mean |         95% CI
-----------------------------
0    | -0.83 | [-1.18, -0.46]
5    |  0.10 | [-0.23,  0.45]
9    |  0.85 | [ 0.48,  1.20]

Thanks for your time on this :)

[mattansb]

To do so, I would ideally like to have an estimate that takes into account somehow the variability coming from the random effects (or at least it shouldn't be setting all the random offsets to 0),

Note that these are two different kinds of estimates and so have different kinds of uncertainties.

Using the maximal random effects structure (Barr 2013) is needed for correct uncertainty about the average effect. So assuming the random effects structure is specified in the model, setting re_form = NA does give proper uncertainty estimates for the average effect.

However if you also want to represent the uncertainty in the individual effects = the heterogeneity of the effect, then setting re_form = NULL and the averaging is not helping you. For this you can do what I suggested above - allow_new_levels = TRUE, sample_new_levels = "uncertainty" - which will give an estimate of the average effect with both sources of uncertainty, or you can look at individual effects.

library(rstanarm)
library(ggplot2)
library(ggdist)
library(tidybayes)
library(posterior)
library(dplyr)

fm1 <- stan_lmer(Reaction ~ Days + (Days | Subject), data = lme4::sleepstudy)

# Mean + uncertainty
slopes1 <- data.frame(Subject = "Doesn't matter", Days = c(0, 1)) |> 
  add_epred_rvars(fm1, re_formula = NA) |> 
  summarize(
    slope = .epred[Days == 1] - .epred[Days == 0]
  )

# Mean + uncertainty + heterogeneity
slopes2 <- data.frame(Subject = "Doesn't matter", Days = c(0, 1)) |> 
  add_epred_rvars(fm1, re_formula = NULL, allow_new_levels = TRUE, sample_new_levels = "uncertainty") |> 
  summarize(
    slope = .epred[Days == 1] - .epred[Days == 0]
  )

# Individual slopes
slopes3 <- expand.grid(Subject = levels(lme4::sleepstudy$Subject), Days = c(0, 1)) |> 
  add_epred_rvars(fm1, re_formula = NULL, allow_new_levels = TRUE, sample_new_levels = "uncertainty") |> 
  summarize(
    slope = .epred[Days == 1] - .epred[Days == 0],

    .by = Subject
  )

ggplot(mapping = aes(xdist = slope)) + 
  stat_slabinterval(data = slopes1, aes(y = "Mean + uncertainty", fill = "Mean + uncertainty")) + 
  stat_slabinterval(data = slopes2, aes(y = "Mean + uncertainty + heterogeneity", fill = "Mean + uncertainty + heterogeneity")) +
  stat_slabinterval(data = slopes3, aes(y = "Individuals", fill = "Individuals"), slab_color = "grey", slab_alpha = 0.5)

Then we can answer questions like:

What is the probability the average slope is positive?

Pr(slopes1$slope > 0)
#> [1] 1

What is the probability that a person sampled at random will have a positive slope? (This is the PPD for the slopes)

Pr(slopes2$slope > 0)
#> [1] 0.93925

What is the probability all individual slopes are positive?

Pr(Reduce("&", slopes3$slope > 0))
#> [1] 0.32775

Etc...

So these different methods can be used to answer different questions.