Mixed-Effects Models

Overview

Mixed-effects models sit at one of those rare junctions where four threads of statistical theory converge on the same algebra. The frequentist GLS-with-known-covariance estimator gives the cleanest answer when the marginal covariance is given. Henderson’s penalized-likelihood derivation produces the same answer without inverting the marginal covariance at all, by treating the random effects as quantities to be solved for jointly with the fixed effects. The REML / Bayesian-marginalization story corrects the maximum-likelihood downward bias on variance components by integrating the fixed effects out under a flat prior. And the full Bayesian fit puts a prior on everything and runs NUTS. When the data is informative enough — which is to say, in the regime where mixed-effects models are most often appropriate — the four routes give numerically equivalent answers. That convergence is the structural reason this model class has been a workhorse of applied statistics for fifty years, and the reason it deserves its own topic in T5.

We develop the four routes in turn on a single running example: six classrooms, $J = 6$, $N = 60$ student-level observations spread unevenly across them. §1 motivates partial pooling geometrically as the third response to grouped data, after complete pooling and no pooling. §2 writes the linear mixed model in matrix form and identifies the marginal covariance $\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}^\top + \mathbf{R}$ that drives generalized least squares. §3 derives Henderson’s mixed-model equations and the closed-form BLUP shrinkage formula $\lambda_j = \tau^2/(\tau^2 + \sigma^2/n_j)$ — the single number that characterizes how mixed-effects partial-pool. §4 derives REML by marginalizing the fixed effects out under a flat improper prior, transparent about the connection to the Gaussian-process marginal-likelihood story. §5 fits the same six-classroom data through both statsmodels.MixedLM (REML) and PyMC (full Bayesian), shows the two views agree numerically, and re-encounters the centered-vs-non-centered story from probabilistic programming in its native habitat. §6 closes with where mixed-effects shine, where they break down, and the connections that motivate the rest of T5.

The reader who has internalized partial pooling from formalStatistics: Hierarchical Bayes and Partial Pooling on the cross-site companion will find this topic the engineering counterpart: same theory, more matrix algebra, code that runs.

1. Three responses to grouped data

1.1 Six classrooms

Six teachers give the same standardized exam to their classes at the end of a term. Class sizes vary — one is a tiny seminar of four students, one is a survey lecture of twenty, the rest fall in between. Alongside each student’s exam score $y$, the school district records a standardized index of preparation hours $x \in [0, 5]$; call it the prep-hours index, with $x = 0$ meaning no out-of-class preparation and $x = 5$ near the upper observed end. We have $N = 60$ student-level observations spread unevenly across $J = 6$ classrooms, with $(n_1, \ldots, n_6) = (4, 6, 8, 10, 12, 20)$.

Panel (a) of Figure 1 shows all sixty points, color-coded by classroom. Some classrooms sit consistently higher (more effective teaching, more prepared students, kinder grading curves — we don’t know, and the data doesn’t tell us); others sit lower. Within each classroom, scores trend upward with prep hours, more or less linearly, more or less consistently. The question is how to summarize the relationship between prep hours and exam score in a way that respects the classroom structure.

There are three natural answers, each tempting in its own way.

1.2 Three classical responses

The first answer is complete pooling: ignore the classroom labels entirely, fit a single line through all sixty points, and report one intercept and one slope. Operationally, this is just OLS on the $(x_i, y_i)$ pairs with no notion of classroom. The fit is clean — one slope, one intercept, narrow standard errors thanks to the full $N = 60$. The cost is also clean: the line is wrong for any classroom whose typical score sits far from the grand mean. A consistently high-scoring classroom sits above the line, a consistently low one sits below, and the line itself is averaging across classroom signal we know is there. We’ve thrown away the structure the data was offering us.

The second answer is no pooling: fit each classroom separately. Concretely, encode classroom identity with six dummy variables and run OLS on the joint design matrix, holding the slope on prep hours shared across classrooms (the prep-hours-to-score relationship is presumably the same across teachers; we’ll relax even this in §2). We get six intercepts, one per classroom, and one shared slope. The fit is now sensitive to classroom — the largest classroom’s intercept is well-estimated, with a tight standard error proportional to $\sigma / \sqrt{n_6} = \sigma/\sqrt{20}$. But the small classrooms? Classroom 1 has four students. A single unusually weak or strong student moves its intercept dramatically; the standard error on that intercept is $\sigma / \sqrt{4} = \sigma/2$, more than twice that of the largest classroom. We’ve gone from one too-aggressive average to six estimates whose reliability varies by a factor of $\sqrt{20/4} \approx 2.2$ across classrooms.

The third answer is partial pooling. Don’t treat the six classroom intercepts as six independent parameters; treat them as six draws from a common population of classroom-level intercepts, and estimate that population’s mean and spread alongside the per-classroom values. The estimate for Classroom $j$‘s intercept becomes a compromise: a weighted average of the no-pooling estimate $\hat\alpha_j^{\mathrm{no\text{-}pool}}$ (what the classroom’s own data says) and the global complete-pooling intercept $\hat\alpha^{\mathrm{pool}}$ (what the entire population says), with the weight depending on how much data the classroom has and how much classrooms vary at the population level. Lots of data plus large between-classroom variation — trust the classroom’s own data, pull little. Little data plus small between-classroom variation — distrust the classroom’s own data, pull hard toward the grand mean.

Panel (c) of Figure 1 shows the partial-pooling fit on the same scatterplot. The lines for the largest classrooms barely move from where panel (b)‘s no-pooling fits placed them. The lines for the smallest classrooms drift visibly toward the grand-mean line that panel (a) drew. The magnitude of the drift is the shrinkage we’ll derive in §3 — for now, it’s a visual fact about the figure.

