Tutorials › Nonlinear Mixed-Effects Models (NLME) › FO With Multiple Observations Per Individual

The Marginal Likelihood and Its Approximations · Tutorial 150 of 797

FO With Multiple Observations Per Individual

Learn how the first-order linearization turns shared random effects into within-individual covariance and how to use that covariance in an FO likelihood.

Advanced 12 min read

What You'll Learn

  • Write the FO linearization for a vector of observations from one individual
  • Derive the marginal covariance, including the condition on cross-covariance
  • Calculate induced correlations when residual errors are uncorrelated
  • Assemble and evaluate a multivariate-normal FO subject likelihood
  • Recognize when linearization or covariance assumptions may be inadequate
  • Distinguish FO covariance from an exact nonlinear marginal distribution

Why Multiple Observations Change the FO Calculation

With one observation per individual, a random effect contributes variability to that observation. With several observations from the same individual, the same random effect enters each prediction. The resulting deviations across observations are therefore generally dependent, even if the observation-level residual errors are uncorrelated. This dependence is central to constructing the first-order (FO) subject likelihood.

The key calculation is vector-valued: linearize all of an individual’s predictions with respect to that individual’s random effects, then propagate the random-effect covariance through the resulting Jacobian. The covariance is an approximation to the nonlinear model’s marginal covariance; it is not, by itself, a guarantee that the marginal distribution is exactly normal.

FO Linearization and the Subject Covariance

Definition: For individual \(i\), let \(y_i\) be the \(n_i\)-vector of observations, \(\eta_i\) the \(q\)-vector of random effects, and \(f_i(\theta,\eta_i)\) the \(n_i\)-vector of model predictions. A first-order Taylor expansion about \(\eta_i=0\) is $$ f_i(\theta,\eta_i)\approx f_i(\theta,0)+J_i\eta_i, $$ where \(J_i\) is the \(n_i\times q\) Jacobian matrix of derivatives of the prediction vector with respect to the random effects, evaluated at \(\eta_i=0\). If \(E(\eta_i)=0\), the residual vector \(\epsilon_i\) has mean zero and covariance \(R_i\), and the cross-covariance between \(J_i\eta_i\) and \(\epsilon_i\) is zero, then the linearized observation vector has covariance $$ V_i=J_i\Omega J_i^{\mathsf T}+R_i, $$ where \(\Omega=\operatorname{Var}(\eta_i)\). Independence of \(\eta_i\) and \(\epsilon_i\) is sufficient for the zero cross-covariance condition, but is not necessary.

To see the source of this result, write the linearized deviation from the prediction at zero random effects as \(u_i=J_i\eta_i+\epsilon_i\). Its covariance expands as $$ \operatorname{Var}(u_i)=J_i\Omega J_i^{\mathsf T}+R_i+ \operatorname{Cov}(J_i\eta_i,\epsilon_i)+ \operatorname{Cov}(\epsilon_i,J_i\eta_i). $$ Thus the stated covariance formula follows when the two cross-covariance terms are zero. This condition is weaker than independence. If those terms are not zero, they must be included; simply adding the two marginal covariance matrices is then incorrect.

The \((j,k)\) entry of the random-effect contribution is $$ \left(J_i\Omega J_i^{\mathsf T}\right)_{jk} =J_{ij}\Omega J_{ik}^{\mathsf T}, $$ where \(J_{ij}\) is row \(j\) of \(J_i\). This expression links observations \(j\) and \(k\) through their sensitivities to the same random effects. If \(R_i\) is diagonal, its off-diagonal entries are zero, but \(V_i\) need not be diagonal. In other words, the marginal deviations \(y_i-f_i(\theta,0)\) can be correlated even when the residual errors themselves are uncorrelated.

The covariance calculation only requires the stated means, covariances, linearization, and zero cross-covariance. To use a multivariate-normal FO subject likelihood, an additional distributional condition is needed: it is sufficient for \(\eta_i\) and \(\epsilon_i\) to be Gaussian and independent, but that is not necessary. More generally, the likelihood is justified as a multivariate normal when the linearized sum \(J_i\eta_i+\epsilon_i\) is Gaussian. Without such a condition, \(V_i\) still describes the covariance of the linearized sum, but does not by itself establish its distribution.

Worked Example: A Shared Random Effect Across Three Observations

Suppose three hypothetical observations from one individual have predictions at \(\eta_i=0\) of \(10\), \(8\), and \(6\) concentration units. There is one random effect with variance \(\Omega=0.25\), and the derivative of each prediction with respect to that effect is \(4\), \(2\), and \(1\) concentration units per random-effect unit. Assume additive residual errors are independent, have variance \(1\) concentration-unit squared, and are independent of the random effect. These assumptions make the cross-covariance zero and give \(R_i=I_3\).

