Marginal Covariance From Linearized Individual Predictions

Why Individual Predictions Create Covariance

Repeated observations from one individual are linked even when their residual errors are independent. The individual’s random effects enter multiple predictions, so a subject whose random effects shift one prediction will often have other predictions shifted as well. The resulting covariance depends on how each prediction responds to each random effect, not just on the random-effect variances themselves.

This tutorial derives a first-order approximation to that covariance. It is a local approximation: predictions are linearized around zero random effects, and the covariance of that linear approximation is combined with residual covariance. The result is useful for understanding the structure of population observations and for checking what features of a model contribute to their variability.

Definition: For individual \(i\), let \(\boldsymbol{\eta}_i\) be a vector of mean-zero random effects with covariance matrix \(\Omega\). Let \(\mathbf{f}_i(\boldsymbol{\eta}_i)\) be the vector of model predictions at that individual’s observation times and covariate values, and let \(\boldsymbol{\epsilon}_i\) be the residual-error vector. The observation vector is \(\mathbf{y}_i=\mathbf{f}_i(\boldsymbol{\eta}_i)+\boldsymbol{\epsilon}_i\). The Jacobian \(J_i\) is the matrix of first derivatives of the prediction vector with respect to the random effects, evaluated at \(\boldsymbol{\eta}_i=\mathbf{0}\).

The subscript \(i\) holds the individual’s design fixed: the observation times, covariates, and fixed-effect parameters are treated as given. Thus, the covariance derived here describes variation across individuals with the same design, including covariance among repeated observations on one individual. It is not a covariance between observations from two independent individuals.

Linearization and the Covariance Formula

Suppose \(\mathbf{f}_i\) is differentiable near \(\boldsymbol{\eta}_i=\mathbf{0}\). A first-order Taylor expansion gives

\mathbf{f}_i(\boldsymbol{\eta}_i) \approx \mathbf{f}_i(\mathbf{0})+J_i\boldsymbol{\eta}_i.

Here, \(J_i\) has one row for each observation and one column for each random effect. Its entry in row \(j\), column \(k\) is the derivative of prediction \(j\) with respect to random effect \(k\), evaluated at zero random effects. The random-effect covariance \(\Omega\) is therefore transformed into observation space by multiplication on both sides by the Jacobian.

Result: If \(E(\boldsymbol{\eta}_i)=\mathbf{0}\), \(\operatorname{Cov}(\boldsymbol{\eta}_i)=\Omega\), \(E(\boldsymbol{\epsilon}_i\mid\boldsymbol{\eta}_i)=\mathbf{0}\), and the residual covariance is the fixed matrix \(R_i\), then the first-order approximation is

\operatorname{Cov}(\mathbf{y}_i) \approx J_i\Omega J_i^{T}+R_i.

The assumptions in this result are substantive. The prediction function must be adequately represented by its first-order expansion over the random-effect values that matter. The residual vector has conditional mean zero, its covariance is \(R_i\), and it is uncorrelated with the random effects under the stated conditional-mean assumption. If residual covariance varies with random effects, it cannot in general be replaced by a fixed \(R_i\) without further approximation.

Why the Formula Holds

Substitute the linearized prediction into the observation model. The approximation becomes \(\mathbf{y}_i\approx \mathbf{f}_i(\mathbf{0})+J_i\boldsymbol{\eta}_i+\boldsymbol{\epsilon}_i\). A fixed vector such as \(\mathbf{f}_i(\mathbf{0})\) does not contribute to covariance. Expanding the covariance of the remaining terms gives

\begin{aligned} \operatorname{Cov}(\mathbf{y}_i) &\approx \operatorname{Cov}(J_i\boldsymbol{\eta}_i) +\operatorname{Cov}(\boldsymbol{\epsilon}_i)\\ &\quad+\operatorname{Cov}(J_i\boldsymbol{\eta}_i,\boldsymbol{\epsilon}_i) +\operatorname{Cov}(\boldsymbol{\epsilon}_i,J_i\boldsymbol{\eta}_i). \end{aligned}