The geometric punchline: panel (c) sits between panels (a) and (b), and the location depends on group size. Partial pooling is the answer that admits both the within-group signal and the between-group signal exist, and weighs them against each other. The rest of the topic spells out how.

1.3 What makes an effect “random”?

The terminology that names this third response — “random effects,” “mixed-effects models” — has a long history of confusing people, in part because the word random is doing two jobs at once. Andrew Gelman’s 2005 paper “Analysis of variance: why it is more important than ever” surveys five competing definitions in the literature and argues for a pragmatic one we’ll adopt here:

A coefficient is fixed if we model it as a single unknown number to be estimated. A coefficient is random if we model it as a draw from a distribution whose own parameters are estimated.

The slope on prep hours in our setup is fixed: there’s one unknown number $\beta$, and the data estimates it. The classroom intercepts are random: the six $\alpha_j$ values are six draws from a common $\mathcal{N}(\alpha, \tau^2)$, and the data estimates the population mean $\alpha$ and the population spread $\tau$ as well as the realized $\alpha_j$. The terminology is unfortunate — both kinds of coefficients are unknown and both are estimated — but the modeling commitment is real and consequential.

That commitment is what enables shrinkage. Six independent intercepts have no reason to inform one another; they’re separate parameters in a flat parameter space. Six draws from a parent distribution share that parent — the same data that informs the parent’s mean and spread also pulls each child toward its siblings. This is the structural reason small classrooms shrink hard and large classrooms shrink little: a small-$n$ classroom has weak within-group signal that the parent can dominate, while a large-$n$ classroom has strong within-group signal that overpowers the parent’s pull.

Three remarks before we move on. First, fixed and random are modeling choices, not properties of the world. We treat the classroom intercepts as random because we believe the six teachers were sampled from some population of teachers and we want to generalize. If we cared specifically about these six teachers and had no interest in extrapolating to others, fixed-effects with classroom dummies would be the right call. Second, the same coefficient can be modeled differently in different studies — replication studies routinely promote a treatment effect from fixed (one study) to random (across studies), because the across-study generalization is the whole point. Third, the random-effects framing brings in a hyperprior. Frequentist mixed-effects models (Henderson, REML) treat the hyperprior parameters $(\alpha, \tau^2)$ as point estimates; Bayesian mixed-effects models (PyMC, Stan) put a prior on them and integrate. We’ll see in §5 that for the typical case, the two views give numerically similar answers. For now, both views agree on the structural point: classroom-level intercepts are tied together by a common parent.

The rest of the topic makes the pieces precise. §2 writes the linear mixed model in matrix form and identifies the marginal covariance $\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}^\top + \mathbf{R}$ that drives the GLS solution. §3 derives Henderson’s mixed-model equations and the closed-form BLUP shrinkage formula. §4 introduces REML for estimating the variance components $(\tau^2, \sigma^2)$. §5 fits §1’s six-classroom data both ways — statsmodels.MixedLM with REML and a PyMC GLMM — and shows the two views agree.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats

# Six-classroom DGP: J=6 classrooms, N=60 students, varying group sizes
rng = np.random.default_rng(20260429)
J, n_per_group = 6, np.array([4, 6, 8, 10, 12, 20])
N = n_per_group.sum()
classroom = np.repeat(np.arange(J), n_per_group)
alpha_t, beta_t, tau_t, sigma_t = 50.0, 5.0, 5.0, 8.0
a = rng.normal(0.0, tau_t, size=J)
x = rng.uniform(0.0, 5.0, size=N)
y = alpha_t + beta_t * x + a[classroom] + rng.normal(0.0, sigma_t, size=N)
data = pd.DataFrame({"y": y, "x": x, "classroom": classroom})

# Three pooling regimes
pooled = stats.linregress(x, y)                                  # complete pooling
mlm = smf.mixedlm("y ~ x", data, groups="classroom").fit(        # partial pooling
    reml=True, method=["powell"]                                 # Powell — see note
)

The explicit method=["powell"] matters. statsmodels’ default optimizer chain (bfgs / lbfgs / cg) silently returns converged=False on this small-$J$ data with the realized $\hat\tau^2$ near the boundary; Powell converges cleanly to the same REML optimum that §4’s custom implementation finds independently. The shrinkage factors at the REML estimates are $\hat\lambda_j \approx (0.222, 0.298, 0.359, 0.412, 0.451, 0.588)$ across the six classrooms — small classrooms shrink hard, the largest shrinks visibly less.

2. The linear mixed model

§1 gave us three procedures and three pictures. We’re going to fold all three into a single matrix equation, the linear mixed model, and then notice something about its marginal distribution that turns out to drive everything in §3 and §4. The translation costs us nothing geometrically — the same six lines from §1’s panel (c) come out the other side — but it gains us the algebraic apparatus we need to derive Henderson’s equations and to estimate variance components.

2.1 The model in matrix form

Stack the $N = 60$ student-level observations into $\mathbf{y} \in \mathbb{R}^N$, with the row order matching §1’s classroom layout (rows 1–4 are Classroom 1, rows 5–10 are Classroom 2, and so on). The linear mixed model writes $\mathbf{y}$ as a sum of three pieces:

$\mathbf{y} \;=\; \mathbf{X}\boldsymbol{\beta} \;+\; \mathbf{Z}\mathbf{u} \;+\; \boldsymbol{\varepsilon}, \quad\quad (2.1)$

where each piece has a job:

