Approximate Population Means and Covariances

Why Approximate Population Moments Matter

A population design concerns measurements across subjects, not just the predictions for one subject with typical parameters. Subjects vary in their random effects, and each observation also contains residual error. To describe the resulting population-level variation in a planned observation vector, we need to combine both sources of variability while accounting for how the model responds to random effects.

The earlier material on first-order Taylor expansion around zero random effects supplies the local linear approximation used here. This tutorial applies that approximation to obtain the marginal mean and covariance of a subject’s repeated observations. These moments describe the approximate population distribution without requiring the nonlinear model to be integrated exactly over the random-effects distribution.

Definition: For a subject measured at a fixed design, let \(\mathbf{Y}_i\) be the \(q\)-dimensional observation vector, \(\boldsymbol{\eta}_i\) the \(r\)-dimensional random-effects vector, and \(\boldsymbol{\theta}\) the fixed-effect parameters. Write the conditional mean as \(\mathbf{f}(\boldsymbol{\theta},\boldsymbol{\eta}_i)\). Assume \(E(\boldsymbol{\eta}_i)=\mathbf{0}\), \(\operatorname{Cov}(\boldsymbol{\eta}_i)=\boldsymbol{\Omega}\), and that the residual vector has conditional mean zero and conditional covariance \(\mathbf{R}(\boldsymbol{\eta}_i)\). The approximate population mean and covariance are \(\boldsymbol{\mu}_{\mathrm{pop}}\approx\mathbf{f}_0\) and \(\boldsymbol{\Sigma}_{\mathrm{pop}}\approx\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}+\mathbf{R}_0\), where \(\mathbf{f}_0=\mathbf{f}(\boldsymbol{\theta},\mathbf{0})\), \(\mathbf{G}\) is the derivative of \(\mathbf{f}\) with respect to \(\boldsymbol{\eta}\) at zero, and \(\mathbf{R}_0=\mathbf{R}(\mathbf{0})\).

From Linearized Predictions to Marginal Moments

For the chosen design, the \(j\)th component of \(\mathbf{f}\) is the predicted mean at observation \(j\). The \(q\)-by-\(r\) sensitivity matrix is \(\mathbf{G}=\left.\partial\mathbf{f}/\partial\boldsymbol{\eta}^{T}\right|_{\boldsymbol{\eta}=\mathbf{0}}\). Its entry \(G_{jk}\) describes the local change in the mean of observation \(j\) for a unit change in random effect \(k\), with other random effects fixed. The first-order approximation replaces the nonlinear conditional mean by \(\mathbf{f}_0+\mathbf{G}\boldsymbol{\eta}_i\).

\mathbf{Y}_i \approx \mathbf{f}_0+\mathbf{G}\boldsymbol{\eta}_i+\boldsymbol{\varepsilon}_i

Here \(\boldsymbol{\varepsilon}_i\) has mean zero and, at the linearization point, covariance \(\mathbf{R}_0\). For the covariance formula, assume the residual error is uncorrelated with the random effects; conditional mean zero for the residual implies this condition when the conditional model is correctly specified. The residual components need not be independent of one another: their dependence is represented by the off-diagonal entries of \(\mathbf{R}_0\).

Take the expectation.
Since the random effects and residuals have mean zero, \(E(\mathbf{Y}_i)\approx\mathbf{f}_0\).

Center the linearized model.
Subtract its mean to obtain \(\mathbf{Y}_i-E(\mathbf{Y}_i)\approx\mathbf{G}\boldsymbol{\eta}_i+\boldsymbol{\varepsilon}_i\).

Expand the covariance.
The random-effect term contributes \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}\); the residual term contributes \(\mathbf{R}_0\); the cross-covariance terms are zero under the stated assumption.

In detail, the variance of a linear transformation is \(\operatorname{Cov}(\mathbf{G}\boldsymbol{\eta}_i)=\mathbf{G}\operatorname{Cov}(\boldsymbol{\eta}_i)\mathbf{G}^{T}\). Expanding the covariance of the sum gives

\begin{aligned} \operatorname{Cov}(\mathbf{Y}_i) &\approx \mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T} +\mathbf{R}_0 +\operatorname{Cov}(\mathbf{G}\boldsymbol{\eta}_i,\boldsymbol{\varepsilon}_i) +\operatorname{Cov}(\boldsymbol{\varepsilon}_i,\mathbf{G}\boldsymbol{\eta}_i)\\ &=\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}+\mathbf{R}_0. \end{aligned}

The \(j,k\) entry of the resulting covariance matrix describes the approximate covariance between observations \(j\) and \(k\) from the same subject. In particular, \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}\) can produce covariance between different observation times even when residual errors at those times are independent. The diagonal entries are marginal variances, including both random-effect and residual contributions.

Key distinction: \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}\) is the covariance induced by between-subject random effects under the linearization. It is not, by itself, the total observation covariance. The latter also includes \(\mathbf{R}_0\).

Worked Examples

Worked Example: One Random Effect and Two Observations