For a fixed matrix \(J_i\), \(\operatorname{Cov}(J_i\boldsymbol{\eta}_i)=J_i\Omega J_i^T\). The two cross-covariance terms are zero because \(E(\boldsymbol{\epsilon}_i\mid\boldsymbol{\eta}_i)=\mathbf{0}\). The residual covariance is \(R_i\), giving the stated formula. This argument also explains why the whole matrix \(\Omega\) matters: off-diagonal elements can contribute to observation variances and covariances.

The associated first-order mean is \(\mathbf{f}_i(\mathbf{0})\). This is not generally the exact marginal mean of a nonlinear prediction. Curvature can shift the mean as well as change the covariance; both effects are omitted by a first-order expansion.

Worked Example: Two Observations Sharing One Random Effect

Consider a hypothetical individual with two concentration predictions, in mg/L, given by \(f_1(\eta)=10\exp(\eta)\) and \(f_2(\eta)=6\exp(\eta)\). The random effect has mean zero and variance \(\omega=0.09\). Assume independent, mean-zero residual errors with variances \(1\) and \(4\) \((\mathrm{mg/L})^2\). The exponential functions are differentiable, and for this first-order calculation we evaluate their derivatives at \(\eta=0\).

The predictions at zero are \(10\) and \(6\) mg/L. Their derivatives are \(f'_1(0)=10\) and \(f'_2(0)=6\), so

J= \begin{bmatrix} 10\\ 6 \end{bmatrix}, \qquad \Omega=[0.09], \qquad R= \begin{bmatrix} 1&0\\ 0&4 \end{bmatrix}.

The random-effect contribution is

J\Omega J^T =0.09 \begin{bmatrix} 100&60\\ 60&36 \end{bmatrix} = \begin{bmatrix} 9&5.4\\ 5.4&3.24 \end{bmatrix}.

Adding the residual covariance gives

\operatorname{Cov} \begin{pmatrix}Y_1\\Y_2\end{pmatrix} \approx \begin{bmatrix} 10&5.4\\ 5.4&7.24 \end{bmatrix} \;(\mathrm{mg/L})^2.

The approximate covariance is \(5.4\;(\mathrm{mg/L})^2\), even though the residual errors are independent. It arises because the same \(\eta\) affects both predictions. The corresponding approximate correlation is \(5.4/\sqrt{10(7.24)}=0.635\). Thus, the predictions tend to move together through their shared random effect.

Correlated Effects and Prediction Sensitivity

When there are multiple random effects, each prediction may respond to several of them. The derivatives express that sensitivity, while \(\Omega\) describes their variances and covariances. A negative derivative for one effect does not by itself guarantee a negative contribution to the prediction variance: the complete quadratic form must be evaluated.

Worked Example: Correlated Effects Can Reduce Variance

Suppose a hypothetical concentration prediction, in mg/L, is \(f(\eta_1,\eta_2)=20\exp(\eta_1-0.5\eta_2)\). The random effects have covariance matrix

\Omega= \begin{bmatrix} 0.04&0.03\\ 0.03&0.09 \end{bmatrix}.

This is a valid covariance matrix: its diagonal entries are positive and its determinant is \(0.04(0.09)-0.03^2=0.0027>0\). The variances are \(0.04\) and \(0.09\), and the covariance is \(0.03\). At \((\eta_1,\eta_2)=(0,0)\), the gradient is \(\mathbf{g}=[20,-10]\), in mg/L per unit of the dimensionless random effects.

For one prediction, the random-effect variance is the scalar \(\mathbf{g}\Omega\mathbf{g}^T\). Substitution gives

\begin{aligned} \mathbf{g}\Omega\mathbf{g}^T &=20^2(0.04)+(-10)^2(0.09) +2(20)(-10)(0.03)\\ &=16+9-12=13\;(\mathrm{mg/L})^2. \end{aligned}

The cross term is negative because the prediction responds in opposite directions to the two effects while the effects have positive covariance. If the residual variance at this observation is \(2.25\;(\mathrm{mg/L})^2\), the approximate total variance is \(13+2.25=15.25\;(\mathrm{mg/L})^2\). Dropping the covariance in \(\Omega\) would instead give \(16+9+2.25=27.25\;(\mathrm{mg/L})^2\), a substantially different result.

When the Local Approximation Is Weak