• $\mathbf{X} \in \mathbb{R}^{N \times p}$ is the fixed-effects design matrix. Each column is a covariate that we model with a single shared coefficient. For our six-classroom data, $p = 2$: column one is a column of ones (the intercept), column two is the prep-hours predictor $x_i$.
• $\boldsymbol{\beta} \in \mathbb{R}^p$ is the fixed-effects vector. Here $\boldsymbol{\beta} = (\alpha, \beta)^\top$ — the global intercept and the shared slope on prep hours.
• $\mathbf{Z} \in \mathbb{R}^{N \times q}$ is the random-effects design matrix. Each column corresponds to one entry of the random-effects vector $\mathbf{u}$, and its rows mark which observations are affected by that entry. For a random-intercept-only model with $J = 6$ classrooms, $q = J = 6$ and $\mathbf{Z}$ is the indicator matrix that puts a 1 in row $i$, column $j$ if observation $i$ belongs to Classroom $j$.
• $\mathbf{u} \in \mathbb{R}^q$ is the random-effects vector. Here $\mathbf{u} = (a_1, \ldots, a_6)^\top$ — the six classroom-level intercept deviations.
• $\boldsymbol{\varepsilon} \in \mathbb{R}^N$ is the residual vector — student-level noise.

The distributional assumptions complete the model:

$\mathbf{u} \;\sim\; \mathcal{N}(\mathbf{0}, \mathbf{G}), \qquad \boldsymbol{\varepsilon} \;\sim\; \mathcal{N}(\mathbf{0}, \mathbf{R}), \qquad \mathbf{u} \;\perp\; \boldsymbol{\varepsilon}, \quad\quad (2.2)$

where $\mathbf{G} \in \mathbb{R}^{q \times q}$ is the random-effects covariance and $\mathbf{R} \in \mathbb{R}^{N \times N}$ is the residual covariance. For our running example, both covariances are scaled identities:

$\mathbf{G} \;=\; \tau^2 \mathbf{I}_J, \qquad \mathbf{R} \;=\; \sigma^2 \mathbf{I}_N. \quad\quad (2.3)$

The first equality says the six classroom-level deviations are mutually independent draws from a common $\mathcal{N}(0, \tau^2)$ — six iid pulls from one parent distribution, the operational meaning of “the classrooms are exchangeable draws from a population.” The second says student-level noise is iid Gaussian. We’ll relax both in §6’s remarks but keep them throughout §3–§5; they cover the bulk of standard mixed-effects practice and let us write the cleanest proofs.

The trio (2.1)–(2.3) is the linear mixed model.

Notice what Definition 1 does not commit to: it does not say what $\mathbf{X}$ and $\mathbf{Z}$ look like, or how $\mathbf{G}$ and $\mathbf{R}$ are structured. Different choices encode different scientific assumptions — random intercepts only versus random slopes too, crossed versus nested classrooms, autoregressive errors within a longitudinal subject, and so on. The model in (2.1) is a chassis; the assumptions about $\mathbf{X}, \mathbf{Z}, \mathbf{G}, \mathbf{R}$ are what we drop into it for a particular study.

2.2 Building $\mathbf{X}$ and $\mathbf{Z}$ on the six-classroom data

Concretely, on §1’s dataset:

$\mathbf{X} \;=\; \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_{60} \end{pmatrix} \in \mathbb{R}^{60 \times 2}, \qquad \mathbf{Z} \;=\; \begin{pmatrix} \mathbf{1}_4 & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{1}_6 & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1}_8 & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1}_{10} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1}_{12} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1}_{20} \end{pmatrix} \in \mathbb{R}^{60 \times 6}, \quad\quad (2.4)$

where $\mathbf{1}_k$ denotes a column of $k$ ones and $\mathbf{0}$ is a zero column of the matching height. Reading $\mathbf{Z}$ row by row: row 1 (a Classroom-1 student) is $(1, 0, 0, 0, 0, 0)$, so $\mathbf{Z}\mathbf{u}$ at that row contributes $a_1$ — the deviation of Classroom 1’s intercept from the global $\alpha$. Row 5 (the first Classroom-2 student) is $(0, 1, 0, 0, 0, 0)$, contributing $a_2$. Each row picks out exactly one classroom random effect, which is the matrix encoding of “every student inherits exactly one classroom’s deviation.”

The fixed-effects piece $\mathbf{X}\boldsymbol{\beta} = \alpha \mathbf{1}_N + \beta \mathbf{x}$ is the global linear trend, the same for every student. The random-effects piece $\mathbf{Z}\mathbf{u}$ adds a classroom-specific intercept shift on top. The residual $\boldsymbol{\varepsilon}$ adds student-level noise. Three additive layers of structure, one matrix equation. Panels (a) and (b) of Figure 2 visualize $\mathbf{X}$ and $\mathbf{Z}$ explicitly; the reader should be able to read the block-of-ones pattern in $\mathbf{Z}$ as the matrix encoding of classroom membership.

2.3 The marginal distribution

If we never observed individual students and only saw the data after marginalizing over the random effects, what distribution would we see? This is the marginal distribution of $\mathbf{y}$, and it’s the engine driving everything that follows.

The factor decomposition of (2.5) is worth pausing on. Two students in the same classroom share the same random intercept $a_j$, so they’re positively correlated: their joint covariance is $\tau^2$ from the shared $a_j$, with no additional contribution from the independent residuals. Two students in different classrooms are independent, because their two random intercepts are independent draws from $\mathcal{N}(0, \tau^2)$ and their residuals are also independent. The same student paired with itself contributes $\tau^2 + \sigma^2$ — once for the random intercept and once for the residual.

