BLUP Derivation for Subject Effects

Why Predict a Subject Effect?

A longitudinal mixed model describes both a population-level mean and subject-to-subject variation around that mean. The random effects are not directly observed: for example, a subject’s deviation from the population intercept is inferred from the subject’s measurements, while allowing for measurement error. The best linear unbiased predictor, or BLUP, formalizes this inference. It combines the subject’s data with the model’s covariance assumptions rather than treating a noisy subject-specific estimate as if it were measured exactly.

The calculation depends on three covariance ingredients: the covariance of the random effects, the design matrix that links those effects to observations, and the residual covariance. The marginal covariance of the observations combines the latter two sources of variation. This tutorial derives the predictor under a Gaussian linear mixed model and clarifies which quantities are assumed known and which are estimated in practice.

Model and Assumptions

Definition: For subject $i$, let $y_i$ be the vector of $n_i$ longitudinal responses. The Gaussian linear mixed model is $$ y_i=X_i\beta+Z_i b_i+\varepsilon_i, $$ where $X_i$ is the fixed-effects design matrix, $\beta$ is the vector of population-level coefficients, $Z_i$ is the random-effects design matrix, $b_i$ is the subject’s random-effects vector, and $\varepsilon_i$ is the residual vector. Assume $b_i\sim N(0,D)$, $\varepsilon_i\sim N(0,R_i)$, and $b_i$ is independent of $\varepsilon_i$. Subjects are independent of one another. The covariance matrices $D$ and $R_i$ are treated as known when deriving the predictor.

The model implies the marginal mean and covariance

E(y_i)=X_i\beta,\qquad V_i=\operatorname{Var}(y_i)=Z_i D Z_i^\mathsf{T}+R_i.

Here $V_i$ is the marginal covariance decomposition: it accounts for both shared variation induced by the random effects and residual variation. In the common special case $R_i=\sigma^2 I$, residual errors have a common variance and are uncorrelated within a subject. More generally, $R_i$ can represent a specified residual covariance structure. The derivation requires $V_i$ to be invertible; the usual formulas also use invertible $D$ and $R_i$, although limiting or generalized-inverse formulations are possible for singular covariance structures.

The Gaussian assumption makes the conditional mean of $b_i$ given the observations a linear function of those observations. The BLUP itself is defined through linear prediction and unbiasedness, and its best-linear property depends on the first two moments; normality is what additionally identifies it as the conditional expectation under the stated model.

Deriving the Predictor

First suppose $\beta$, $D$, and $R_i$ are known. The joint distribution of $b_i$ and $y_i$ is Gaussian. Their covariance is

\operatorname{Cov}(b_i,y_i) =\operatorname{Cov}(b_i,Z_i b_i+\varepsilon_i) =D Z_i^\mathsf{T}.

For jointly Gaussian vectors, the conditional mean of one vector given the other is its mean plus the cross-covariance multiplied by the inverse covariance of the conditioned vector and the centered observed value. Applying that result gives

E(b_i\mid y_i) =D Z_i^\mathsf{T}V_i^{-1}(y_i-X_i\beta).

This is the BLUP when the fixed-effects mean is known. The expression uses the subject’s residual vector relative to the population mean, $y_i-X_i\beta$, and weights it according to the covariances. It is not generally the same as fitting a separate unregularized regression to each subject.

When $\beta$ is unknown, it is estimated by generalized least squares using all subjects. Stack the response vectors and design matrices as $y$, $X$, and $Z$, and let $V=ZDZ^\mathsf{T}+R$ be the marginal covariance of the stacked observations. Then

\hat\beta=(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}y,\qquad \hat b_i=D Z_i^\mathsf{T}V_i^{-1}(y_i-X_i\hat\beta).

The subject blocks of $V$ are independent under the assumptions above, so $V_i$ is the block for subject $i$. The estimate $\hat\beta$ uses data across subjects; replacing $\beta$ with an estimate calculated from only subject $i$ is not generally the BLUP.

A useful justification comes from the mixed-model equations. For the stacked model, the joint estimates solve

\begin{pmatrix} X^\mathsf{T}R^{-1}X & X^\mathsf{T}R^{-1}Z\\ Z^\mathsf{T}R^{-1}X & Z^\mathsf{T}R^{-1}Z+D^{-1} \end{pmatrix} \begin{pmatrix}\hat\beta\\\hat b\end{pmatrix} = \begin{pmatrix}X^\mathsf{T}R^{-1}y\\Z^\mathsf{T}R^{-1}y\end{pmatrix}.