A first-order approximation is most credible when predictions are close to linear over the random-effect values that occur in the population. Large random-effect variability, strong curvature, or predictions that change sharply over a small range can make the derivative at zero unrepresentative. In those settings, \(J_i\Omega J_i^T+R_i\) remains the covariance of the linearized model, but may be a poor approximation to the covariance of the nonlinear model.

Worked Example: Checking Curvature Against an Exact Calculation

Let a hypothetical scalar prediction be \(Y=\exp(\eta)+\epsilon\), where \(\eta\sim N(0,0.25)\) and the independent residual has mean zero and variance \(0.04\). The derivative of \(\exp(\eta)\) at zero is \(1\), so the first-order approximation gives prediction variance \(1^2(0.25)=0.25\) and total variance \(0.25+0.04=0.29\), in squared response units.

For comparison, the exact lognormal moments are \(E[\exp(\eta)]=\exp(0.25/2)=1.133\) and \(\operatorname{Var}[\exp(\eta)]=\exp(0.25)\{\exp(0.25)-1\}=0.3647\). Therefore the exact total variance is \(0.3647+0.04=0.4047\). The first-order result is lower by \(0.1147\) squared response units. This difference illustrates that a valid Taylor expansion at zero need not be accurate across the full distribution of \(\eta\).

A practical diagnostic is to compare the linearized covariance with predictions evaluated at a representative set of random-effect draws from the specified distribution. Large discrepancies, marked asymmetry in prediction changes, or strong variation in local derivatives suggest that the first-order approximation needs caution. This comparison is a diagnostic, not a proof that a particular approximation is accurate.

The main alternatives are to integrate the nonlinear model over the random-effect distribution or to use simulation to approximate its marginal moments. These approaches can capture curvature that the first-order formula omits, but require more computation and careful numerical evaluation. Linearization is transparent and often useful for interpretation; it should not be treated as exact merely because its matrix calculation is straightforward.

Interpretation and Common Pitfalls

The matrix \(J_i\Omega J_i^T\) describes approximate covariance induced by shared random effects, and \(R_i\) describes residual covariance. If residual errors are independent across observations, \(R_i\) is diagonal, so off-diagonal observation covariance in this approximation comes from shared random effects. If residual covariance is not diagonal, its off-diagonal terms must also be included. For observations from different individuals with independent random effects and residuals, the cross-individual covariance is zero under this model.

Using random-effect variances alone. This discards covariances between effects. Use the full \(\Omega\), including its off-diagonal entries.
Adding residual variance to every matrix entry. Add the residual covariance matrix \(R_i\); under independent residual errors, only its diagonal entries are nonzero.
Calling the result exact for a nonlinear model. The result is the covariance of a first-order approximation. Check whether curvature is important.
Confusing a shared-effect covariance with residual correlation. Two observations can covary because they share random effects even when their residual errors are independent.

Key takeaway: Linearizing individual predictions gives the approximate marginal covariance \(J_i\Omega J_i^T+R_i\). The Jacobian translates random-effect variation into observation variation; the residual covariance is then added under the stated assumptions.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A prediction has gradient \([3,2]\), random-effect covariance matrix \(\begin{bmatrix}0.04&0.01\\0.01&0.09\end{bmatrix}\), and residual variance \(0.5\). What is its approximate total variance?
Two observations share one random effect with variance \(0.16\). Their derivatives with respect to that effect are \(5\) and \(-2\), and their residual errors are independent. What is their approximate covariance?
State two assumptions needed to add a fixed residual covariance matrix to the linearized random-effect covariance.
Why can an accurate derivative at zero still produce a poor approximation to marginal covariance?

Solutions. (1) The random-effect contribution is \(3^2(0.04)+2^2(0.09)+2(3)(2)(0.01)=0.36+0.36+0.12=0.84\). Adding residual variance gives \(0.84+0.5=1.34\) squared response units. (2) The covariance is \(5(-2)(0.16)=-1.6\) squared response units; independent residual errors add no off-diagonal covariance. (3) For example, the residual vector must have conditional mean zero given the random effects, and its covariance must be the specified fixed matrix \(R_i\); the first-order prediction expansion must also be adequate for using this approximation. (4) The prediction may have substantial curvature across the random-effect distribution, so its local slope at zero does not represent its behavior over the values that occur.