Plugging the random-intercept-only $\mathbf{G} = \tau^2 \mathbf{I}_J$ and $\mathbf{R} = \sigma^2 \mathbf{I}_N$ into (2.5) makes this concrete. With rows ordered by classroom as in §1, $\mathbf{V}$ becomes block-diagonal, with one $n_j \times n_j$ block per classroom:

$\mathbf{V} \;=\; \mathrm{blkdiag}(\mathbf{V}_1, \ldots, \mathbf{V}_J), \qquad \mathbf{V}_j \;=\; \tau^2 \mathbf{1}_{n_j}\mathbf{1}_{n_j}^\top \;+\; \sigma^2 \mathbf{I}_{n_j}. \quad\quad (2.6)$

Each block has $\tau^2 + \sigma^2$ on its diagonal and $\tau^2$ on every off-diagonal entry — what the longitudinal-data literature calls compound symmetry and what panel (c) of Figure 2 displays as a heatmap. The off-block entries are exactly zero, which is the matrix encoding of “different classrooms are independent.” The visual takeaway: $\mathbf{V}$ is sparse where the modeling commits to independence and dense where it commits to shared structure, and the block pattern is the geometry of the hierarchy itself.

The intra-classroom correlation,

$\rho \;=\; \frac{\mathrm{Cov}(y_{ij}, y_{i'j})}{\sqrt{\mathrm{Var}(y_{ij}) \mathrm{Var}(y_{i'j})}} \;=\; \frac{\tau^2}{\tau^2 + \sigma^2}, \quad\quad (2.7)$

is called the intraclass correlation coefficient (ICC). It’s the fraction of total variance explained by the grouping — $\rho \to 0$ means classrooms barely matter, $\rho \to 1$ means classrooms explain almost everything. For our DGP, $\rho = 25 / (25 + 64) \approx 0.281$, so about 28% of the total variance lives at the classroom level. The ICC is the most-quoted single number in mixed-effects practice, and (2.7) is its full derivation.

2.4 GLS and the chicken-and-egg problem

Theorem 1 hands us a marginal distribution that’s exactly the setup of generalized least squares: a Gaussian linear model with mean $\mathbf{X}\boldsymbol{\beta}$ and a known, non-spherical covariance $\mathbf{V}$. Aitken’s classical theorem says the best linear unbiased estimator (BLUE) of $\boldsymbol{\beta}$ in this setup is the GLS estimator.

We’ll use the GLS covariance in §3 to read off standard errors for the fixed effects, and in §5 to compare frequentist Wald intervals against PyMC posterior intervals.

So far so good — except for one inconvenient fact. We don’t know $\mathbf{V}$. Equation (2.5) gives $\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}^\top + \mathbf{R}$, which depends on $\mathbf{G}$ and $\mathbf{R}$, which depend on the variance components $(\tau^2, \sigma^2)$, which are unknown parameters we have to estimate from the same data we’re using to fit $\boldsymbol{\beta}$. Plugging the true $\mathbf{V}$ into (2.8) gives a clean estimator of $\boldsymbol{\beta}$; plugging the OLS estimate of $\boldsymbol{\beta}$ into a residuals-based estimate of $\mathbf{V}$ would give a biased mess. We need a way to estimate everything jointly.

Two paths through this circularity, and they’re the subjects of §3 and §4 respectively. §3 (Henderson’s path) does not try to invert $\mathbf{V}$ at all; it treats the random effects $\mathbf{u}$ as unknown quantities to be solved for jointly with $\boldsymbol{\beta}$, and writes down the normal equations for the joint penalized log-density. The resulting mixed-model equations solve for $(\hat{\boldsymbol{\beta}}, \hat{\mathbf{u}})$ simultaneously without ever forming $\mathbf{V}^{-1}$ explicitly, and the prediction $\hat{\mathbf{u}}$ has its own name — the best linear unbiased predictor, BLUP — and a clean shrinkage interpretation. §4 (REML’s path) estimates the variance components $(\tau^2, \sigma^2)$ first, then plugs them into (2.8). The naive route — maximum likelihood on the marginal distribution — turns out to be biased downward for the variance components. Restricted maximum likelihood (REML) fixes the bias by marginalizing $\boldsymbol{\beta}$ out before estimating the variance components.

The two paths produce the same $\hat{\boldsymbol{\beta}}$ and the same $\hat{\mathbf{u}}$ when run with the same variance-component estimates, so they’re complementary rather than competitive.

# Build X and Z explicitly, then form V at the truth and solve GLS by hand
import numpy as np
from scipy.linalg import solve

X = np.column_stack([np.ones(N), data["x"].to_numpy()])             # 60 × 2
Z = np.zeros((N, J)); Z[np.arange(N), classroom] = 1.0              # 60 × 6
G_t, R_t = (tau_t ** 2) * np.eye(J), (sigma_t ** 2) * np.eye(N)
V_t = Z @ G_t @ Z.T + R_t                                           # 60 × 60

# GLS at the truth — never form V^{-1}, solve V x = ... instead
beta_gls = solve(X.T @ solve(V_t, X), X.T @ solve(V_t, y), assume_a="pos")
icc_t = tau_t ** 2 / (tau_t ** 2 + sigma_t ** 2)                    # ≈ 0.281

The GLS-with-true-$\mathbf{V}$ fit gives $\hat\alpha \approx 50.99$ and $\hat\beta \approx 4.86$, both within sampling noise of the truths $(50, 5)$. The OLS standard errors (which assume $\mathbf{R} = \sigma^2\mathbf{I}_N$ rather than the block-diagonal $\mathbf{V}$) are systematically too narrow on the intercept by about a factor of $\sqrt{1 + (n_j - 1)\rho}$ averaged over students — the standard “ignore-the-clustering” anticonservatism that mixed-effects practice routinely confronts.