The Jacobian is \(J_i=(4,2,1)^{\mathsf T}\). The random-effect covariance contribution is $$ J_i\Omega J_i^{\mathsf T} =0.25 \begin{pmatrix}4\\2\\1\end{pmatrix} \begin{pmatrix}4&2&1\end{pmatrix} = \begin{pmatrix} 4&2&1\\ 2&1&0.5\\ 1&0.5&0.25 \end{pmatrix}. $$ Adding the residual covariance gives $$ V_i= \begin{pmatrix} 5&2&1\\ 2&2&0.5\\ 1&0.5&1.25 \end{pmatrix}. $$

For example, the covariance between observations 1 and 2 is \(4(0.25)(2)=2\) concentration-unit squared. Their approximate marginal correlation is $$ \frac{2}{\sqrt{5(2)}}=\frac{2}{\sqrt{10}}\approx 0.6325. $$ The residual errors are uncorrelated by assumption, but the two observation deviations are positively correlated because an upward shift in the shared random effect moves both predictions upward. This conclusion concerns the linearized, marginal deviations, not the residual errors alone.

How Random-Effect Structure Shapes the Covariance

A random intercept often gives a particularly transparent pattern: if every observation has the same derivative with respect to that intercept, then the random-effect contribution adds the same covariance to every pair of observations. A random slope, by contrast, contributes differently according to the observation times or covariates. With several random effects, their covariance—including any correlation among them—determines how these contributions combine.

Worked Example: Random Intercept and Random Slope

Consider hypothetical observations at times \(t=0,1,2\). Suppose the linearized prediction has sensitivities \(1\) to a random intercept and \(t\) to a random slope, so the Jacobian rows are \((1,t)\). Let $$ \Omega= \begin{pmatrix} 0.16&0.03\\ 0.03&0.09 \end{pmatrix}, \qquad R_i=0.04I_3. $$ The random effects have positive variances and \(\det(\Omega)=0.16(0.09)-0.03^2=0.0135>0\), so \(\Omega\) is positive definite. Assume zero-mean errors and zero cross-covariance between the random-effect contribution and residual vector; independent random effects and residuals would be one sufficient way to meet the latter condition.

The Jacobian is $$ J_i= \begin{pmatrix} 1&0\\ 1&1\\ 1&2 \end{pmatrix}. $$ The random-effect variances at the three times are \(0.16\), \(0.16+2(0.03)+0.09=0.31\), and \(0.16+4(0.03)+4(0.09)=0.64\). The corresponding covariances are \(0.16+0.03=0.19\) between times 0 and 1, \(0.16+2(0.03)=0.22\) between times 0 and 2, and \(0.16+3(0.03)+2(0.09)=0.43\) between times 1 and 2. Therefore, $$ V_i= \begin{pmatrix} 0.20&0.19&0.22\\ 0.19&0.35&0.43\\ 0.22&0.43&0.68 \end{pmatrix}. $$

The residual variance \(0.04\) has been added to each diagonal entry, not to the off-diagonal entries. For example, the approximate marginal correlation between observations at times 0 and 2 is \(0.22/\sqrt{0.20(0.68)}\approx 0.5966\). The correlation is induced by their shared intercept and slope, even though the residual covariance is diagonal.

From Covariance to an FO Subject Likelihood

Under the FO linearization, the approximate mean vector is \(f_i(\theta,0)\), and the approximate covariance is \(V_i\). If the linearized observation vector is multivariate normal, the individual contribution to the likelihood is its multivariate-normal density: $$ L_i(\theta,\Omega,R_i)= \frac{\exp\!\left[-\frac{1}{2} \{y_i-f_i(\theta,0)\}^{\mathsf T}V_i^{-1} \{y_i-f_i(\theta,0)\}\right]} {(2\pi)^{n_i/2}|V_i|^{1/2}}. $$ The quadratic form measures the observation vector’s discrepancy from the FO mean while accounting for its covariance; the determinant term accounts for the volume or scale of that covariance. Treating each observation as an independent univariate contribution would generally omit the off-diagonal terms and produce a different likelihood.

Worked Example: Evaluating a Two-Observation FO Likelihood

Suppose a hypothetical individual has two observations, each with FO prediction \(10\) concentration units. One random effect has variance \(0.25\) and derivative \(2\) concentration units per random-effect unit for both observations. Assume independent additive Gaussian residual errors with variance \(1\) concentration-unit squared, independent of the random effect. The linearized sum is Gaussian under these assumptions, and $$ J_i\Omega J_i^{\mathsf T} =0.25 \begin{pmatrix}2\\2\end{pmatrix} \begin{pmatrix}2&2\end{pmatrix} = \begin{pmatrix}1&1\\1&1\end{pmatrix}, \qquad V_i= \begin{pmatrix}2&1\\1&2\end{pmatrix}. $$