The second block equation rearranges to $\hat b=D Z^\mathsf{T}R^{-1}(y-X\hat\beta-Z\hat b)$. Substitution into the first block equation, or equivalently using the marginal covariance $V=ZDZ^\mathsf{T}+R$, yields the generalized least-squares estimate of $\beta$ and the marginal-covariance form of the predictor above. Under the model, the GLS estimate is unbiased: $E(\hat\beta)=\beta$. Consequently, $E\{y_i-X_i\hat\beta\}=0$, and the predicted random effect has mean zero, matching $E(b_i)=0$. The covariance weighting gives the minimum prediction error variance among linear unbiased predictors; under Gaussian assumptions, the known-$\beta$ version is also the conditional mean.

In applications, $D$ and $R_i$ are usually not known. Their fitted estimates are substituted into the formulas, producing an empirical BLUP. The term “empirical” matters: the displayed best-linear property is exact conditional on the covariance parameters, while estimating those parameters adds uncertainty not represented by simply plugging in their estimates.

Worked Calculations

Worked Example: A Random Intercept

Suppose, hypothetically, that a subject has two measurements from a random-intercept model $y_{ij}=\mu+b_i+\varepsilon_{ij}$. Take the known population mean to be $\mu=10$ response units, random-intercept variance $D=4$ response-units squared, and independent residual variance $R=1$ response-units squared. The observed responses are 13 and 11.

The design vector is $Z_i=(1,1)^\mathsf{T}$, and the residuals from the population mean are $y_i-X_i\beta=(3,1)^\mathsf{T}$. The marginal covariance is

V_i=4 \begin{pmatrix}1&1\\1&1\end{pmatrix} + \begin{pmatrix}1&0\\0&1\end{pmatrix} = \begin{pmatrix}5&4\\4&5\end{pmatrix}.

For two observations with a shared random intercept, the predictor simplifies to

\hat b_i=\frac{D}{R+n_iD}\sum_{j=1}^{n_i}(y_{ij}-\mu) =\frac{4}{1+2(4)}\{(13-10)+(11-10)\} =\frac{16}{9}=1.778\text{ response units}.

The subject’s estimated intercept is therefore $\mu+\hat b_i=11.778$ response units. The estimate is not the raw mean deviation of 2 units: it is reduced because the observed deviations include residual error. With the covariance parameters and mean treated as known, the conditional variance of $b_i$ is $(1/D+n_i/R)^{-1}=4/9$ response-units squared.

Worked Example: Random Intercept and Slope

Consider a hypothetical subject measured at times 0 and 1. The fixed-effects mean is 10 response units at time 0 and 12 at time 1. Let $b_i=(b_{0i},b_{1i})^\mathsf{T}$, with $D=\operatorname{diag}(4,1)$ in the corresponding squared units, and let $R_i=I$, so each residual variance is 1 and the residuals are independent. Suppose the observed responses are 12 and 15. The random-effects design matrix and residual vector are

Z_i=\begin{pmatrix}1&0\\1&1\end{pmatrix},\qquad r_i=y_i-X_i\beta=\begin{pmatrix}2\\3\end{pmatrix}.

Using the equivalent conditional-precision form, the covariance of the random effects conditional on these observations is

C_i=(D^{-1}+Z_i^\mathsf{T}R_i^{-1}Z_i)^{-1} =\left[ \begin{pmatrix}0.25&0\\0&1\end{pmatrix} + \begin{pmatrix}2&1\\1&1\end{pmatrix} \right]^{-1} =\frac{1}{3.5}\begin{pmatrix}2&-1\\-1&2.25\end{pmatrix}.

The conditional mean is $C_i Z_i^\mathsf{T}R_i^{-1}r_i$. Since $Z_i^\mathsf{T}r_i=(5,3)^\mathsf{T}$,

\hat b_i =\frac{1}{3.5} \begin{pmatrix}2&-1\\-1&2.25\end{pmatrix} \begin{pmatrix}5\\3\end{pmatrix} =\begin{pmatrix}2\\0.5\end{pmatrix}.

Thus the predicted random intercept is 2 response units and the predicted random slope is 0.5 response units per time unit. The off-diagonal entries in the conditional covariance calculation arise from the design: both measurements contain the intercept, and the time-1 measurement also contains the slope. The calculation requires the stated fixed mean, $D$, and $R_i$; changing those covariance inputs can change both predictions.

Worked Example: Estimating the Fixed Mean Across Subjects

