First-Order Estimation by Linearization

Why Linearize Individual Predictions?

In a population model, each subject’s prediction depends on population parameters and that subject’s random effects. The likelihood contribution for a subject is therefore obtained by averaging the conditional data density over the random-effects distribution. When predictions are nonlinear in the random effects, this integral will generally not have a closed form. First-order (FO) estimation makes the integral tractable by approximating the predictions with a linear function of the random effects near zero.

The approximation is local: it uses the prediction and its slope at the population-typical random effect, not the full shape of the nonlinear prediction over the random-effects distribution. That distinction matters when the prediction has substantial curvature or the random effects are widely dispersed.

Definition: For subject $i$, let $y_i$ be an $m_i$-vector of observations, $\eta_i$ a $q$-vector of random effects, and $f_i(\theta,\eta_i)$ the $m_i$-vector of model predictions. The first-order approximation expands $f_i$ in a Taylor series around $\eta_i=0$ and retains the constant and linear terms: $$ f_i(\theta,\eta_i)\approx f_i(\theta,0)+J_i(\theta)\eta_i, \qquad J_i(\theta)=\left.\frac{\partial f_i(\theta,\eta)}{\partial\eta^\mathsf{T}}\right|_{\eta=0}. $$ Here $J_i$ is an $m_i\times q$ sensitivity matrix.

From the Linearization to a Marginal Likelihood

Assume the random effects have distribution $\eta_i\sim N(0,\Omega)$, where $\Omega$ is a positive semidefinite covariance matrix. Assume also that the residual error vector is independent of $\eta_i$ and follows $\epsilon_i\sim N(0,R_i)$, with $R_i$ fixed with respect to $\eta_i$. The observation model is $y_i=f_i(\theta,\eta_i)+\epsilon_i$. For the basic FO derivation, $R_i$ may depend on model parameters but does not vary with the random effects.

After replacing the nonlinear prediction with its linearization, the approximate observation model is $y_i\approx f_i(\theta,0)+J_i\eta_i+\epsilon_i$. A linear transformation of a normal random vector is normal, and the sum of independent normal vectors is normal. Thus the approximate marginal distribution has mean $f_i(\theta,0)$ and covariance $V_i=J_i\Omega J_i^\mathsf{T}+R_i$. The random effects contribute covariance through their effect on predictions; they do not shift the approximate mean away from the zero-random-effect prediction.

y_i\ \dot{\sim}\ N\!\left(f_i(\theta,0),V_i\right), \qquad V_i=J_i\Omega J_i^\mathsf{T}+R_i.

The notation $\dot{\sim}$ emphasizes that this is the distribution implied by the linearized model, not necessarily the exact marginal distribution of the original nonlinear model. For $N$ independent subjects, the resulting FO log-likelihood is the sum of the subject contributions:

\ell_{\mathrm{FO}}(\theta,\Omega,\ldots) = -\frac{1}{2}\sum_{i=1}^{N} \left[ m_i\log(2\pi)+\log|V_i| +\{y_i-f_i(\theta,0)\}^{\mathsf T}V_i^{-1} \{y_i-f_i(\theta,0)\} \right].

This expression requires $V_i$ to be positive definite for each subject’s contribution to have an ordinary nonsingular Gaussian density. The quadratic form is dimensionless. The determinant and density do depend on measurement units: $|V_i|$ has units equal to the product of the outcome units squared across the $m_i$ observations, and the density has reciprocal outcome-volume units. If outcomes are converted by a fixed change of units, the log-likelihood changes by an additive Jacobian constant. That constant does not change parameter estimates when the conversion is fixed, but the units do not cancel between the determinant and quadratic-form terms.

A Short Derivation and Its Implications

The Taylor expansion provides the approximation directly. For a small vector $\eta_i$, the prediction is its value at zero plus its first derivative at zero times $\eta_i$, with terms involving second and higher derivatives omitted. Substituting that expression into the observation model gives a normal random vector because both $\eta_i$ and $\epsilon_i$ are normal and the approximation is linear in $\eta_i$. Its mean and covariance follow from the mean and covariance rules:

\begin{aligned} E[y_i]_{\mathrm{FO}} &= f_i(\theta,0)+J_iE[\eta_i]+E[\epsilon_i] = f_i(\theta,0),\\ \operatorname{Var}(y_i)_{\mathrm{FO}} &=J_i\operatorname{Var}(\eta_i)J_i^\mathsf{T} +\operatorname{Var}(\epsilon_i) =J_i\Omega J_i^\mathsf{T}+R_i. \end{aligned}

