From a Nonlinear Prediction to a Gaussian Approximation
In a nonlinear mixed-effects model, an individual’s prediction depends on random effects. The marginal likelihood for that individual averages the conditional observation model over the random-effect distribution. For a nonlinear prediction, this integration is generally difficult. The first-order (FO) approximation makes the prediction locally linear in the random effects, using its slope at zero, so that the resulting Gaussian model can be integrated analytically.
This approximation is local: it uses the prediction and its derivatives at the population value of the random effects. It does not replace the population model with an arbitrary linear model, nor does it simply set the random effects to zero. The local slopes determine how random-effect variation contributes to the approximate observation covariance.
The assumptions for the integration step are explicit: the random effects have a multivariate normal distribution with mean zero and covariance \(\Omega\); conditional residuals are Gaussian with covariance \(R_i\); and residuals are independent of random effects. In this formulation, \(R_i\) does not depend on \(\eta_i\). The prediction itself may be nonlinear, but it must be differentiable with respect to the random effects at zero for the stated linearization to be defined.
These assumptions describe the model being approximated. The FO result is exact if the prediction is already affine in \(\eta_i\). With a nonlinear prediction, the Taylor expansion is an approximation, and its quality depends on how well the tangent at zero represents predictions over the random-effect values that matter.
Linearizing at Zero Random Effects
Let \(m_i=f_i(\theta,0)\) be the prediction at zero random effects. Let \(J_i\) be the \(n_i\times q\) Jacobian matrix of prediction derivatives with respect to the random effects, evaluated at zero:
The first-order Taylor expansion is
Thus the approximate observation model is \(y_i\approx m_i+J_i\eta_i+\epsilon_i\). The random-effect contribution is not discarded: it enters through \(J_i\eta_i\). Because this expression is a sum of independent Gaussian vectors under the stated assumptions, its marginal distribution is Gaussian.
Why the Integration Gives This Covariance
After linearization, define \(z_i=J_i\eta_i+\epsilon_i\), so that \(y_i\approx m_i+z_i\). Its mean is zero because both \(\eta_i\) and \(\epsilon_i\) have mean zero. Independence gives
To verify that the result is Gaussian, consider the characteristic function of \(y_i\), evaluated at any \(n_i\)-vector \(t\). Independence makes the characteristic function of the sum a product:
This is the characteristic function of a multivariate normal distribution with mean \(m_i\) and covariance \(V_i=R_i+J_i\Omega J_i^\mathsf{T}\), establishing the stated marginal result for the linearized model. Equivalently, it is the integral of the conditional Gaussian observation density against the Gaussian random-effect density. In a nonlinear model, FO approximates that integral by first changing the prediction inside the integral; it does not make the original nonlinear integral exact.
Worked Example: One Observation and One Random Effect
Suppose a hypothetical prediction for one concentration is \(f(\eta)=10\exp(\eta)\), where \(\eta\sim\mathcal N(0,0.04)\). The additive residual variance is \(4\) concentration-units squared, and the residual is independent of \(\eta\). Find the FO marginal distribution and the approximate density at an observed concentration of \(13\).
At zero random effect, the prediction is \(m=f(0)=10\). Its derivative is \(f'(\eta)=10\exp(\eta)\), so the FO slope is \(J=f'(0)=10\). The random-effect contribution to variance is \(J^2\Omega=10^2(0.04)=4\) concentration-units squared. Therefore
For \(y=13\), the approximate density is
This is an approximate density per concentration unit, not a probability that the observation equals exactly \(13\). The calculation includes residual variance and the extra variance propagated through the local slope.
Worked Example: Two Observations and Correlated Random Effects
Consider a hypothetical individual with two observations, both in the same measurement units. Let the predictions be \(f_1(\eta)=10\exp(\eta_1)\) and \(f_2(\eta)=20+4\eta_2\), where \(\Omega=\begin{pmatrix}0.04&0.01\\0.01&0.09\end{pmatrix}\). Assume independent residuals with covariance \(R=\begin{pmatrix}1&0\\0&4\end{pmatrix}\), independent of the random effects. Calculate the FO mean and covariance.
At \(\eta=0\), the prediction vector is \(m=(10,20)^\mathsf{T}\). The derivative of the first prediction with respect to \(\eta_1\) is \(10\) at zero, and the second prediction has slope \(4\) with respect to \(\eta_2\). Hence
For example, the off-diagonal entry is \(10(0.01)(4)=0.4\); it carries the covariance between the two random effects through the two slopes. Adding the residual covariance gives
The FO approximation is therefore a bivariate normal distribution with mean \((10,20)^\mathsf{T}\) and covariance \(V\). The determinant is \(5(5.44)-0.4^2=27.04>0\), consistent with a positive-definite covariance matrix. The nonzero covariance between observations comes from correlated random effects, not from residual correlation in this example.
Worked Example: A Case Where FO Is Exact
Suppose a hypothetical scalar prediction is \(f(\eta)=8+2\eta\), with \(\eta\sim\mathcal N(0,0.25)\) and independent additive residual variance \(1\). Find the marginal distribution and determine whether FO introduces approximation error.
Here \(f(0)=8\) and \(f'(0)=2\). Since the prediction is affine, \(f(\eta)=f(0)+f'(0)\eta\) for every \(\eta\), not just near zero. The marginal mean and variance are
Thus \(Y\sim\mathcal N(8,2)\) exactly. There is no Taylor remainder because all second and higher derivatives of this prediction are zero. This example separates the accuracy of the Gaussian integration itself, which is exact for the linear-Gaussian model, from the linearization error that may occur for nonlinear predictions.
Assumptions, Diagnostics, and Alternatives
FO is most credible when the prediction is close to linear over the random-effect values with appreciable probability. A direct diagnostic is to compare the nonlinear prediction with its tangent at zero across representative random-effect values, for example values spanning the central range implied by \(\Omega\). Large discrepancies indicate that the local approximation may be poor. A small derivative in one direction also deserves attention: FO assigns little observation-scale variance to that direction, even if curvature makes the nonlinear prediction change substantially away from zero.
The covariance calculation provides additional checks. \(J_i\Omega J_i^\mathsf{T}\) must be positive semidefinite when \(\Omega\) is positive semidefinite; it cannot subtract variance. If the random-effect covariance is diagonal, the resulting contribution may still have off-diagonal entries when observations depend on shared random effects. Conversely, omitting those entries when they are present changes the implied marginal model. A non-positive-definite \(V_i\), when the residual covariance is positive definite, is a sign of an implementation or parameterization problem.
FO depends on the specified Gaussian assumptions as well as on the Taylor expansion. If residual variance depends on random effects, the conditional covariance also varies across \(\eta_i\); the formula \(R_i+J_i\Omega J_i^\mathsf{T}\) does not account for that dependence. Non-Gaussian random effects or residuals likewise do not generally produce the stated Gaussian marginal distribution. These are model changes, not merely numerical details.
| Approach | How it handles marginalization | Main consideration |
|---|---|---|
| FO | Linearizes predictions at zero random effects, then integrates the Gaussian approximation. | Computationally direct, but local linearization can be inaccurate for strongly nonlinear predictions. |
| Exact integration | Integrates the original conditional model over the random-effect distribution without linearizing. | Available in special cases, such as affine predictions with Gaussian effects; nonlinear integrals may not have a closed form. |
| Numerical integration | Approximates the original integral by evaluating it at selected random-effect values. | Can represent nonlinear behavior more directly, but computational cost and numerical accuracy depend on dimension and integration strategy. |
A common mistake is to evaluate the prediction at \(\eta_i=0\), use only the residual covariance \(R_i\), and call the result the FO marginal model. That calculation omits random-effect variability. The correction is to retain the local slope and add \(J_i\Omega J_i^\mathsf{T}\) to the residual covariance. A different mistake is to assume that this added covariance corrects every consequence of nonlinearity: it is the covariance of the linearized model, not a guarantee that the original nonlinear marginal distribution is Gaussian or has the same mean and covariance.
Check Your Understanding
Try each question, then check the solutions below.
- A scalar prediction is \(f(\eta)=7\exp(\eta)\), with random-effect variance \(0.01\) and residual variance \(0.36\). What are the FO mean and variance?
- Why can FO produce a covariance between two observations when their residual errors are independent?
- Under what prediction form is FO exact, given Gaussian random effects and independent Gaussian residuals?
- If \(f(\eta)=12+\eta^2\), what does the FO derivative at zero imply, and what limitation does this illustrate?
Solutions. (1) \(f(0)=7\) and \(f'(0)=7\), so the FO mean is \(7\) and the variance is \(0.36+7^2(0.01)=0.85\). (2) The term \(J\Omega J^\mathsf{T}\) can have off-diagonal entries when observations depend on correlated or shared random effects; residual independence only makes \(R\) diagonal in this setup. (3) FO is exact when the prediction is affine in the random effects, because its first-order Taylor expansion then has no remainder and the Gaussian sum is exactly Gaussian. (4) The derivative at zero is \(0\), so FO adds no random-effect variance for this prediction and approximates it by the constant \(12\). Yet \(12+\eta^2\) changes away from zero. This illustrates that a zero local slope does not establish that the nonlinear prediction is constant over the random-effect distribution.