3. Henderson’s mixed-model equations and BLUP

§2 ended on a circularity: GLS gives the cleanest fixed-effects estimator when the marginal covariance $\mathbf{V}$ is known, but $\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}^\top + \mathbf{R}$ depends on variance components we have to estimate from the same data. There are two ways forward, and §3 takes the first one — Henderson’s. Don’t try to invert $\mathbf{V}$ at all. Treat the random effects $\mathbf{u}$ as quantities to be solved for jointly with $\boldsymbol{\beta}$, paying a Gaussian “penalty” for each $u_j$ that pulls it toward zero, and let the penalty’s strength encode the prior $\mathbf{u} \sim \mathcal{N}(\mathbf{0}, \mathbf{G})$. The resulting block linear system is Henderson’s mixed-model equations; its solution simultaneously gives $\hat{\boldsymbol{\beta}}_{\mathrm{GLS}}$ and the best linear unbiased predictor of $\mathbf{u}$, the BLUP.

What we get out of this section is the closed-form shrinkage formula that drove §1’s panel (c) — the BLUP for classroom $j$ is the no-pooling residual mean times a shrinkage factor $\lambda_j = \tau^2 / (\tau^2 + \sigma^2/n_j)$ — together with the structural understanding of why that formula is what falls out of the general theory. We’re still assuming the variance components $(\tau^2, \sigma^2)$ are known throughout §3; §4 fills in their estimation.

3.1 The joint penalized log-density

Henderson’s 1950 / 1975 trick is to write down the joint density of $(\mathbf{y}, \mathbf{u})$ given $\boldsymbol{\beta}$, treat $\mathbf{u}$ as if it were a parameter to be optimized over rather than a latent variable to be marginalized, and maximize the joint log-density over $(\boldsymbol{\beta}, \mathbf{u})$ together. The result is a clean linear system whose left-hand side has dimension $(p + q) \times (p + q)$ — small even when $N$ is enormous — and whose right-hand side reads off the data through $\mathbf{X}^\top \mathbf{R}^{-1} \mathbf{y}$ and $\mathbf{Z}^\top \mathbf{R}^{-1} \mathbf{y}$.

The structure of (3.2) is worth pausing on. The first term is the negative log-likelihood of $\mathbf{y}$ given fixed and random effects — a generalized least-squares residual norm under the residual covariance $\mathbf{R}$. The second term is a quadratic penalty on $\mathbf{u}$: the further $\mathbf{u}$ strays from $\mathbf{0}$ (the random-effects prior mean), the larger the penalty, with the penalty’s strength governed by $\mathbf{G}^{-1}$. The penalty is the algebraic encoding of “we believe a priori that the random effects are small, with the meaning of small set by $\mathbf{G}$.” Larger $\mathbf{G}$ (looser prior) → softer penalty → more freedom for $\mathbf{u}$ to take large values. Smaller $\mathbf{G}$ (tighter prior) → harder penalty → $\mathbf{u}$ is pulled aggressively toward zero. The whole shrinkage story sits inside that one penalty term.

The reader who has seen Bayesian ridge regression will recognize (3.2) immediately: it’s the maximum a posteriori (MAP) objective for $\mathbf{u}$ under a Gaussian prior, with $\boldsymbol{\beta}$ playing the role of an unpenalized nuisance parameter. Henderson didn’t derive his equations as Bayesian MAP — he derived them via a clever marginal-likelihood argument that we’ll see in §4 — but the algebraic form is identical, and the Bayesian reading is what §5 will use to bridge to the PyMC fit.

3.2 The mixed-model equations

Maximizing (3.2) jointly in $(\boldsymbol{\beta}, \mathbf{u})$ — equivalently, minimizing the negative log-density — is a problem in linear algebra. The objective is quadratic in both arguments, so its gradients are linear, and setting both gradients to zero gives a block linear system: Henderson’s mixed-model equations.

The first form of $\hat{\mathbf{u}}$ in (3.4) — $\mathbf{G}\mathbf{Z}^\top\mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}})$ — has a clean probabilistic reading: it is the conditional mean of $\mathbf{u}$ given $\mathbf{y}$ under the joint Gaussian model. Standard multivariate-Gaussian conditioning gives

$\mathbb{E}[\mathbf{u} \mid \mathbf{y}] \;=\; \mathrm{Cov}(\mathbf{u}, \mathbf{y})\,\mathrm{Cov}(\mathbf{y})^{-1}(\mathbf{y} - \mathbb{E}[\mathbf{y}]) \;=\; \mathbf{G}\mathbf{Z}^\top\mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}), \quad\quad (3.9)$

using $\mathrm{Cov}(\mathbf{u}, \mathbf{y}) = \mathrm{Cov}(\mathbf{u}, \mathbf{Z}\mathbf{u} + \boldsymbol{\varepsilon}) = \mathbf{G}\mathbf{Z}^\top$ and $\mathrm{Cov}(\mathbf{y}) = \mathbf{V}$. So the BLUP equals the posterior mean of $\mathbf{u}$ given the data — under the prior $\mathbf{u} \sim \mathcal{N}(\mathbf{0}, \mathbf{G})$ that the model already commits us to. We’ll return to this in §5: the Bayesian PyMC fit doesn’t compute Henderson’s equations, but its posterior mean for the random effects is exactly $\hat{\mathbf{u}}$ (up to the difference between plugged-in REML estimates of $(\mathbf{G}, \mathbf{R})$ and properly-marginalized Bayesian estimates), which is why the two views agree numerically.