This derivation also explains two practical features. First, subjects with different observation times or covariates can have different sensitivity matrices and therefore different contributions from the same $\Omega$. Second, within-subject observations can be correlated in the marginal model even when residual errors are independent: a shared random effect influences more than one prediction.

Worked Example: One Observation from an Exponential Random-Effects Model

Consider the hypothetical model $f(\eta)=\theta\exp(\eta)$, with $\eta\sim N(0,\omega^2)$ and independent residual error $\epsilon\sim N(0,\sigma^2)$. Take $\theta=14.62$ outcome units, $\omega=0.20$, $\sigma=1.00$ outcome unit, and an observed value $y=15.00$ outcome units.

At zero random effect, $f(0)=14.62$. The derivative is $f'(\eta)=\theta\exp(\eta)$, so $J=f'(0)=14.62$ outcome units per unit of $\eta$. The FO prediction is therefore $14.62+14.62\eta$, and the approximate marginal variance is

V=(14.62)^2(0.20)^2+(1.00)^2 =213.7444(0.04)+1 =9.5498\ \text{outcome units}^2.

The residual from the zero-random-effect prediction is $r=15.00-14.62=0.38$ outcome units. The approximate log-likelihood contribution is

\begin{aligned} \ell_{\mathrm{FO}} &=-\frac12\left[\log(2\pi)+\log(9.5498)+\frac{(0.38)^2}{9.5498}\right]\\ &\approx-\frac12[1.8379+2.2565+0.0151] \approx-2.0548. \end{aligned}

This is the likelihood contribution under the FO Gaussian approximation. It is not the exact marginal likelihood for the exponential model: the exponential prediction is curved in $\eta$, and linearization omits that curvature.

Worked Example: Shared Random Effect Creates Within-Subject Covariance

For a hypothetical subject, suppose two predictions are $f_1(\eta)=10\exp(\eta)$ and $f_2(\eta)=8\exp(\eta)$, with a shared $\eta\sim N(0,0.10^2)$. Let the two residual errors be independent, each with variance $0.25$ outcome units squared. The observed values are $y=(11.0,8.2)^\mathsf{T}$ outcome units.

At zero, the prediction vector is $(10,8)^\mathsf{T}$, and the sensitivity matrix is the column vector $J=(10,8)^\mathsf{T}$. Thus the marginal covariance is

V=0.01 \begin{pmatrix} 100&80\\ 80&64 \end{pmatrix} + \begin{pmatrix} 0.25&0\\ 0&0.25 \end{pmatrix} = \begin{pmatrix} 1.25&0.80\\ 0.80&0.89 \end{pmatrix}.

The off-diagonal covariance, $0.80$ outcome units squared, arises from the common random effect. The residual covariance alone is diagonal, but that does not make the marginal covariance diagonal. The residual vector is $r=(1.0,0.2)^\mathsf{T}$ outcome units. The determinant is $1.25(0.89)-(0.80)^2=0.4725$, and the quadratic form is

r^\mathsf{T}V^{-1}r = \frac{0.89(1.0)^2-2(0.80)(1.0)(0.2)+1.25(0.2)^2}{0.4725} = \frac{0.6200}{0.4725} \approx1.3122.

Therefore,

\ell_{\mathrm{FO}} =-\frac12\left[2\log(2\pi)+\log(0.4725)+1.3122\right] \approx-2.119.

The covariance structure matters to the likelihood: treating these observations as independent marginally would discard the effect of the shared random effect and yield a different quadratic form and determinant.

Curvature, Diagnostics, and Failure Modes

A first-order approximation is most credible when prediction changes are close to linear over the random-effects values with appreciable probability. A useful local check is to inspect second derivatives or compare predictions at zero with predictions at plausible positive and negative random-effect values. For a scalar random effect with standard deviation $\omega$, values such as $-2\omega$, $0$, and $2\omega$ give a simple check of the region containing much of the random-effects distribution. This is a diagnostic, not a guarantee: the joint distribution and the sensitivity of the likelihood also matter.

The approximation can be poor when random-effect variability is large relative to the scale over which the prediction changes, when the zero point is near a boundary or steep transition, or when nonlinear transformations produce substantial curvature. A first-order model also fixes the approximate marginal mean at $f_i(\theta,0)$. For a curved function, the exact expectation $E[f_i(\theta,\eta_i)]$ need not equal that value. This discrepancy can affect estimation, not merely the shape of the distribution.

