V = ZGZ′ + R: The One Equation Behind Every Mixed Model

Abstract

Three posts of arithmetic collapse into one line of matrix notation: the marginal covariance of a participant's measurements is Z G Z-transpose plus R. This post stops adding arithmetic and draws the map. It names each object in a sentence, substitutes the running simulated example, and does the multiplication in print, showing that a column of four ones times the between-person variance 6.25 times a row of four ones gives 6.25 in every cell, and that adding 3.75 down the diagonal reproduces the compound-symmetric matrix built by hand in Posts 1 to 3. The equation has exactly two mutable blocks, so every covariance decision ahead is a change to Z and G or a change to R. Expressibility is never defensibility.

---

One Number Per Person: Linear Mixed Models from the Ground Up — Post 4 of 8.

Where we are in the series. Post 3 — Conditional or Marginal? Two Questions, One Mixed Model — finished the arithmetic by writing out both covariance matrices for the same simulated participants: independent once the person is known, correlated at 0.625 once the person is hidden. This post stops adding arithmetic and starts drawing the map, because everything you have read so far is one special case of a single line of matrix notation. Post 5 — Random Slopes: When Each Person Carries a Trajectory — is the first time we deliberately change something inside it.

Every number in this series comes from simulated teaching data; later posts introduce clearly labelled variant simulations. It is not an empirical finding, it describes no real cohort, and it must not be cited as evidence about skin physiology.

The matrix you have now built three times

Here is the object again — the marginal covariance matrix of the four simulated transepidermal water loss (TEWL) measurements that one participant contributes, in the fixed site order forearm, hand, shin, back:

\[ V=\begin{pmatrix} 10.00 & 6.25 & 6.25 & 6.25\\ 6.25 & 10.00 & 6.25 & 6.25\\ 6.25 & 6.25 & 10.00 & 6.25\\ 6.25 & 6.25 & 6.25 & 10.00 \end{pmatrix} \]

You have arrived at this matrix by three different roads. In Post 1 you derived a single cell, $\operatorname{Cov}(Y_{ij},Y_{ik}) = \operatorname{Var}(b_i) = \tau^2 = 6.25$, by cancelling everything that was not shared between two rows of the same person. In Post 2 you filled in the remaining fifteen cells by the observation that no site is special, which delivered compound symmetry and an intraclass correlation of $6.25/10.00 = 0.625$. In Post 3 you rebuilt the same numbers a third time from the law of total variance and the law of total covariance, and watched the off-diagonal appear the moment the person's offset was hidden again.

When one object arrives by three independent routes, it is usually a special case of something more general. It is. Every linear mixed model ever fitted — one random intercept or five random effects, four body sites or forty visits, balanced or ragged — has a marginal covariance matrix of the form $Z_i G Z_i' + R_i$, and the matrix printed above is what that expression returns when the model contains a single shared offset together with independent, equal-variance residuals. One qualification belongs here rather than in a footnote: writing that covariance as one block per person, as this post does throughout, additionally assumes that the rows partition into independent participant clusters. The general form itself is untouched by crossed random effects; the per-person block $V_i$ is not.

The practical payoff is not elegance. It is that the equation has exactly two places where a covariance decision can be lodged, so a reader who owns it can classify every mixed-model covariance option they will ever meet — random slopes, heterogeneous residual variances, autoregressive errors, unstructured covariance — by pointing at one of two blocks.

The objects in the equation, one sentence each

Before the hard part: matrix notation looks intimidating mainly because it is shorthand, not because it is new. One person's four readings are simply stacked into a column, and each object named below is a label for one part of what produced that column: what the population does, how this person departs from it, and what is left over afterwards. Read the list the way you would read the legend under a figure, rather than as algebra waiting to be solved.

The general linear mixed model, in the two-stage form given by Laird and Ware, writes the whole vector of one person's measurements at once:

\[ Y_i = X_i\beta + Z_i b_i + e_i, \qquad b_i \sim N(0, G), \quad e_i \sim N(0, R_i), \quad b_i \perp e_i \]

Read it left to right, one object at a time.

1. $Y_i$ is the stack of measurements contributed by person $i$ — in the running example a $4\times1$ vector holding the forearm, hand, shin and back readings for that participant.
2. $X_i\beta$ is what the population does: the design matrix of fixed-effect columns multiplied by the coefficients that are common to everybody.
3. $Z_i b_i$ is how this person departs from the population, built from a design matrix $Z_i$ that says which departures are allowed and a vector $b_i$ that says how large they are for this individual.
4. $e_i$ is what is left over after both — the measurement-to-measurement noise that neither the population structure nor the person's own offset explains.
5. $G$ is the covariance matrix of the personal departures, and $R_i$ is the covariance matrix of the leftovers.

The dimensions matter more than they look, because they are what makes the multiplication in the next section work:

Object	Dimensions (running example)	What it holds
$Y_i$	4 × 1	four TEWL values, one per site
$X_i$	4 × p	the fixed-effect columns for this person's four rows
$\beta$	p × 1	the population coefficients
$Z_i$	4 × 1	a column of ones — one shared offset, applied to every row
$b_i$	1 × 1	this person's offset
$G$	1 × 1	$\tau^2$, the variance of those offsets
$e_i$	4 × 1	four residuals
$R_i$	4 × 4	the residual covariance matrix

Two statements, one model

Before the hard part: the two lines below are the same two ward-round positions as the previous post, written in the language of distributions. The first is what you expect from a patient whose baseline you already know. The second is what you expect from the next patient, whose baseline you do not. Everything else in this section is bookkeeping about which of the two your software has printed.

With the objects named, the model makes two distinct probabilistic statements, and the whole of Post 3 lives in the gap between them.

Conditional on the person's own random effect — the sealed envelope opened — the four measurements are normal around a mean that includes that person's offset, with covariance $R_i$ and nothing else:

\[ Y_i \mid b_i \sim N\!\left(X_i\beta + Z_i b_i,\ R_i\right) \]

Averaged over persons — the envelope sealed again, because for the next participant walking into the simulated clinic you do not know the offset — the mean loses the offset (its expectation is zero) and the covariance gains it:

\[ Y_i \sim N\!\left(X_i\beta,\ Z_i G Z_i' + R_i\right) \]

This second line is not a new derivation. It is the integration you performed by hand in Post 3 for $k = 4$, written once for every $Z$, every $G$ and every $R$ at the same time, and it is the form in which Verbeke and Molenberghs present the marginal covariance of the linear mixed model. Note what moved and what did not: the random effect changed the covariance and left the mean alone. That is a property of the identity link, not a universal law, and Post 7 collects on the difference when the outcome turns binary.

Substituting the running example

Before the hard part: the multiplication below looks like linear algebra, but all it really does is copy one number into a table. One shared offset per person means every row that person contributes receives the same number. Ask how two of those rows co-vary, and that same number turns up in every off-diagonal cell — the column of ones is only the instruction "give this to all four rows". Watch the 6.25 spread across the whole table, and then watch the 3.75 land on the diagonal alone.

Now do the substitution in print, because the point of a general equation is worthless until the reader has seen it reproduce a number they already trust.

The simulated study contributes $k = 4$ sites per participant and one shared offset per participant, so $Z_i$ is a $4\times1$ column of ones and $G$ is the $1\times1$ matrix holding $\tau^2 = 6.25$. The residuals are independent with a single pooled variance, so $R_i = \sigma^2 I_4 = 3.75\,I_4$. The left-hand block multiplies out as a $4\times1$ times a scalar times a $1\times4$:

\[ Z_i G Z_i' = \begin{pmatrix}1\\1\\1\\1\end{pmatrix} (6.25) \begin{pmatrix}1&1&1&1\end{pmatrix} = 6.25\begin{pmatrix}1&1&1&1\\1&1&1&1\\1&1&1&1\\1&1&1&1\end{pmatrix} = 6.25\,J_4 \]

Every cell of that product is $6.25$ — including the diagonal — because every cell is the same $1 \times 6.25 \times 1$. Adding the residual block puts $3.75$ on the diagonal only:

\[ V_i = Z_i G Z_i' + R_i = 6.25\,J_4 + 3.75\,I_4 = \begin{pmatrix} 10.00 & 6.25 & 6.25 & 6.25\\ 6.25 & 10.00 & 6.25 & 6.25\\ 6.25 & 6.25 & 10.00 & 6.25\\ 6.25 & 6.25 & 6.25 & 10.00 \end{pmatrix} \]

The matrix from Post 2 has reappeared, digit for digit, without anyone deciding what the correlation should be. This is the sentence worth carrying out of the post: compound symmetry is not a setting the software offered and we accepted; it is what a column of ones plus a single pooled residual variance produces. Both halves of that sentence carry weight, and the series is strict about them elsewhere. The all-ones structure of $J_4$ comes from $Z_i$ being all ones, and that is what fixes every off-diagonal cell at $6.25$; the equal diagonal comes from $R_i$ being a single variance times the identity; and the intraclass correlation of $6.25/10.00 = 0.625$ is the arithmetic consequence of the two together. Keep the column of ones but let the residual variances differ from site to site, and the off-diagonals stay flat at $6.25$ while compound symmetry itself is formally dead — which is exactly the demonstration Post 6 is built on. Change either input and the output changes with it, which is precisely the programme of Posts 5 and 6.

What this means when you read papers: a methods section that says only "a random intercept for participant was included" has already told you the covariance structure it assumed — compound symmetry — even if the phrase itself never appears.

Reading the design matrix

Before the hard part: a design matrix is the data sheet the model actually reads: one row per measurement, one column per thing you asked it to account for. The only distinction that matters here is a clinical one. Some columns change as you walk down a single patient's own rows — which site was measured — while others are printed identically in every one of them, because they are properties of the patient rather than of the measurement. That distinction, and nothing more sophisticated, is what decides which comparisons come cheap and which come expensive.

The right-hand side of the model has a second matrix, and it is the one most readers skip. For the running model — TEWL regressed on site, group, their interaction, and the person-level covariates — $X_i$ has four rows, one per site, and its columns fall into two families:

Column of $X_i$	Number of columns	Varies within a person?
Intercept	1	no
Site indicators (hand, shin, back; forearm is the reference)	3	yes
Group indicator (AD vs control)	1	no
Site × group interaction	3	yes
age	1	no
sex	1	no
phototype (Fitzpatrick, as a factor)	one fewer than the number of levels present	no

Walk down one participant's four rows and only two families move: the site indicators and the interaction terms. Age, sex, phototype and the group label are printed identically in all four rows, because they are properties of the person rather than of the measurement. This is not a bookkeeping detail. It is the seed of the argument in Post 8, where four standard errors come out identical: a contrast that lives entirely in the moving columns can be estimated within a person and never has to pay for the between-person variance, whereas a contrast that needs the constant columns cannot.

The simulated numbers make the point concrete, provided each one is tied to the estimand it belongs to. The mean of one group at one site is an average over the 40 participants of that group, each contributing a single row, and it has a standard error of $\sqrt{10.00/40} = 0.50$. The AD − control contrast at that same site sets two disjoint sets of 40 participants against each other, and it has a standard error of $\sqrt{2\times10.00/40} = 0.71$. The within-person difference between two sites, such as forearm − back, has a standard error of $\sqrt{2\times3.75/40} = 0.43$, smaller than either site mean's own standard error. All three are closed-form values for the balanced, unadjusted cell-means model. Adding the covariate columns of $X_i$ to the fit moves the first two slightly off those figures, whereas the within-person difference is untouched by them, because every column that is constant within a person cancels out of it — the same cancellation, seen from the design matrix, that Post 8 spends a whole post defending.

What this means when you read papers: a within-person comparison and a between-group comparison printed side by side in the same table are not equally precise, and the design matrix is where that difference becomes visible.

The map: two places, and only two

Before the hard part: the table below is the practical use of the whole equation. Whenever you meet an unfamiliar covariance option — in a manual, in a reviewer's comment, in somebody else's methods section — it is one of only two kinds of statement: how people differ from one another, or what the leftover noise within a single person looks like. There is no third place for it to live.

Here is the map the series has been working toward. Every modelling decision available to you changes either the left term or the right term, and each purchase has a price.

You change…	What it buys	What it costs
$Z_i$ and $G$ — add columns to $Z$, enlarge $G$	Person-specific departures beyond a constant offset: a trajectory, a site-specific tendency, a treatment-response term. The marginal variance stops being constant and the correlation stops being the same for every pair.	More variance and covariance parameters, a harder likelihood surface, and estimates that can collapse to the boundary when the data cannot support them. Post 5.
$R_i$ — relax the shape of the residual covariance	Residual variances that differ by unit, residual correlations that decay with distance, or an unstructured matrix that assumes nothing at all.	Parameters spent fast, an identifiability ceiling that arrives sooner than most people expect, and the loss of $\tau^2$ as a separable quantity when the ceiling is hit. Post 6.

The map is a claim about notation, not a licence to model. Any covariance structure you can imagine can be written in this form; that a structure is expressible is no evidence at all that it is defensible for your data. Which $G$ and which $R$ to fit remains a substantive decision, argued from the design that generated the measurements and from the diagnostics of the fitted model, and it should never be argued from which output looks tidier.

Where the software prints $G$ and $R$

Both terms are printed by every standard implementation, just under different headings. The Spec A model, fitted two ways in R:

# Spec A — random intercept, pooled residual variance
library(nlme)
fitA <- lme(tewl ~ site * group + age + sex + phototype,
            random = ~ 1 | id, data = skin, method = "REML")

# equivalent in lme4
library(lme4)
fitA2 <- lmer(tewl ~ site * group + age + sex + phototype + (1 | id), data = skin)

* Stata equivalent of Spec A
* mixed estimates by ML unless reml is requested; lme() and lmer() above are REML
mixed tewl i.site##i.group age i.sex i.phototype || id:, reml

In R, VarCorr() reports the random-effects panel — that is $G$, here the single number $\tau^2$ — together with the residual term of $R$. Read the column label before you read the number: what sits inside $R$ is the residual variance $\sigma^2 = 3.75$, while several printouts show only its square root, the residual standard deviation of approximately $1.94$, and the two are not interchangeable. The function getVarCov() will hand you an assembled matrix, but not by default the one you want here: its default is type = "random.effects", which returns $G$, so the $V_i$ printed above has to be requested explicitly with getVarCov(fitA, individuals = "S01", type = "marginal"). In Stata, the random-effects table and the residual panel of mixed carry the same two objects in the same order, reported as variances unless stddeviations is asked for. Reassembling $V_i$ by hand from those panels once, early in a project, is a small habit that may save a great deal of confusion later, because from then on every line of the output has a known address inside the equation.

What this post does not license

The generality of $V_i = Z_i G Z_i' + R_i$ is easy to over-read, so four boundaries are worth stating explicitly.

First, the equation is a container, not a recommendation: it accommodates compound symmetry exactly as comfortably as it accommodates a structure that contradicts it, and it offers no opinion about which one your measurements deserve.

Second, $G$ and $R$ are not always separately recoverable from the data — a sufficiently free $R$ can absorb everything $G$ was carrying, and the parameter counting that shows exactly when this happens is the business of Post 6.

Third, $b_i$ is not a per-person parameter that the model estimated; the predicted values the software will give you are shrunken predictions of a random quantity, and they should be read as such.

Finally, the tidy separation between "the random effect moves the covariance" and "the fixed effects carry the mean" is a property of the linear model with an identity link, and it does not survive the move to a binary outcome.

Deliberately out of scope

This series makes six omissions on purpose, and each is named once rather than quietly skipped.

1. REML versus ML estimation mechanics — the two differ in how the variance components are estimated and in which likelihoods may be compared, and the practical consequence is stated where it bites, in Post 5, rather than derived here.
2. Denominator degrees of freedom (Satterthwaite, Kenward–Roger) — an approximation problem about testing, not about the covariance structure this series is built around; Post 8 gives it one paragraph as a reporting item.
3. Bayesian mixed models — the same $Z$, $G$ and $R$ appear, but priors and posterior computation are a separate curriculum.
4. Crossed and partially crossed random effects — the general form $V = ZGZ' + R$ extends without modification, but the per-person block $V_i$ does not, because the rows no longer partition into independent participant clusters; the intuition of "one number per person" does not survive either, and the case deserves its own treatment.
5. Generalised and non-linear mixed models beyond the single logistic case used in Post 7 — the conditional-versus-marginal divergence is the only aspect this series needs.
6. Missing-data mechanisms beyond the balance argument in Post 8 — a full treatment would require the missingness taxonomy in its own right.

What to do in your own analysis

First, when you next read a mixed-model output, locate $G$ and $R$ before you read a single coefficient, because the covariance assumption determines every standard error printed underneath it.

Second, write your own $Z_i$ out by hand for one participant — four rows and however many columns your random-effects specification implies — as this is usually enough to reveal whether you have asked for a structure the design can actually support.

Third, whenever you meet an unfamiliar mixed-model option in a manual or a reviewer's comment, ask which of the two blocks it modifies; that single question could reasonably replace most of the confusion around covariance-structure terminology.

Fourth, note in your analysis plan which columns of $X_i$ are constant within a person, since contrasts built from those columns and contrasts built from the within-person columns will not have comparable precision, and it may be worth saying so in the manuscript before a reviewer says it for you.

Finally, state the assumed structure of $G$ and $R$ explicitly in the methods section rather than naming the software command alone — a command name is not a covariance assumption, and readers cannot reconstruct one from the other.

Key takeaways

• Every linear mixed model, however complicated, has the marginal covariance $Z_i G Z_i' + R_i$, and the compound-symmetric matrix built by hand in Posts 1 to 3 is what that expression returns when $Z_i$ is a column of ones and $R_i$ is a single variance times the identity.
• The multiplication is worth doing once in print: a $4\times1$ column of ones times $\tau^2 = 6.25$ times a $1\times4$ row of ones gives $6.25\,J_4$, and adding $3.75\,I_4$ reproduces the 10.00 diagonal and 6.25 off-diagonal exactly.
• Compound symmetry is therefore not an option the software selected on your behalf; it is the arithmetic consequence of one shared offset plus independent, equal-variance residuals.
• The equation has exactly two mutable blocks, so every covariance decision ahead of you is either a change to $Z$ and $G$ (Post 5) or a change to $R$ (Post 6).
• In $X_i$, the site columns move within a person while age, sex, phototype and group do not, and that asymmetry is what will explain the four identical standard errors in Post 8.
• Expressibility is not defensibility: the general form will happily write down a covariance structure that your design cannot support, and only the design and the diagnostics can tell you which one to fit.

Next in the series. In Random Slopes: When Each Person Carries a Trajectory we give $Z$ a second column, watch $G$ grow into a 2×2 matrix, and see the variance stop being constant.

References

1. Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982;38(4):963–974.
2. Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. Springer; 2000.