The conditional covariance is $\mathrm{Cov}(\mathbf{u} \mid \mathbf{y}) = \mathbf{G} - \mathbf{G}\mathbf{Z}^\top\mathbf{V}^{-1}\mathbf{Z}\mathbf{G}$, a useful identity for the BLUP standard errors below.

A practical note. The MMEs in (3.3) have dimension $(p + q) \times (p + q) = 8 \times 8$ on our running example, while the GLS form (3.4) requires $\mathbf{V}^{-1}$, an $N \times N = 60 \times 60$ matrix. For larger problems the gap widens dramatically — a study with $N = 10^5$ observations across $J = 100$ groups and $p = 5$ fixed effects has MMEs of dimension $105 \times 105$ versus a GLS that nominally inverts a $10^5 \times 10^5$ matrix. The savings is what made Henderson’s formulation practical at all in the pre-computer era and is what makes industrial-scale mixed-effects fits tractable today; modern implementations exploit additional sparsity in $\mathbf{Z}$ (which is highly sparse for grouped data) and in the bottom-right block of (3.3) to bring the cost down further. statsmodels.MixedLM and lme4’s lmer both ultimately solve a sparse version of (3.3).

3.3 Shrinkage in the random-intercept special case

The general form (3.4) is what runs in software. The pedagogical payoff comes from specializing it to the random-intercept-only model, where every quantity in (3.3) collapses to something an undergraduate can compute by hand.

The three identities (3.12)–(3.14) deserve to be read as a unit. The shrinkage factor $\lambda_j$ depends only on the variance ratio $\tau^2 / \sigma^2$ and the group size $n_j$ — nothing else. It is the number that characterizes how mixed-effects models partial-pool. Read it as a competition between two precisions: $1/\tau^2$ is the prior precision of the random effect (how much the population-level prior pulls toward zero), and $n_j/\sigma^2$ is the data precision for group $j$ (how much group $j$‘s own data informs its random effect). The shrinkage factor

$\lambda_j \;=\; \frac{n_j/\sigma^2}{n_j/\sigma^2 + 1/\tau^2} \;=\; \frac{\text{data precision}}{\text{data precision} + \text{prior precision}}$

is the data-precision share of the total. Three regimes deserve names:

• $n_j \to \infty$ at fixed $\tau^2$. Data precision dominates; $\lambda_j \to 1$; $\hat{u}_j \to \bar{r}_j(\hat{\boldsymbol{\beta}})$, the no-pooling estimate. The classroom’s own data overwhelms the prior, and the BLUP equals the no-pooling residual mean.
• $n_j \to 0$ at fixed $\tau^2$. Data precision vanishes; $\lambda_j \to 0$; $\hat{u}_j \to 0$, the prior mean. With no within-classroom data, the BLUP can do nothing but report the prior — every group looks like the population.
• $\tau^2 \to \infty$ at fixed $n_j$. Prior precision vanishes; $\lambda_j \to 1$; $\hat{u}_j \to \bar{r}_j$. An uninformative random-effects prior recovers no pooling. $\tau^2 \to 0$ has the symmetric effect: $\lambda_j \to 0$ and $\hat{u}_j \to 0$, complete pooling.

Equation (3.14) is the conditional-uncertainty companion: when the data dominates ($\lambda_j \to 1$), the conditional variance of $u_j$ collapses to zero; when the prior dominates ($\lambda_j \to 0$), the conditional variance returns to the prior variance $\tau^2$. The two limits are dual — high information content shrinks both the bias and the variance of the BLUP toward zero.

For our six-classroom data with the truth $\tau^2 = 25$, $\sigma^2 = 64$, the shrinkage factors are $\lambda_1 \approx 0.610, \lambda_2 \approx 0.701, \lambda_3 \approx 0.758, \lambda_4 \approx 0.796, \lambda_5 \approx 0.824, \lambda_6 \approx 0.887$, and the conditional variances $(1 - \lambda_j)\tau^2$ run from $9.75$ (Classroom 1) down to $2.83$ (Classroom 6). The smallest classroom has more than three times the BLUP uncertainty of the largest, and shrinks more than three times as hard — both in the same direction, both governed by the same single number $n_j$.

The partial-pooled intercept for Classroom $j$, which is what we plotted in §1’s panel (c), is the global intercept plus the BLUP:

$\hat{\alpha}_j^{\mathrm{partial}} \;=\; \hat{\alpha} + \hat{u}_j \;=\; \hat{\alpha} + \lambda_j \bar{r}_j(\hat{\boldsymbol{\beta}}). \quad\quad (3.15)$

Rewriting $\bar{r}_j$ in terms of the no-pooling intercept $\hat{\alpha}_j^{\mathrm{no\text{-}pool}} = \hat{\alpha} + \bar{r}_j(\hat{\boldsymbol{\beta}})$ — the per-classroom intercept from §1’s panel (b) — gives the equivalent weighted-average form

$\hat{\alpha}_j^{\mathrm{partial}} \;=\; \lambda_j\,\hat{\alpha}_j^{\mathrm{no\text{-}pool}} + (1 - \lambda_j)\,\hat{\alpha}, \quad\quad (3.16)$

which makes panel (c)‘s geometry algebraic: the partial-pooled intercept is a convex combination of the no-pooling estimate (panel b) and the global intercept (panel a’s line), with the convex weight $\lambda_j$ set entirely by group size and the variance ratio.

3.4 BLUE, BLUP, and the predictive frame