Worked Example: Linearizing a Saturable Prediction

Consider the hypothetical prediction $f(\eta)=100\exp(\eta)/(5+\exp(\eta))$, with $\eta\sim N(0,0.30^2)$. Let the residual standard deviation be $2.00$ outcome units and the observation be $20.00$ outcome units. At zero, the prediction is $100/6=16.6667$ outcome units. Differentiation gives

f'(\eta)=\frac{100(5)\exp(\eta)}{(5+\exp(\eta))^2}, \qquad J=f'(0)=\frac{500}{36}=13.8889.

The FO variance and residual are

V=(13.8889)^2(0.30)^2+(2.00)^2 \approx17.3611+4.0000=21.3611\ \text{outcome units}^2, \qquad r=20.00-16.6667=3.3333\ \text{outcome units}.

The approximate log-likelihood contribution is

\ell_{\mathrm{FO}} =-\frac12\left[\log(2\pi)+\log(21.3611)+\frac{(3.3333)^2}{21.3611}\right] \approx-2.710.

Here the derivative describes the local sensitivity at $\eta=0$. The prediction is saturable rather than linear, so this calculation alone cannot establish that the approximation is adequate across the spread of $\eta$. Checking predictions over plausible random-effect values, or comparing with a more direct marginal-likelihood calculation, is warranted if curvature is consequential.

Other warning signs include estimates that change materially when starting values or model parameterizations change, an approximate covariance matrix that is nearly singular, and poor agreement between FO likelihood evaluations and evaluations using numerical integration. These signs do not prove that linearization is the cause, but they motivate checking the approximation and the model implementation. Also verify that the derivative matrix corresponds to the exact prediction used by the model and that $R_i$ and $\Omega$ have the intended units and structure.

Choosing and Interpreting the Approximation

FO is attractive because it replaces a nonlinear random-effects integral with a multivariate normal density whose mean and covariance are explicit. Direct numerical integration or quadrature instead evaluates the nonlinear conditional density across random-effect values, avoiding the particular Taylor approximation but requiring suitable numerical integration. The later material on FOCE conditional likelihood treats expansion around a subject-specific conditional value rather than zero; it is a distinct approximation, not a synonym for FO. Laplace approximation and stochastic estimation are further alternatives, each with their own assumptions and computational trade-offs.

FO is not automatically preferable simply because it is computationally convenient, nor is it automatically invalid because the structural model is nonlinear. Its adequacy depends on how nonlinear the predictions are over the random-effects distribution and how that approximation interacts with the data and model parameters. A defensible workflow records the model, the expansion point, the random-effects and residual assumptions, and diagnostics or sensitivity checks supporting the approximation.

Common pitfall: A common error is to interpret $f_i(\theta,0)$ as the exact marginal mean, or to assume independent observations remain marginally independent when they share random effects. Under FO, $f_i(\theta,0)$ is the approximate mean, and $J_i\Omega J_i^\mathsf{T}$ can induce covariance between observations. Both statements follow from the approximation and assumptions; neither is a general property of every nonlinear mixed-effects model.

Check Your Understanding

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

For a scalar prediction $f(\eta)$, what two quantities from the expansion at zero determine its FO approximation?
If a subject has two observations and one shared random effect, can the approximate marginal covariance be nondiagonal when residual errors are independent? Explain.
For $f(\eta)=\theta\exp(\eta)$, give the FO mean and variance when $\eta\sim N(0,\omega^2)$ and residual variance is $\sigma^2$.
Why does a fixed conversion from one outcome unit to another change the log-likelihood by a constant rather than changing the parameter estimates?
Name one diagnostic that can reveal a potentially poor FO approximation.

Solutions. (1) The value $f(0)$ and derivative $f'(0)$; together they give $f(\eta)\approx f(0)+f'(0)\eta$. (2) Yes. The shared random effect contributes $J\Omega J^\mathsf{T}$, which can have nonzero off-diagonal entries even if $R_i$ is diagonal. (3) The approximate mean is $\theta$, and the approximate variance is $\theta^2\omega^2+\sigma^2$, because $f(0)=\theta$ and $f'(0)=\theta$. (4) The density transforms by a fixed Jacobian factor under a fixed unit conversion, so its log adds a parameter-independent constant; maximizing it gives the same parameter estimates. (5) One can compare predictions at zero with predictions at plausible random-effect values, or compare FO likelihood evaluations with a direct numerical integration when available.