Constructing the FO Subject Likelihood

From a Linearized Prediction to a Subject Likelihood

A subject’s observations are linked: the same latent random effects contribute to each observation from that subject. Consequently, the likelihood for the subject cannot generally be built by treating those observations as independent after the random effects have been integrated out. The first-order (FO) construction makes this dependence explicit by translating random-effect variation through the locally linear prediction and into an approximate covariance matrix.

The earlier material on the first-order approximation establishes the key local step: approximate the individual prediction by a Taylor expansion around zero random effects. Here, that approximation is assembled into a subject-level likelihood. The central calculation is to combine the covariance induced by the shared random effects with the residual covariance, then evaluate a multivariate normal density.

Definition: For subject $i$, let $y_i$ be the $n_i$-vector of observations, $f_i(\theta,\eta_i)$ the corresponding vector of individual predictions, and $\eta_i$ a $q$-vector of random effects with $\eta_i\sim N(0,\Omega)$. Define $f_{i0}=f_i(\theta,0)$ and the sensitivity matrix $Z_i=\left.\partial f_i(\theta,\eta_i)/\partial\eta_i^\mathsf{T}\right|_{\eta_i=0}$. The FO linearization is $f_i(\theta,\eta_i)\approx f_{i0}+Z_i\eta_i$.

The matrix $Z_i$ has $n_i$ rows and $q$ columns. Its entry in row $j$, column $k$ is the local change in observation $j$’s prediction per unit change in random effect $k$, evaluated at $\eta_i=0$. Let the residual vector be $\epsilon_i$, with mean zero and covariance $R_i$. In the FO construction used here, $R_i$ is evaluated at the zero-random-effect prediction. Assume that subjects are independent, $\eta_i$ is normally distributed with covariance $\Omega$, residuals are normally distributed, and residuals are independent of random effects. The residual covariance may be diagonal or may represent a specified within-subject residual correlation structure.

FO subject likelihood: Under the stated assumptions and linearization, $y_i$ is approximated by a multivariate normal vector with mean $f_{i0}$ and covariance $V_i=Z_i\Omega Z_i^\mathsf{T}+R_i$. Its approximate likelihood contribution is $$ L_i^{\mathrm{FO}}(\theta,\Omega)= \frac{\exp\left\{-\frac12(y_i-f_{i0})^\mathsf{T}V_i^{-1}(y_i-f_{i0})\right\}} {(2\pi)^{n_i/2}|V_i|^{1/2}}. $$ The approximate log likelihood is $$ \ell_i^{\mathrm{FO}}=-\frac12\left[n_i\log(2\pi)+\log|V_i|+(y_i-f_{i0})^\mathsf{T}V_i^{-1}(y_i-f_{i0})\right]. $$

The quadratic form measures the subject’s residual vector relative to its covariance, while $\log|V_i|$ accounts for the volume or scale of that covariance. Both terms matter: a covariance matrix cannot be used only to rescale the residuals while omitting its determinant. With independent subjects, the total FO log likelihood is the sum of the subject contributions, $\ell^{\mathrm{FO}}=\sum_i\ell_i^{\mathrm{FO}}$.

Why Repeated Observations Become Correlated

For observations $j$ and $k$ on the same subject, the random-effect contribution to their approximate covariance is $$ \operatorname{Cov}_{\eta}(y_{ij},y_{ik})\approx z_{ij}\Omega z_{ik}^\mathsf{T}, $$ where $z_{ij}$ and $z_{ik}$ are the corresponding rows of $Z_i$. If residuals are independent across observations, the off-diagonal entry of $R_i$ is zero, but the random-effect contribution need not be. This is the mathematical expression of a shared subject effect: a random effect that shifts one prediction may also shift the others.

More generally, the full covariance entry is $V_{i,jk}=z_{ij}\Omega z_{ik}^\mathsf{T}+R_{i,jk}$. Thus residual correlation, when included in $R_i$, adds to rather than replaces covariance induced by random effects. For a single observation, $Z_i$ is a row vector and the same formula reduces to a scalar variance.

Worked Example: A Shared Random Intercept Across Two Observations

Hypothetical data. Suppose one subject has two observations, each with zero-random-effect prediction 5.0 concentration units. A single additive random intercept has variance $\Omega=0.25$ concentration-units squared, and independent additive residual errors each have variance $0.04$ concentration-units squared. The observed values are 5.1 and 4.9 concentration units. Under the FO linearization, the sensitivity of each prediction to the intercept is 1.

Here $Z_i=(1,1)^\mathsf{T}$, and $R_i$ is diagonal. Therefore $$ V_i=Z_i\Omega Z_i^\mathsf{T}+R_i = \begin{pmatrix}0.25&0.25\\0.25&0.25\end{pmatrix} + \begin{pmatrix}0.04&0\\0&0.04\end{pmatrix} = \begin{pmatrix}0.29&0.25\\0.25&0.29\end{pmatrix}. $$