The terminology Best Linear Unbiased Predictor — BLUP — deserves a remark. We’ve estimated $\boldsymbol{\beta}$ throughout as a fixed unknown parameter; the GLS estimator $\hat{\boldsymbol{\beta}}$ is the Best Linear Unbiased Estimator (BLUE) under Gauss–Markov conditions, a standard fact about linear models. We’ve also computed $\hat{\mathbf{u}}$, but the random effects $\mathbf{u}$ are not parameters in the same sense — they’re realized values of a random vector with a known prior distribution. The convention in the mixed-effects literature is to call the corresponding estimator a predictor rather than an estimator, on the grammatical grounds that we estimate parameters and predict random variables. The mean-squared-error optimality result is the same in both cases — minimum MSE among linear unbiased rules — but the framing matters for two reasons.

First, the standard error of $\hat{u}_j$ is not the standard error of estimation of a fixed parameter; it is the standard error of prediction of a random quantity, and (3.14) is what governs it. The frequentist prediction interval $\hat{u}_j \pm 1.96 \sqrt{(1 - \lambda_j)\tau^2}$ covers the true realized $u_j$ with the stated frequentist probability under repeated sampling of $(\mathbf{y}, \mathbf{u})$ pairs, not the standard “what is the probability that $u_j$ lies in this interval” — the latter requires a Bayesian posterior, and §5 shows that for this model the posterior interval and the prediction interval numerically coincide.

Second, the unbiasedness of the BLUP is unbiasedness in the predictive sense: $\mathbb{E}[\hat{\mathbf{u}} - \mathbf{u}] = \mathbf{0}$, where the expectation is over both $\mathbf{y}$ and $\mathbf{u}$. Crucially, $\hat{\mathbf{u}}$ is not unbiased for any particular realized $\mathbf{u}$ — it shrinks toward zero, so for a $\mathbf{u}$ that happens to be far from zero, $\hat{\mathbf{u}}$ is systematically biased toward zero in the conditional sense $\mathbb{E}[\hat{\mathbf{u}} \mid \mathbf{u}] \ne \mathbf{u}$. The unbiasedness in the predictive average integrates over the prior, where the expected $\mathbf{u}$ is zero. This is the Stein-estimator family at work — biased point estimates can have lower MSE than unbiased ones — and it is the structural reason mixed-effects models systematically beat both complete pooling and no pooling on out-of-sample prediction.

In modern Bayesian language we’d say nothing fancy at all about any of this: we have a joint Gaussian model, $\mathbb{E}[\mathbf{u} \mid \mathbf{y}]$ is the posterior mean, that’s the MMSE estimator, and the predictive interval is the posterior credible interval. The frequentist BLUP framing was historically prior to the Bayesian framing in this context — Henderson 1950, Patterson and Thompson 1971 — but the two views give the same point estimates and (asymptotically) the same intervals, which is one of the nicer correspondences in applied statistics.

# Henderson MMEs at the truth, then verify against §1's REML fit
G_t, R_t = (tau_t ** 2) * np.eye(J), (sigma_t ** 2) * np.eye(N)
R_inv, G_inv = np.linalg.inv(R_t), np.linalg.inv(G_t)
M = np.block([[X.T @ R_inv @ X,           X.T @ R_inv @ Z],          # (p+q)×(p+q)
              [Z.T @ R_inv @ X,           Z.T @ R_inv @ Z + G_inv]])  # = 8×8 here
rhs = np.concatenate([X.T @ R_inv @ y, Z.T @ R_inv @ y])
sol = solve(M, rhs, assume_a="pos")
beta_mme, u_mme = sol[:2], sol[2:]                                    # MME solve
lam = (tau_t ** 2) / (tau_t ** 2 + sigma_t ** 2 / n_per_group)        # shrinkage

The MME solve gives $\hat\alpha \approx 50.99$, $\hat\beta \approx 4.86$ (matching §2’s GLS-with-true-$\mathbf{V}$ fit to machine precision), and BLUPs $\hat{u}_j$ that agree with mlm.random_effects from §1’s REML fit to within the difference between truth-plug-in and REML-plug-in variance components. The shrinkage factors at the truth are $(0.610, 0.701, 0.758, 0.796, 0.824, 0.887)$ — exactly Proposition 2’s prediction, and the curve panel (a) draws.

4. REML for variance components

§3 assumed the variance components $(\tau^2, \sigma^2)$ were known. They aren’t. Estimating them is where the structural difficulty of mixed-effects fitting actually lives — the fixed effects and random effects come along for free once $(\tau^2, \sigma^2)$ are pinned down via §3’s MMEs, but pinning them down is a job in its own right and the obvious approach (maximum likelihood on the marginal $\mathbf{y} \sim \mathcal{N}(\mathbf{X}\boldsymbol{\beta}, \mathbf{V})$) is biased downward in a way that’s significant for the small-to-moderate group counts that motivate mixed-effects in the first place.

REML — restricted or residual maximum likelihood — is the standard fix. The original Patterson–Thompson 1971 derivation transformed $\mathbf{y}$ to a vector of error contrasts orthogonal to $\mathbf{X}$ and computed the likelihood of those. We’ll take the equivalent and tighter route: marginalize $\boldsymbol{\beta}$ out under a flat improper prior, leaving a residual likelihood that depends only on $(\mathbf{G}, \mathbf{R})$, and maximize it. The structural answer is the same as what you’d get from §1’s textbook degrees-of-freedom correction, generalized to arbitrary mixed models. The mechanics turn out to mirror the Gaussian-process §5 marginal-likelihood-based kernel-hyperparameter learning almost exactly — same Bayesian-marginalization machinery, same Occam’s-razor-flavored objective, different model class.

