Chapter 13: Comparison with lme4 and glmmTMB

This chapter is about scope. lmebayes does a narrower set of models than lme4 (Bates et al. 2015) or glmmTMB (Brooks et al. 2017), and does something different with them. Knowing where the line falls is the difference between using it well and fighting it.

Because glmmTMB is only suggested, the comparison is precomputed by data-raw/make_vignette_outputs.R.

Same model, three packages

data(sleepstudy, package = "lme4")
dat <- sleepstudy
dat$Days_c <- dat$Days - mean(dat$Days)
form <- Reaction ~ Days_c + (1 + Days_c || Subject)

m_lmer <- lme4::lmer(form, data = dat, REML = TRUE)
m_tmb  <- glmmTMB::glmmTMB(form, data = dat, family = gaussian())

ps <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01)
set.seed(2024)
fit <- lmerb(form, data = dat,
             pfamily_list     = pfamily_list(ps),
             dispersion_ranef = ps$group.dispersion,
             n = 1000L, progbar = FALSE)

All three agree to two decimal places. With a prior weight of one percent they should: the prior is nearly flat, so the posterior mean lands on the maximum-likelihood estimate.

The agreement is worth stating precisely, because it is easy to over-interpret. It shows that the Bayesian machinery reproduces the right answer in the case where a right answer is agreed on. It does not show the three approaches are interchangeable.

What is actually different

lmer() and glmmTMB() report a point estimate and a standard error derived from the curvature of the likelihood at the optimum. That standard error is asymptotic: it is exact as the number of groups goes to infinity, and approximate otherwise. With 18 subjects, “otherwise” is where you are.

lmebayes reports draws from the posterior. Uncertainty comes from the spread of those draws, with no asymptotic step:

#> Call:
#>   lmerb(formula = form, data = dat, pfamily_list = pfamily_list(ps), 
#>     dispersion_ranef = ps$group.dispersion, n = 1000L, progbar = FALSE)
#> Bayesian linear mixed model  [1000 draws, exact iid]
#> Formula: Reaction ~ Days_c + (1 + Days_c || Subject) 
#> Random effects (variance components fixed at lmer estimates):
#>  Groups    Name        Std.Dev.
#>  Subject   (Intercept) 37.533  
#>  Subject.1 Days_c       5.922  
#>  Residual              25.592  
#> Number of obs: 180,  groups: Subject, 18
#> --- Posterior means (ICM exact, under fixed variance components) ---
#>   RE component  parameter      fixef.mode   fixef.means    draws SD
#>   ------------  -----------  ------------  ------------  ----------
#>   (Intercept)   (Intercept)      298.5079      298.6709      8.7483
#>   Days_c        (Intercept)       10.3748       10.4519      1.4520

For this model the two happen to agree closely, because the posterior is Gaussian and the sample is large enough. Where they diverge is small numbers of groups, variance components near zero, and non-Gaussian likelihoods — precisely the situations where the asymptotic approximation is weakest.

The second difference is the one that has no counterpart at all. When the header says exact iid, the 1000 draws are independent samples from the exact posterior. When it says two-block Gibbs, each draw is certified within tv_tol in total variation. Neither lme4 nor glmmTMB makes a distributional claim of that kind, because they are not sampling. Compared against other Bayesian packages the claim is still unusual: Stan and JAGS give correlated draws with diagnostics that are heuristic (Gelman and Rubin (1992)), not certificates.

The feature matrix

The three “No” rows in the middle are the cost of the two-block structure. Crossed and nested factors break the conditional independence that makes Block 1 factorize; correlated random effects break the coefficient-wise factorization of Block 2; and the fixed-effect rule follows from the centered parameterization (Chapter 02).

Variance components

lmer() gives one number per component. In Chapter 07 those same components got full posteriors, with credible intervals that were noticeably wide — 18 subjects is little information about how subjects vary. A REML point estimate does not carry that uncertainty forward into the fixed-effect standard errors; the Bayesian fit does, which is the main practical reason to prefer it here.

Choosing

Use lme4 for exploration, for crossed or nested designs, and for any model with correlated random effects. It is fast, well understood, and the reference implementation.

Use glmmTMB for zero-inflation, for families outside the log-concave class, and for complex dispersion models.

Use lmebayes when the model fits its class — one grouping factor, uncorrelated random effects, a log-concave family — and you want exact finite-sample uncertainty rather than an asymptotic approximation. It is strongest exactly where the alternatives are weakest: few groups, small groups, variance components near the boundary.

A reasonable workflow uses both. Fit in lme4 first: it is fast, it tells you whether the random-effect correlation you are about to discard is large, and its estimates are what Prior_Setup_GLMM() calibrates against anyway. Then refit in lmebayes for inference.

m <- lme4::lmer(Reaction ~ Days_c + (1 + Days_c | Subject), data = dat)
print(lme4::VarCorr(m), comp = "Std.Dev.")   # is the correlation large?

If the correlation is near zero, the || restriction costs nothing. If it is large, that is a substantive finding and this package cannot represent it.

Interface differences to remember

## lme4
m <- lmer(y ~ x + (1 + x | g), data = d)
fixef(m); ranef(m); coef(m)

## lmebayes
ps  <- Prior_Setup_GLMM(y ~ x + (1 + x || g), data = d, pop.pwt = 0.01)
fit <- lmerb(y ~ x + (1 + x || g), data = d,
             pfamily_list     = pfamily_list(ps),
             dispersion_ranef = ps$group.dispersion,
             n = 1000)
fixef(fit, type = "hyper.mean")   # fixef(fit) returns the lmer reference
ranef(fit)                        # full coefficients, not deviations

Three traps, all covered earlier: ranef() returns full coefficients rather than mean-zero deviations (Chapter 02); fixef() defaults to the reference lmer fit rather than the posterior (Chapter 05); and the prior is a required argument, not an option (Chapter 03).

See also

• Chapter 02 — why the restrictions exist.
• Chapter 07 — variance components with real uncertainty.
• vignette("Chapter-B02", package = "lmebayesCore") — the convergence theory behind the certified-draw claim.