Let the observed values be \(11\) and \(9\) concentration units. The deviation vector is \(d_i=y_i-f_i(\theta,0)=(1,-1)^{\mathsf T}\) concentration units. Since $$ |V_i|=2(2)-1(1)=3, \qquad V_i^{-1}=\frac{1}{3} \begin{pmatrix}2&-1\\-1&2\end{pmatrix}, $$ the quadratic form is $$ d_i^{\mathsf T}V_i^{-1}d_i =(1,-1)\begin{pmatrix}1\\-1\end{pmatrix}=2. $$ Consequently, the log likelihood contribution is $$ \log L_i=-\frac{2}{2}\log(2\pi)-\frac{1}{2}\log(3)-\frac{1}{2}(2) \approx -3.387. $$ This is the log density for this individual’s two observations under the specified FO Gaussian model. It is not a probability that the model is true.

Assumptions, Diagnostics, and Alternatives

The FO covariance should be read as a local approximation around \(\eta_i=0\). It is most credible when the prediction is adequately represented by its first-order expansion over the random-effect values that materially contribute to the marginal distribution. Strong curvature can make the linearized mean and covariance poor summaries of the nonlinear model. A small random-effect variance can help, but does not by itself guarantee a good approximation; the relevant issue is also how rapidly the predictions change with the random effects.

In practice, check that \(\Omega\), \(R_i\), and the resulting \(V_i\) are valid covariance matrices, and inspect the Jacobian and the induced correlations. An unexpectedly large correlation may be mathematically valid, but it can reveal that observations share a strong random-effect contribution or that the specified random-effect structure deserves scrutiny. A non-positive-definite \(V_i\), or unstable calculations of its inverse and determinant, signals a problem with the covariance specification or numerical computation. Large residual patterns after fitting can also indicate that the FO mean, covariance structure, or structural prediction is inadequate; covariance calculations do not repair structural misspecification.

If residual errors are correlated, their covariance belongs in \(R_i\); it should not be attributed automatically to shared random effects. Conversely, setting \(R_i\) diagonal does not justify setting \(V_i\) diagonal. These are different statements about different sources of variation. The earlier material on distinguishing residual correlation from occasion variability is useful context for keeping those sources conceptually separate.

ApproachWhat it usesKey distinction
FO covariance and likelihoodA Taylor expansion at zero random effects and the induced covariance \(J_i\Omega J_i^{\mathsf T}+R_i\)Efficiently represents shared-effect dependence, but may be inaccurate when curvature matters
Exact marginal likelihood in a linear-Gaussian modelThe specified linear model and Gaussian distributionsMarginalization is exact under those assumptions; FO is not needed to approximate a nonlinear prediction
Numerical integration of the nonlinear modelThe nonlinear prediction across the random-effect distributionCan represent nonlinear marginal behavior more directly, at greater computational cost

A frequent mistake is to reason, “The residual errors are independent, so the observations from this individual are independent.” Independence of the residual errors says only that \(R_i\) has no off-diagonal entries under that residual model. The observations still share \(\eta_i\), so \(J_i\Omega J_i^{\mathsf T}\) can have nonzero off-diagonal entries. The correction is to form the full matrix \(V_i\) before constructing the subject likelihood.

Key takeaway: In FO with multiple observations per individual, each row of the Jacobian describes how one observation responds to the shared random effects. The matrix product \(J_i\Omega J_i^{\mathsf T}\) converts those shared effects into within-individual covariance; add \(R_i\) only after specifying the residual covariance, and use a Gaussian likelihood only when the linearized sum is Gaussian.

Check Your Understanding

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

  1. Two observations have derivative rows \(2\) and \(3\) with respect to a scalar random effect of variance \(0.5\). Residual errors are uncorrelated. What is the random-effect contribution to their covariance?
  2. Does a diagonal residual covariance imply a diagonal FO covariance? Explain.
  3. What condition permits adding \(J_i\Omega J_i^{\mathsf T}\) and \(R_i\) without cross-covariance terms? Is independence required?
  4. For what purpose does the determinant \(|V_i|\) appear in the multivariate-normal log likelihood?

Solutions. 1. The contribution is \(2(0.5)(3)=3\) in the product units of the two observations, because the same random effect enters both. 2. No. A diagonal \(R_i\) has no residual cross-covariances, but \(J_i\Omega J_i^{\mathsf T}\) can have nonzero off-diagonal entries. 3. The cross-covariance between \(J_i\eta_i\) and \(\epsilon_i\) must be zero. Independence is sufficient, not necessary. 4. The determinant term normalizes the multivariate-normal density and accounts for the scale or volume of the covariance; omitting it would not give the stated likelihood.