Suppose a hypothetical subject has two planned concentration measurements. A random effect \(\eta_i\) acts as a multiplicative shift on the typical concentration, so the conditional mean at time \(j\) is \(f_j(\eta_i)=m_j\exp(\eta_i)\). Let the typical predictions be \(m_1=10\) mg/L and \(m_2=5\) mg/L, and let \(\eta_i\) have mean zero and variance \(0.04\). Assume independent additive residual errors, each with variance \(0.25\) (mg/L)\(^2\).

At \(\eta_i=0\), the mean vector is \(\mathbf{f}_0=(10,5)^T\) mg/L. The derivative of \(m_j\exp(\eta_i)\) with respect to \(\eta_i\), evaluated at zero, is \(m_j\). Thus \(\mathbf{G}=(10,5)^T\) mg/L and \(\boldsymbol{\Omega}=(0.04)\). The residual covariance is \(\mathbf{R}_0=0.25\mathbf{I}\) (mg/L)\(^2\).

The random-effect contribution is \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}=0.04(10,5)^T(10,5)\). Applying the formula gives

\boldsymbol{\Sigma}_{\mathrm{pop}} \approx \begin{pmatrix}4&2\\2&1\end{pmatrix} + \begin{pmatrix}0.25&0\\0&0.25\end{pmatrix} = \begin{pmatrix}4.25&2\\2&1.25\end{pmatrix} \;(\mathrm{mg/L})^2.

The approximate marginal means are 10 and 5 mg/L. The marginal standard deviations are \(\sqrt{4.25}=2.062\) mg/L and \(\sqrt{1.25}=1.118\) mg/L. The covariance between the two measurements is \(2\;(\mathrm{mg/L})^2\), arising entirely from the shared random effect because the residual errors are independent. Their approximate correlation is \(2/\sqrt{4.25(1.25)}=0.8677\). This is a moment calculation: it does not establish that the exact nonlinear observation distribution is normal.

Worked Example: Two Correlated Random Effects

Consider a hypothetical two-observation design with zero-random-effect mean vector \(\mathbf{f}_0=(8,4)^T\) response units. Suppose two random effects have covariance \(\boldsymbol{\Omega}=\begin{pmatrix}0.09&0.03\\0.03&0.04\end{pmatrix}\), and the sensitivity matrix is \(\mathbf{G}=\begin{pmatrix}2&1\\1&3\end{pmatrix}\) response units per random-effect unit. Let the residual covariance be \(\mathbf{R}_0=\begin{pmatrix}0.16&0.04\\0.04&0.25\end{pmatrix}\) response-units squared.

First multiply the sensitivity matrix by the random-effect covariance: \(\mathbf{G}\boldsymbol{\Omega}=\begin{pmatrix}0.21&0.10\\0.18&0.15\end{pmatrix}\). Therefore,

\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T} = \begin{pmatrix}0.21&0.10\\0.18&0.15\end{pmatrix} \begin{pmatrix}2&1\\1&3\end{pmatrix} = \begin{pmatrix}0.52&0.51\\0.51&0.63\end{pmatrix}.

Adding the residual covariance gives \(\boldsymbol{\Sigma}_{\mathrm{pop}}\approx\begin{pmatrix}0.68&0.55\\0.55&0.88\end{pmatrix}\) response-units squared. The total off-diagonal covariance is \(0.51+0.04=0.55\), combining covariance induced by the random effects with residual covariance. The approximate marginal means remain \((8,4)^T\). The positive off-diagonal entries in \(\boldsymbol{\Omega}\) matter because the random effects co-vary, but their influence on observation covariance also depends on the sensitivities in \(\mathbf{G}\).

Worked Example: Proportional Residual Error

Suppose two hypothetical typical concentrations are \(m_1=12\) mg/L and \(m_2=6\) mg/L, with a shared multiplicative random effect having variance \(0.01\). Use the same local sensitivity form as in the first example, so \(\mathbf{G}=(12,6)^T\) mg/L. Assume independent residual errors with variance \(\sigma_a^2+\sigma_p^2m_j^2\) at observation \(j\), where the additive residual standard deviation is \(\sigma_a=0.5\) mg/L and the proportional residual standard deviation is \(\sigma_p=0.1\).

Evaluate the residual variance at each typical prediction, as required by \(\mathbf{R}_0\). At the first observation it is \(0.5^2+(0.1)^2(12)^2=0.25+1.44=1.69\) (mg/L)\(^2\). At the second it is \(0.5^2+(0.1)^2(6)^2=0.25+0.36=0.61\) (mg/L)\(^2\). Hence \(\mathbf{R}_0=\operatorname{diag}(1.69,0.61)\) (mg/L)\(^2\).

The random-effect contribution is \(0.01(12,6)^T(12,6)=\begin{pmatrix}1.44&0.72\\0.72&0.36\end{pmatrix}\) (mg/L)\(^2\). Adding the residual contribution gives \(\boldsymbol{\Sigma}_{\mathrm{pop}}\approx\begin{pmatrix}3.13&0.72\\0.72&0.97\end{pmatrix}\) (mg/L)\(^2\). The proportional component makes residual variance larger at the higher typical concentration. It contributes no off-diagonal covariance here because residual errors were assumed independent.

