Validating predictInterval() against brms

Why compare to brms?

merTools::predictInterval() was written to put prediction intervals on mixed-effect models that are too large or too slow to bootstrap. It does this by simulating from the estimated distribution of the fixed and random effects – fast, but an approximation. Since then, brms has become a standard for fully Bayesian multilevel modeling, and its posterior predictive distribution is about as close to a gold standard as we have. This vignette asks a simple question: how close are predictInterval()’s intervals to the ones brms produces?

We look at three things: a Gaussian random-slopes model, what happens for an entirely unobserved group, and a binomial GLMM. This vignette is precompiled – reproducing it requires the brms package and a working Stan toolchain.

A random-slopes model

We use the bundled hsb data (7,185 students in 160 schools) and model math achievement from student SES, school-mean SES, and a school-varying SES slope. We hold out 20% of students for testing, plus six whole schools to stand in for “new groups” later.

data(hsb)
hsb$schid <- factor(hsb$schid)
new_schools <- sample(levels(hsb$schid), 6)
hsb$.set <- "train"
hsb$.set[hsb$schid %in% new_schools] <- "test_new"
seen <- which(hsb$.set == "train")
hsb$.set[sample(seen, round(0.2 * length(seen)))] <- "test_seen"
train     <- droplevels(hsb[hsb$.set == "train", ])
test_seen <- hsb[hsb$.set == "test_seen" & hsb$schid %in% levels(train$schid), ]
test_new  <- hsb[hsb$.set == "test_new", ]

f1 <- mathach ~ ses + meanses + (ses | schid)
m_lme <- lmer(f1, data = train)
brms_fit_seconds <- system.time(
  m_brm <- fit_brm(f1, train, gaussian(), "brms_hsb_meanses"))[["elapsed"]]

We draw a predictive distribution for each held-out student from each method – predictInterval(returnSims = TRUE) and brms::posterior_predict() – and compute everything (point estimates, intervals, coverage) from those draws.

mt_seen <- mt_draws(m_lme, test_seen)
pp_seen <- posterior_predict(m_brm, newdata = test_seen, allow_new_levels = TRUE)

lme_mean <- predict(m_lme, newdata = test_seen)
brm_mean <- colMeans(posterior_epred(m_brm, newdata = test_seen, allow_new_levels = TRUE))
c(correlation = cor(lme_mean, brm_mean),
  mean_abs_diff = mean(abs(lme_mean - brm_mean)),
  response_sd = sd(test_seen$mathach))

#>   correlation mean_abs_diff   response_sd 
#>    0.99984045    0.03670134    6.71745425

ggplot(data.frame(merTools = lme_mean, brms = brm_mean), aes(merTools, brms)) +
  geom_abline(slope = 1, intercept = 0, linetype = 2, color = "grey50") +
  geom_point(alpha = 0.3, size = 0.8, color = "#1B9E77") + coord_equal() +
  labs(x = "lme4 predict()", y = "brms posterior_epred()") +
  theme_minimal(base_size = 12)

b_mt <- qband(mt_seen, 0.90); b_brm <- qband(pp_seen, 0.90)
c(width_merTools = median(b_mt$upr - b_mt$lwr),
  width_brms = median(b_brm$upr - b_brm$lwr),
  cor_lower = cor(b_mt$lwr, b_brm$lwr),
  cor_upper = cor(b_mt$upr, b_brm$upr))

#> width_merTools     width_brms      cor_lower      cor_upper 
#>     20.1789145     20.2344075      0.9901934      0.9899334

endpts <- rbind(data.frame(bound = "lower", merTools = b_mt$lwr, brms = b_brm$lwr),
                data.frame(bound = "upper", merTools = b_mt$upr, brms = b_brm$upr))
ggplot(endpts, aes(merTools, brms, color = bound)) +
  geom_abline(slope = 1, intercept = 0, linetype = 2, color = "grey50") +
  geom_point(alpha = 0.3, size = 0.8) + coord_equal() +
  labs(title = "90% prediction-interval endpoints", x = "merTools", y = "brms") +
  theme_minimal(base_size = 12)

data.frame(nominal = NOMINAL,
           merTools = coverage(mt_seen, test_seen$mathach, NOMINAL),
           brms     = coverage(pp_seen, test_seen$mathach, NOMINAL))

#>   nominal  merTools      brms
#> 1    0.50 0.4569776 0.4598698
#> 2    0.80 0.7895879 0.7932032
#> 3    0.90 0.9168474 0.9168474
#> 4    0.95 0.9660159 0.9703543

t_lme <- system.time(lmer(f1, data = train))[["elapsed"]]
t_mt  <- system.time(mt_draws(m_lme, test_seen))[["elapsed"]]
t_pp  <- system.time(posterior_predict(m_brm, newdata = test_seen,
                                       allow_new_levels = TRUE))[["elapsed"]]
data.frame(  # brms_fit_seconds was captured when the model was fit, above
  step = c("lme4 fit", "predictInterval", "brms fit (compile+sample)", "posterior_predict"),
  seconds = round(c(t_lme, t_mt, brms_fit_seconds, t_pp), 2))

#>                        step seconds
#> 1                  lme4 fit    0.05
#> 2           predictInterval    0.53
#> 3 brms fit (compile+sample)  398.97
#> 4         posterior_predict    1.52

What about an entirely new group?

For a school that was never seen during fitting there is a genuine modeling choice: do you drop the school’s effect and predict from the fixed effects only, or sample a plausible effect for it from the estimated school-level distribution? Both packages support either, and the comparison is only fair when the choices match:

Matched up, the two tools agree on both the width and the coverage of the new-school intervals:

nd <- test_new
pp_pop  <- posterior_predict(m_brm, newdata = nd, re_formula = NA)
pp_draw <- posterior_predict(m_brm, newdata = nd, allow_new_levels = TRUE)
mt_zero <- mt_draws(m_lme, nd, new.levels = "zero")
mt_draw <- mt_draws(m_lme, nd, new.levels = "draw")
w <- function(d) median(qband(d, 0.90)$upr - qband(d, 0.90)$lwr)
data.frame(
  unseen_group   = c("drop effect", "sample effect"),
  merTools_width = c(w(mt_zero), w(mt_draw)),
  brms_width     = c(w(pp_pop),  w(pp_draw)),
  merTools_cov90 = c(coverage(mt_zero, nd$mathach, 0.90), coverage(mt_draw, nd$mathach, 0.90)),
  brms_cov90     = c(coverage(pp_pop,  nd$mathach, 0.90), coverage(pp_draw, nd$mathach, 0.90)))

#>    unseen_group merTools_width brms_width merTools_cov90 brms_cov90
#> 1   drop effect       19.96801   20.02855      0.8708487  0.8745387
#> 2 sample effect       20.67910   20.76924      0.8892989  0.8856089

So the interval width for an unseen group is set by the choice, not the method. Earlier I compared predictInterval()’s default against brms’s default, which is an apples-to-oranges pairing – the apparent “gap” was the two packages making different default assumptions, not a deficiency in either.

How much the choice matters depends on how much group variance is left after the fixed effects. Here meanses already absorbs most of the between-school variation, so dropping vs. sampling the school effect barely changes the width. Remove meanses and the school SD roughly doubles, so the two options diverge:

m0 <- lmer(mathach ~ ses + (ses | schid), data = train)
data.frame(
  model = c("with meanses", "without meanses"),
  school_SD = c(attr(VarCorr(m_lme)$schid, "stddev")[["(Intercept)"]],
                attr(VarCorr(m0)$schid,    "stddev")[["(Intercept)"]]),
  width_zero = c(w(mt_draws(m_lme, nd, new.levels = "zero")),
                 w(mt_draws(m0,    nd, new.levels = "zero"))),
  width_draw = c(w(mt_draws(m_lme, nd, new.levels = "draw")),
                 w(mt_draws(m0,    nd, new.levels = "draw"))))

#>             model school_SD width_zero width_draw
#> 1    with meanses  1.598503   19.96801   20.67910
#> 2 without meanses  2.125770   19.94993   21.22576

The practical guidance: predictInterval()’s default (new.levels = "zero") answers “what is the expected outcome for a typical member of an unseen group,” while new.levels = "draw" answers “what is the outcome for an unseen group drawn from the population” – the same quantity brms returns by default. Pick the one that matches your question.

A generalized linear mixed model

The story carries over to GLMMs. We fit a binomial model for whether grouse had ticks and compare predicted probabilities – predictInterval(type = "probability") against brms::posterior_epred().

data(grouseticks)
grouseticks$TICKS_BIN <- as.integer(grouseticks$TICKS >= 1)
grouseticks$cHEIGHT   <- as.numeric(scale(grouseticks$HEIGHT))
gk <- grouseticks
gk$.set <- ifelse(runif(nrow(gk)) < 0.2, "test", "train")
gk_tr <- droplevels(gk[gk$.set == "train", ])
gk_te <- gk[gk$.set == "test" & gk$BROOD %in% levels(gk_tr$BROOD) &
            gk$LOCATION %in% unique(gk_tr$LOCATION), ]

gf <- TICKS_BIN ~ YEAR + cHEIGHT + (1 | BROOD) + (1 | LOCATION)
g_lme <- glmer(gf, family = binomial, data = gk_tr,
               control = glmerControl(optimizer = "bobyqa"))
g_brm <- fit_brm(gf, gk_tr, bernoulli(), "brms_grouseticks")

mt_pr <- mt_draws(g_lme, gk_te, type = "probability", include.resid.var = FALSE)
if (max(mt_pr) > 1 || min(mt_pr) < 0) mt_pr <- plogis(mt_pr)
brm_pr <- posterior_epred(g_brm, newdata = gk_te, allow_new_levels = TRUE)
mt_p <- apply(mt_pr, 2, median); brm_p <- apply(brm_pr, 2, median)
c(correlation = cor(mt_p, brm_p), mean_abs_diff = mean(abs(mt_p - brm_p)))

#>   correlation mean_abs_diff 
#>    0.99488587    0.02704796

ggplot(data.frame(merTools = mt_p, brms = brm_p), aes(merTools, brms)) +
  geom_abline(slope = 1, intercept = 0, linetype = 2, color = "grey50") +
  geom_point(alpha = 0.4, color = "#7570B3") + coord_equal(xlim = c(0, 1), ylim = c(0, 1)) +
  labs(title = "Predicted P(tick): merTools vs brms",
       x = "predictInterval(type = 'probability')", y = "posterior_epred()") +
  theme_minimal(base_size = 12)

Takeaways

• For observed groups, predictInterval() reproduces brms’s point estimates and prediction intervals almost exactly, with the same out-of-sample coverage, at a tiny fraction of the computational cost.
• For new groups, you choose how to treat the unseen effect: new.levels = "zero" (the default) drops it and predicts from the fixed effects, while new.levels = "draw" samples it from the estimated covariance – matching brms’s re_formula = NA and allow_new_levels = TRUE respectively. Either way the two packages agree; pick the convention that fits your question.
• The same agreement holds for GLMMs on the probability scale.

In short, when bootstrapping or full MCMC is impractical, predictInterval() is a fast, well-calibrated stand-in for the cases it was designed to handle.