4.1 Why ML is biased downward

The cleanest illustration is the elementary one: ordinary linear regression. With $\mathbf{y} \sim \mathcal{N}(\mathbf{X}\boldsymbol{\beta}, \sigma^2 \mathbf{I}_N)$, $p$ fixed-effects parameters, no random effects, the maximum-likelihood estimator of $\sigma^2$ is

$\hat{\sigma}^2_{\mathrm{ML}} \;=\; \frac{1}{N}\,\sum_{i=1}^N (y_i - \mathbf{x}_i^\top \hat{\boldsymbol{\beta}})^2 \;=\; \frac{\mathrm{RSS}(\hat{\boldsymbol{\beta}})}{N},$

where $\mathrm{RSS}(\hat{\boldsymbol{\beta}}) = \|\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}\|^2$ is the residual sum of squares at the OLS fit. The textbook unbiased estimator divides by $N - p$ instead of $N$:

$\hat{\sigma}^2_{\mathrm{unbiased}} \;=\; \frac{\mathrm{RSS}(\hat{\boldsymbol{\beta}})}{N - p}, \qquad \mathbb{E}[\hat{\sigma}^2_{\mathrm{unbiased}}] \;=\; \sigma^2. \quad\quad (4.1)$

The factor $N - p$ in the denominator is the residual degrees of freedom: there are $N$ raw observations, but $p$ of them have been used to fit $\hat{\boldsymbol{\beta}}$, leaving $N - p$ “effective” observations to estimate $\sigma^2$. The ML estimator, by dividing by $N$ rather than $N - p$, fails to deduct the cost of estimating $\hat{\boldsymbol{\beta}}$ and consequently underestimates $\sigma^2$ by a factor of $(N - p)/N$. For $N = 100, p = 5$, that’s a 5% downward bias — small but real. For $N = 12, p = 5$, that’s a 58% downward bias — catastrophic.

Mixed-effects fits suffer from the same problem, in a more elaborate form. The marginal log-likelihood we’d plausibly maximize in $(\boldsymbol{\beta}, \mathbf{G}, \mathbf{R})$ jointly is

$\ell_{\mathrm{ML}}(\boldsymbol{\beta}, \mathbf{G}, \mathbf{R}) \;=\; -\tfrac{N}{2}\log 2\pi - \tfrac{1}{2}\log\det \mathbf{V} - \tfrac{1}{2}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}), \quad\quad (4.2)$

with $\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}^\top + \mathbf{R}$. Maximizing in $\boldsymbol{\beta}$ at fixed $(\mathbf{G}, \mathbf{R})$ gives the GLS estimator from §2; maximizing in $(\mathbf{G}, \mathbf{R})$ at the resulting $\hat{\boldsymbol{\beta}}_{\mathrm{GLS}}$ — profiling $\boldsymbol{\beta}$ out — gives the profile log-likelihood in $(\mathbf{G}, \mathbf{R})$:

$\ell_{\mathrm{prof}}(\mathbf{G}, \mathbf{R}) \;=\; -\tfrac{N}{2}\log 2\pi - \tfrac{1}{2}\log\det \mathbf{V} - \tfrac{1}{2}\mathrm{RSS}_\mathbf{V}(\hat{\boldsymbol{\beta}}_{\mathrm{GLS}}), \quad\quad (4.3)$

where $\mathrm{RSS}_\mathbf{V}(\hat{\boldsymbol{\beta}}_{\mathrm{GLS}}) = (\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}_{\mathrm{GLS}})^\top \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}_{\mathrm{GLS}})$ is the GLS-weighted residual sum of squares. Maximizing $\ell_{\mathrm{prof}}$ over $(\mathbf{G}, \mathbf{R})$ gives the ML estimates $(\hat{\mathbf{G}}_{\mathrm{ML}}, \hat{\mathbf{R}}_{\mathrm{ML}})$, and these are biased downward for exactly the same reason (4.1) is: $\ell_{\mathrm{prof}}$ has paid no penalty for the $p$ degrees of freedom consumed by $\hat{\boldsymbol{\beta}}_{\mathrm{GLS}}$.

For our six-classroom example with $p = 2$, the bias is mild — about $(N - p)/N = 58/60 \approx 96.7\%$ of the truth, a 3.3% downward bias. For a study with more covariates or fewer observations the bias scales linearly in $p$. With six covariates and twenty observations, ML would underestimate $\sigma^2$ by 30%. The variance-component estimates inherit that bias, the fixed-effects standard errors (which depend on $\mathbf{V}^{-1}$) inherit it from there, and confidence intervals based on the ML fit are systematically too narrow.

REML’s job is to clean this up.

4.2 The REML derivation by Bayesian marginalization

Patterson and Thompson 1971 derived REML by transforming $\mathbf{y}$ to a vector of error contrasts — linear combinations $\mathbf{a}^\top \mathbf{y}$ with $\mathbf{a}^\top \mathbf{X} = \mathbf{0}$ — and writing down the likelihood of those contrasts. The contrasts are statistics that depend on the data only through residuals from the fixed-effects fit, so their distribution depends on $(\mathbf{G}, \mathbf{R})$ but not on $\boldsymbol{\beta}$, and maximizing their likelihood gives unbiased variance-component estimates. The construction is geometrically clean but algebraically heavy. We’ll take the equivalent route that’s cleaner algebraically and connects directly to the GP §5 story: integrate $\boldsymbol{\beta}$ out under a flat improper prior $p(\boldsymbol{\beta}) \propto 1$ on $\mathbb{R}^p$.