Assumptions, Checks, and Limits

The formula is useful only when its assumptions and local approximation are appropriate. The random effects must be centered at zero with covariance \(\boldsymbol{\Omega}\), and the conditional observation mean must be differentiable with respect to them near zero. The residual model must have conditional mean zero. The formula as written evaluates residual covariance at zero random effects and assumes no residual–random-effect cross-covariance. If residual variance changes materially across the random-effects distribution, using only \(\mathbf{R}_0\) may be a poor approximation to the marginal residual contribution.

The first-order approximation also omits curvature in the conditional mean. For a sufficiently smooth \(j\)th prediction, a second-order Taylor expansion has a term \(\tfrac12\boldsymbol{\eta}^{T}\mathbf{H}_j\boldsymbol{\eta}\), where \(\mathbf{H}_j\) is its Hessian with respect to the random effects at zero. Its expectation is \(\tfrac12\operatorname{tr}(\mathbf{H}_j\boldsymbol{\Omega})\). Thus a nonlinear model can have a marginal mean different from \(\mathbf{f}_0\), even though \(E(\boldsymbol{\eta})=\mathbf{0}\). The first-order result omits that correction; it is most credible when the prediction is nearly linear over the random-effect values with appreciable probability.

Several numerical checks help identify mistakes. \(\boldsymbol{\Omega}\) and \(\mathbf{R}_0\) should be symmetric positive semidefinite covariance matrices, and the computed \(\boldsymbol{\Sigma}_{\mathrm{pop}}\) should also be symmetric positive semidefinite. Confirm dimensions: if there are \(q\) observations and \(r\) random effects, then \(\mathbf{G}\) is \(q\)-by-\(r\), so \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}\) is \(q\)-by-\(q\). Check that every covariance has squared response units, while entries of \(\mathbf{G}\) have response units per random-effect unit.

Inspect the implied marginal variances and correlations as well. A negative variance, a correlation outside \([-1,1]\), or an unexpected asymmetry usually signals a coding, indexing, or matrix-multiplication error. Conversely, a valid covariance matrix does not prove that the approximation is accurate: compare predictions across plausible random-effect values or compare with numerical integration or simulation when curvature or heteroscedasticity appears substantial.

This moment calculation differs from two common alternatives. A conditional prediction at \(\boldsymbol{\eta}=\mathbf{0}\) gives only a typical-subject mean and says nothing by itself about population variation. A model that ignores random effects but retains residual error uses \(\mathbf{R}_0\) alone and understates marginal variability when \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}\) is appreciable. More elaborate integration or simulation can represent nonlinear effects more faithfully, but costs more computation. The first-order moments are a transparent local approximation, not a guarantee of exact population moments or a complete description of the distribution.

Common Pitfalls

Leaving out residual error. The random-effect term describes between-subject variation in the linearized prediction; add \(\mathbf{R}_0\) to obtain total observation covariance.
Adding variances without sensitivities. Random-effect variances are on the random-effect scale. They must be mapped to observation scale through \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}\).
Assuming zero-mean random effects force the exact marginal mean to equal \(\mathbf{f}_0\). That equality holds for the linearized model. Nonlinear curvature can shift the exact mean.
Confusing correlation sources. Off-diagonal covariance can arise from shared random effects, correlated residuals, or both. State which source produces each term.
Treating the approximation as an exact distribution. Approximate means and covariances do not, by themselves, establish that the original nonlinear observations follow a multivariate normal distribution.

Key takeaway: Linearize the conditional mean at zero random effects, propagate random-effect covariance through the sensitivity matrix, and add residual covariance evaluated at that point. The result is an approximate marginal mean and covariance for the subject’s planned observation vector.

Check Your Understanding

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

If \(\mathbf{G}\) is \(4\)-by-\(2\), what are the dimensions of \(\boldsymbol{\Omega}\), \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}\), and \(\mathbf{R}_0\)?
A one-observation model has sensitivity \(G=3\) response units per random-effect unit, random-effect variance \(0.04\), and residual variance \(0.25\) response-units squared. Find the approximate marginal variance.
For two observations, the random-effect contribution to covariance is \(0.8\;(\mathrm{mg/L})^2\) and the residual covariance is \(-0.2\;(\mathrm{mg/L})^2\). What is the total off-diagonal covariance, and what do the two terms represent?
Why can the exact marginal mean differ from the prediction at zero random effects even when the random effects have mean zero?

Solutions. (1) \(\boldsymbol{\Omega}\) is \(2\)-by-\(2\); \(\mathbf{G}\boldsymbol{\Omega}\mathbf{G}^{T}\) and \(\mathbf{R}_0\) are each \(4\)-by-\(4\). (2) The random-effect contribution is \(3^2(0.04)=0.36\) response-units squared, so the total is \(0.36+0.25=0.61\) response-units squared. (3) The total is \(0.8+(-0.2)=0.6\;(\mathrm{mg/L})^2\). The first term is induced by random effects; the second is residual covariance. (4) A nonlinear prediction can have curvature: its second-order terms can have nonzero expectation under a zero-mean random-effects distribution.