Suppose, hypothetically, two independent subjects each contribute one measurement under a random-intercept model with common population mean $\mu$, random-intercept variance 4, and residual variance 1. The measurements are $y_1=12$ and $y_2=8$. Because each marginal variance is $4+1=5$, generalized least squares gives the equally weighted mean

\hat\mu= \frac{12/5+8/5}{1/5+1/5} =10.

The residuals from this estimated population mean are 2 and $-2$. For one observation per subject, the shrinkage multiplier is $D/(D+R)=4/5$, so

\hat b_1=\frac{4}{5}(12-10)=1.6,\qquad \hat b_2=\frac{4}{5}(8-10)=-1.6.

The calculations show why the fixed-effects estimate is part of the subject-effect prediction: the deviations are defined relative to the fitted population mean. Here the covariance parameters are treated as known, and the same data are used to estimate the common mean and predict both subject effects.

What Covariances Are Needed?

To calculate a BLUP, identify the following quantities and their roles:

$D$: the covariance matrix of the subject’s random effects. It describes their population variation and, when non-diagonal, their covariance.
$Z_i$: the subject’s random-effects design matrix. It maps random effects onto the observed measurement times and covariate values.
$R_i$: the within-subject residual covariance matrix. It specifies residual variances and any modeled residual correlation.
$V_i=Z_iDZ_i^\mathsf{T}+R_i$: the marginal covariance of the observations. It weights the observed deviations in the BLUP formula.
$X_i\beta$: the fixed-effects mean, or $X_i\hat\beta$ when the fixed effects are estimated. The predictor uses deviations from this mean.

A random-intercept-and-slope model, for example, requires a $2\times2$ random-effects covariance matrix, not just separate intercept and slope variances, unless the model explicitly sets their covariance to zero. Likewise, substituting an independent-residual matrix for a model with correlated residuals changes $V_i$ and therefore changes the prediction.

Diagnostics, Alternatives, and Limitations

A BLUP is a model-based predictor, so it can be misleading when the model’s covariance structure is poorly supported. Check that the design matrices encode the intended time scale and random effects, and that the fitted covariance estimates are admissible. A covariance estimate at or near a boundary, a singular random-effects covariance, or a nearly non-invertible marginal covariance can make predictions unstable or reveal that the specified random-effects structure is not supported by the data. Sparse observations for a subject also mean that the prediction relies more heavily on the assumed population covariance than on that subject’s data.

The main alternative is a subject-specific estimate obtained by fitting a separate regression to each subject. Such estimates do not generally account for the population distribution of effects, and may be very noisy when each subject has few observations. The BLUP instead combines subject-level information with the assumed distribution and residual variation. This is useful for prediction, but does not make the predicted effect an observed fact or establish a causal subject-specific effect.

A common mistake is to calculate a separate subject mean or regression coefficient and call it the BLUP. That ignores the covariance weighting. Another is to use $D$ but omit $R_i$, or to treat $V_i$ as the random-effects covariance. The correction is to form $V_i=Z_iDZ_i^\mathsf{T}+R_i$, center the observations at the fixed-effects mean, and apply $D Z_i^\mathsf{T}V_i^{-1}$. With estimated covariance parameters, this is an empirical BLUP; its plug-in formula does not by itself capture all uncertainty from estimating those parameters.

Key takeaway: The subject-effect BLUP is a covariance-weighted prediction from the subject’s deviations from the population mean. Its calculation requires the random-effects covariance $D$, random-effects design $Z_i$, residual covariance $R_i$, marginal covariance $V_i$, and the fixed-effects mean.

Check Your Understanding

Answer each question, then compare with the solutions below.

In the Gaussian mixed model, what is $V_i$, and which two sources of variation does it combine?
For a random-intercept model with known $D=3$, $R=2$, two observations, and a sum of residuals from the fixed mean equal to 5, calculate the BLUP.
Why does the BLUP use $X_i\hat\beta$ rather than simply using the observed responses as deviations?
What covariance information is missing if a random-intercept-and-slope model specifies both variances but not their covariance?

Solutions. (1) $V_i=Z_iDZ_i^\mathsf{T}+R_i$. It combines variation induced by the random effects and residual variation. (2) The random-intercept formula gives $3/(2+2(3))\times5=15/8=1.875$ response units. (3) Random effects are deviations from the population-level mean, so the response must be centered at its fixed-effects mean; otherwise the population mean is incorrectly attributed to the subject effect. (4) The off-diagonal entry of $D$, the intercept-slope covariance, is unspecified. It must be estimated or explicitly constrained, for example to zero, to complete the covariance model.