The residual vector is $r_i=y_i-f_{i0}=(0.1,-0.1)^\mathsf{T}$ concentration units. The determinant is $|V_i|=0.29^2-0.25^2=0.0216$, and $r_i$ lies along the eigenvector $(1,-1)^\mathsf{T}$, whose covariance eigenvalue is $0.29-0.25=0.04$. Hence its quadratic form is $r_i^\mathsf{T}V_i^{-1}r_i=(0.1^2+(-0.1)^2)/0.04=0.5$.

Substitution into the log likelihood gives $$ \ell_i^{\mathrm{FO}}=-\frac12\left[2\log(2\pi)+\log(0.0216)+0.5\right] =-\frac12\left[3.675754-3.835062+0.5\right]\approx -0.1703. $$ The positive off-diagonal covariance, 0.25 concentration-units squared, comes from the common intercept. Ignoring it would describe a different model and produce a different likelihood.

Deriving the Approximate Normal Distribution

The covariance formula follows directly from the linearized model. Write $y_i\approx f_{i0}+Z_i\eta_i+\epsilon_i$. The random vector $Z_i\eta_i+\epsilon_i$ is a sum of independent normal vectors under the assumptions above, so it is normal. Its mean is zero. Its covariance is the sum of the two component covariances:

Propagate random-effect covariance.
$\operatorname{Cov}(Z_i\eta_i)=Z_i\operatorname{Cov}(\eta_i)Z_i^\mathsf{T}=Z_i\Omega Z_i^\mathsf{T}$.

Add residual covariance.
Independence gives $\operatorname{Cov}(Z_i\eta_i+\epsilon_i)=Z_i\Omega Z_i^\mathsf{T}+R_i$, because the cross-covariance terms are zero.

Identify the subject distribution.
Thus $y_i\approx N(f_{i0},V_i)$. Evaluating this multivariate normal density at the observed $y_i$ gives $L_i^{\mathrm{FO}}$.

This derivation is exact for the linearized model, but the linearized model itself is only an approximation to the nonlinear prediction. If the prediction is exactly linear in $\eta_i$ and the residual covariance is fixed at the stated value, the resulting normal distribution and density are exact under the Gaussian assumptions. Otherwise, they define the FO approximate likelihood.

Worked Example: Two Random Effects and Two Observations

Hypothetical data. Suppose a subject has two observations with zero-random-effect predictions 8.0 and 6.0 concentration units, and observed values 8.2 and 5.9. The two random effects are independent with variances 0.16 and 0.09. The sensitivity matrix is $Z_i=\begin{pmatrix}2&1\\1&-1\end{pmatrix}$, and independent residual errors have variance 0.04 at each observation.

The random-effect covariance is $\Omega=\operatorname{diag}(0.16,0.09)$. Multiplying gives $$ Z_i\Omega Z_i^\mathsf{T} = \begin{pmatrix}0.73&0.23\\0.23&0.25\end{pmatrix}, \qquad V_i= \begin{pmatrix}0.77&0.23\\0.23&0.29\end{pmatrix}. $$ For example, the off-diagonal term is $2(1)(0.16)+(1)(-1)(0.09)=0.23$ concentration-units squared.

The residual vector is $r_i=(0.2,-0.1)^\mathsf{T}$ concentration units. The determinant is $0.77(0.29)-0.23^2=0.1704$. For a symmetric matrix with entries $a=0.77$, $b=0.23$, and $c=0.29$, the quadratic form is $$ r_i^\mathsf{T}V_i^{-1}r_i =\frac{c(0.2)^2-2b(0.2)(-0.1)+a(-0.1)^2}{ac-b^2} =\frac{0.0116+0.0092+0.0077}{0.1704} \approx 0.1673. $$

Therefore $$ \ell_i^{\mathrm{FO}}=-\frac12\left[2\log(2\pi)+\log(0.1704)+0.1673\right] \approx -1.0367. $$ Both the random-effect variances and the sensitivity matrix are needed: changing either changes the induced covariance and, consequently, the subject likelihood.

Heterogeneous Residual Variance and a Further Calculation

Residual variance need not be constant across observations. For example, under a proportional residual error model, the variance at a predicted value may increase with the magnitude of that prediction. In the FO construction described here, the residual covariance is evaluated at $f_{i0}$, the prediction at zero random effects. If residuals are independent, this produces a diagonal $R_i$, although its diagonal entries can differ. This evaluation convention is part of the approximation; it does not preserve additional dependence of residual variance on the random effects.

Worked Example: FO Likelihood With Unequal Residual Variances

Hypothetical data. Consider two observations with zero-random-effect predictions 10 and 20 concentration units. A single additive random intercept has variance 0.25 concentration-units squared and sensitivity 1 at both observations. Suppose the residual variance at prediction $f$ is $0.25+(0.10f)^2$, in concentration-units squared. The observations are 11 and 18 concentration units.

