Chapter 02: From lm/glm to mixed models

Chapter 01 fit a model without saying much about what it means. This chapter fills that gap: what partial pooling buys you, how lmebayes parameterizes a mixed model, and why the two model restrictions (one grouping factor, uncorrelated terms) exist.

There is no sampling in this chapter.

Three ways to handle groups

Take reaction time measured repeatedly on each of \(J\) subjects. You have three options.

Complete pooling. lm(Reaction ~ Days_c). One intercept, one slope, everyone the same. Simple, and wrong whenever subjects genuinely differ: the residuals are correlated within subject, so standard errors are too small.

No pooling. lm(Reaction ~ Days_c * Subject), or 18 separate regressions. Each subject gets their own coefficients, estimated from their own 10 observations and nothing else. Honest, but noisy — a subject with an unlucky run of measurements gets an extreme slope, and nothing pulls it back.

Partial pooling. The mixed model. Each subject still gets their own coefficients, but those coefficients are themselves modeled as draws from a population distribution. A subject with little or noisy data is shrunk toward the population mean; a subject with lots of clean data is barely moved.

The amount of shrinkage is not a tuning choice. It falls out of the ratio of within-group information to between-group variance:

\[ \text{shrinkage for group } j \;\approx\; \frac{1/\tau^2}{1/\tau^2 + n_j/\sigma^2}. \]

Small groups (small \(n_j\)) and tight population distributions (small \(\tau^2\)) shrink more. That expression is not an analogy: it is literally the convergence rate \(\lambda^*\) of the two-block Gibbs sampler for the balanced random-intercept model, derived in vignette("Chapter-B02", package = "lmebayesCore"). Models that shrink hard are exactly the models whose sampler mixes slowly, because in both cases the group coefficients and the population mean are tightly coupled.

The parameterization

Here is where lmebayes differs from lme4 in a way you will notice in the output.

lme4 writes the model as a fixed part plus a mean-zero deviation:

\[ y_j = X_j\beta + Z_j b_j + \varepsilon_j, \qquad b_j \sim N(0, G). \]

fixef() gives \(\beta\), ranef() gives \(b_j\) (centered at zero), and coef() gives \(\beta + b_j\).

lmebayes writes the same model in centered, two-stage form (Lindley and Smith 1972):

\[ \begin{aligned} \text{Stage 1:}\quad & y_j \mid \beta_j \sim N(D_j\beta_j,\; \sigma^2 I) \\ \text{Stage 2:}\quad & \beta_j \mid \gamma \sim N(\mathcal{W}_j\gamma,\; \Psi) \end{aligned} \]

The group-level coefficient vector \(\beta_j\) is not mean zero. It is the full coefficient vector for group \(j\), and it enters the likelihood directly. The population coefficients \(\gamma\) appear only in stage 2, as the mean of the \(\beta_j\).

Substituting stage 2 into stage 1 recovers lme4’s form exactly:

\[ y_j = D_j(\mathcal{W}_j\gamma + u_j) + \varepsilon_j = \underbrace{(D_j\mathcal{W}_j)}_{X_j}\gamma + \underbrace{D_j}_{Z_j} u_j + \varepsilon_j, \qquad u_j := \beta_j - \mathcal{W}_j\gamma . \]

So the two are the same model, written differently.

Why it matters in practice. ranef() on an lmerb fit returns \(\beta_j\), not \(u_j\). In Chapter 01, subject 308’s intercept was about 339 ms — an actual expected reaction time, not a +40 offset from the population mean. If you are used to reading ranef() output as deviations that should average to zero, this will look wrong until you recognize the scale.

Why the package does it this way. The centered form is what makes the two-block split natural. Given \(\gamma\), the \(\beta_j\) are conditionally independent across groups, so Block 1 is \(J\) small independent draws. Given all \(\beta_j\), Block 2 is a plain regression of the \(\beta\)’s on the group-level designs. Neither block ever forms a matrix bigger than it has to, and the convergence theory of Nygren (2020) applies directly to this two-block structure.

The level-2 design

\(\mathcal{W}_j\) is the piece with no lme4 counterpart. It maps the population vector \(\gamma\) to group \(j\)’s prior mean.

In the simplest case — no group-level covariates — \(\mathcal{W}_j\) is an identity-like map and every group has the same prior mean. That is the sleepstudy model of Chapter 01, where model_setup() reported both level-2 rows as ~ 1.

Add a group-level covariate and it becomes interesting. If schools have a title1 indicator, the level-2 row for the intercept reads ~ 1 + title1, and each school’s prior mean intercept depends on whether it is a Title I school. Shrinkage is then toward a fitted line rather than a single grand mean. Chapter 04 does exactly this.

Crucially, each random coefficient may carry its own level-2 predictors:

(Intercept)     ~ 1 + private_school + title1
distracted_ppvt ~ 1
distracted_a1   ~ 1 + free_reduced_lunch

An interaction between a group-level variable and an individual-level random slope — free_reduced_lunch:distracted_a1 — is cross-level moderation, and it lands in the level-2 row of distracted_a1. model_setup() prints this mapping, and it is worth reading before every fit: it is the clearest statement of what model you actually specified.

The two restrictions

data(sleepstudy, package = "lme4")
dat <- sleepstudy
dat$Days_c <- dat$Days - mean(dat$Days)
model_setup(Reaction ~ Days_c + (1 + Days_c | Subject), data = dat)
#> Error:
#> ! Correlated random effects are not supported (off-diagonal RE covariance). Use '||' for uncorrelated terms, e.g. (1 + x || group) instead of (1 + x | group).

model_setup(Reaction ~ Days_c + (1 | Subject), data = dat)
#> Error:
#> ! Fixed effects must be constant within Subject (level-2), a population mean slope (matching a random slope name), or cross-level RE moderation (moderator:random_slope); not allowed: Days_c

See also

• Chapter 03 — the prior workflow.
• Chapter 04 — group-level covariates and cross-level moderation.
• vignette("Chapter-B00", package = "lmebayesCore") — the same notation, stated formally.
• vignette("Chapter-B02", package = "lmebayesCore") — where the shrinkage ratio reappears as a convergence rate.