At the two baseline predictions, the residual variances are $0.25+(0.10\times10)^2=1.25$ and $0.25+(0.10\times20)^2=4.25$. Thus $$ R_i=\begin{pmatrix}1.25&0\\0&4.25\end{pmatrix}, \qquad V_i=\begin{pmatrix}1.50&0.25\\0.25&4.50\end{pmatrix}. $$ The off-diagonal value 0.25 is induced by the shared random intercept; the unequal diagonal residual variances reflect the proportional component.

The residual vector is $r_i=(1,-2)^\mathsf{T}$ concentration units. The determinant is $1.50(4.50)-0.25^2=6.6875$. Its quadratic form is $$ r_i^\mathsf{T}V_i^{-1}r_i =\frac{4.50(1)^2-2(0.25)(1)(-2)+1.50(-2)^2}{6.6875} =\frac{11.5}{6.6875} \approx 1.7196. $$

The log likelihood is therefore $$ \ell_i^{\mathrm{FO}}=-\frac12\left[2\log(2\pi)+\log(6.6875)+1.7196\right] \approx -3.6478. $$ The residual variances were evaluated at predictions 10 and 20, not at the observed values 11 and 18. Substituting observed values instead would use a different variance convention from the one specified for this FO calculation.

Assumptions, Checks, and Alternatives

Before using an FO subject likelihood, check that the model components match the calculation. The random effects must have the stated mean and covariance, the residual covariance must represent the assumed residual model at the chosen evaluation point, and the random effects and residuals must be independent for the simple covariance sum above. The normal density also requires $V_i$ to be positive definite, so it has an inverse and a positive determinant. With a positive-definite residual covariance and a positive-semidefinite $\Omega$, $V_i$ is positive definite.

Useful computational checks include verifying that $Z_i$ has dimensions $n_i\times q$, $\Omega$ has dimensions $q\times q$, and $R_i$ and $V_i$ have dimensions $n_i\times n_i$. Check that $V_i$ is symmetric, that its diagonal entries are positive, and that its determinant is positive. For a small example, compute $Z_i\Omega Z_i^\mathsf{T}$ directly and confirm that the off-diagonal terms agree with the row-wise covariance formula.

A common mistake is to calculate a separate marginal variance for each observation and then multiply independent normal densities. That discards off-diagonal covariance from shared random effects. The correction is to form the full $V_i$, retain its determinant, and use the full quadratic form. Another mistake is to omit $\log|V_i|$: doing so fails to account for the density’s normalization and changes how parameter values are evaluated.

The exact marginal likelihood instead integrates the conditional density over the random-effect distribution without replacing the nonlinear prediction by its local linearization. It is generally an integral that may require numerical evaluation. In the linear-Gaussian case, the covariance construction above gives the exact marginal density; with nonlinear predictions, it gives the FO approximation. FO is computationally direct, but its quality depends on whether the local linearization and the chosen residual covariance evaluation adequately represent the model over the random-effect distribution. The likelihood formula itself does not diagnose that adequacy.

In practice, inspect the implied covariance and its components, not just the final scalar objective. Implausible correlations, a non-positive-definite matrix, or a nearly singular covariance can indicate incompatible variance assumptions, a coding or indexing error, or poorly identified components. A successful matrix calculation does not establish that the FO approximation is accurate; it establishes only that the specified approximate likelihood was evaluated consistently.

Key takeaway: For subject $i$, the FO construction uses $V_i=Z_i\Omega Z_i^\mathsf{T}+R_i$ and evaluates the multivariate normal density with mean $f_{i0}$. Repeated observations share random effects, so their likelihood contribution generally contains off-diagonal covariance.

Check Your Understanding

Answer each question, then compare with the solutions below.

If a subject has three observations and two random effects, what are the dimensions of $Z_i$, $\Omega$, and $V_i$?
With independent residual errors, why can two observations from one subject still have nonzero covariance under FO?
For one observation, suppose the sensitivity to a scalar random effect is 2, its variance is 0.09, and the residual variance is 0.16. What is the FO variance?
What two terms appear in the FO log likelihood in addition to $n_i\log(2\pi)$, and what role does each play?
Under what condition does the FO likelihood coincide with the exact marginal likelihood under the stated Gaussian assumptions?

Solutions. (1) $Z_i$ is $3\times2$, $\Omega$ is $2\times2$, and $V_i$ is $3\times3$. (2) The same random effects enter both predictions, contributing $z_{ij}\Omega z_{ik}^\mathsf{T}$, even if the residual covariance off the diagonal is zero. (3) The variance is $2^2(0.09)+0.16=0.52$, in the squared units of the observation. (4) The terms are $\log|V_i|$, which accounts for covariance scale in the density normalization, and the quadratic form, which measures the residual vector relative to that covariance. (5) It is exact when the prediction is linear in the random effects and the residual covariance is fixed as specified, with the Gaussian and independence assumptions holding.