Why FO Bias Grows With Random-Effect Variance

Why the First-Order Approximation Can Distort Estimates

The first-order (FO) approximation replaces an individual prediction that depends nonlinearly on random effects with its local linear approximation. This can be useful, but its accuracy depends on more than whether the model is differentiable: it also depends on how much random-effect variation the model explores and how curved the prediction is over that range. When appreciable curvature is ignored, the population distribution implied by the linearized model can differ from the distribution implied by the original model.

The central issue is that a centered random effect does not make a nonlinear prediction centered at its value when the random effect is zero. Curvature can shift the population-average prediction, and it can add variability beyond the contribution captured by the local slope. A likelihood fit may then compensate by changing fixed population parameters, random-effect variances, or both. The direction and size of those changes depend on the model, observation design, and error structure; there is no universal rule that FO always raises or lowers a particular estimate.

Definition: Let an individual prediction be $f(\theta,\eta)$, where $\theta$ is a population parameter vector and $\eta$ is a centered random-effect vector with covariance matrix $\Omega$. The FO approximation about $\eta=0$ keeps the constant and linear terms: $$ f(\theta,\eta)\approx f(\theta,0)+J(\theta)\eta, $$ where $J(\theta)$ is the derivative of $f$ with respect to $\eta$, evaluated at zero. The terms omitted include quadratic and higher-order terms. Curvature refers, locally, to the second derivatives of $f$ with respect to $\eta$.

What Curvature Adds to the Population Distribution

For a scalar random effect, suppose $\eta\sim N(0,\omega^2)$, and consider a second-order expansion of a scalar prediction:

f(\eta)\approx f_0+g\eta+\frac{1}{2}h\eta^2,

where $f_0=f(0)$, $g$ is the slope at zero, and $h$ is the second derivative at zero. The linear FO prediction has expectation $f_0$, since $E(\eta)=0$. The second-order expansion instead has expectation $f_0+\tfrac{1}{2}h\omega^2$, since $E(\eta^2)=\omega^2$. Thus a centered random effect can produce a nonzero shift in the average prediction whenever curvature and random-effect variance are both nonzero.

The variance also changes. For a normal random effect, $E(\eta^3)=0$ and $\operatorname{Var}(\eta^2)=2\omega^4$. The covariance between the linear and quadratic terms is therefore zero, giving:

\operatorname{Var}\!\left(g\eta+\frac{1}{2}h\eta^2\right) =g^2\omega^2+\frac{1}{2}h^2\omega^4.

The FO variance contribution is $g^2\omega^2$. The second term, $\tfrac{1}{2}h^2\omega^4$, is absent from that linearized prediction. It grows with the square of the random-effect variance and with the square of the curvature. This explains why increasing random-effect variability can make the neglected terms important even if the local linear approximation seems reasonable very close to zero.

These calculations describe the prediction under the stated second-order polynomial. For a general nonlinear function, a Taylor expansion also has higher-order remainder terms; the displayed corrections are useful when the expansion is accurate over the relevant range of $\eta$. They are not a guarantee that the whole nonlinear model is summarized by its second-order expansion.

A Justification of the Multivariate Mean Correction

The same mechanism applies with multiple random effects. Let $H$ be the Hessian matrix of second derivatives of a scalar prediction with respect to $\eta$, evaluated at zero. The quadratic term in the Taylor expansion is $\tfrac{1}{2}\eta^\mathsf{T}H\eta$. For a centered random-effect vector with covariance $\Omega$, its expectation is $\tfrac{1}{2}\operatorname{tr}(H\Omega)$. To see this, write the quadratic form as a sum over its entries:

E(\eta^\mathsf{T}H\eta) =\sum_{a,b}H_{ab}E(\eta_a\eta_b) =\sum_{a,b}H_{ab}\Omega_{ab} =\operatorname{tr}(H\Omega),

where $H$ and $\Omega$ are symmetric. Hence the second-order correction to the mean is $\tfrac{1}{2}\operatorname{tr}(H\Omega)$. Curvature aligned with directions of high random-effect variance can make this correction substantial. Correlations in $\Omega$ can also matter because they contribute to the cross-products in the trace.

Worked Example: Quantifying the Omitted Terms

Consider a hypothetical prediction, in concentration units, given by $f(\eta)=10+2\eta+3\eta^2$, where $\eta\sim N(0,0.25)$ is dimensionless. The slope at zero is $g=2$, and the second derivative is $h=6$.

FO mean and variance. The FO prediction is $10+2\eta$. Its mean is $10$, and its random-effect variance is: $$ g^2\omega^2=2^2(0.25)=1.00 $$ in squared concentration units.

Second-order mean. The quadratic contribution to the mean is: $$ \frac{1}{2}h\omega^2=\frac{1}{2}(6)(0.25)=0.75. $$ Thus the polynomial prediction has mean $10+0.75=10.75$ concentration units, not $10$.

Second-order variance. The variance contribution from the quadratic term is: $$ \frac{1}{2}h^2\omega^4 =\frac{1}{2}(6^2)(0.25^2) =18(0.0625)=1.125. $$ The total random-effect variance is $1.00+1.125=2.125$ squared concentration units. In this example, FO omits a mean shift of $0.75$ concentration units and a variance contribution of $1.125$ squared concentration units. The omitted variance is larger than the FO variance contribution, so a local linearization is not a close summary of this prediction distribution.

How Approximation Error Becomes Parameter Bias

Approximation error is not itself the same thing as parameter bias. Bias arises when parameters are estimated using an approximate likelihood that does not represent the data distribution well enough. The fitting procedure can compensate for a missing mean shift by moving a population parameter, and it can compensate for missing variability by altering a variance parameter or other features of the model. Which compensation is favored depends on the observations and on which parameters can change the predicted distribution.

The following example isolates a mean shift. It assumes that the observations are informative about their population mean and that the only freely adjusted population parameter in the FO fit is an intercept. This is a moment-matching illustration, not a claim about every likelihood fit.

Worked Example: A Population Intercept Compensates for Curvature

Suppose a hypothetical model has prediction $f(\theta,\eta)=\theta+2\eta+3\eta^2$, with $\eta\sim N(0,0.25)$, and the true population intercept is $\theta=10$ concentration units. The mean of the nonlinear prediction is: $$ E\{f(\theta,\eta)\} =\theta+2E(\eta)+3E(\eta^2) =10+0+3(0.25)=10.75 $$ concentration units.

The FO model drops $3\eta^2$, so its population mean is simply $\theta_{\mathrm{FO}}$. If the fit matches the population mean, it must set $\theta_{\mathrm{FO}}=10.75$ concentration units. In this simplified setting, the estimated intercept is $0.75$ concentration units above the true intercept. This shift is a consequence of omitted curvature, not evidence that the true intercept is $10.75$.

In a full model, other parameters may also change the predicted mean, and the fit uses the full likelihood rather than only this moment. Consequently, the same curvature does not imply that the intercept will always be biased upward. The example shows a possible compensation mechanism and why its direction must be assessed in the model context.

How Variance Estimates Can Be Distorted

The omitted quadratic variability can also be absorbed by a random-effect variance estimate. Consider a case where the residual variance is known, the slope $g$ is fixed, and the fit matches the remaining prediction variance using the FO expression $g^2\omega^2$. Under those restrictive conditions, the FO variance parameter must account for variability that the true curved prediction attributes partly to its quadratic term.

Worked Example: Moment Matching Inflates a Variance Estimate

Use the prediction in the first example, $f(\eta)=10+2\eta+3\eta^2$, with true $\omega^2=0.25$. Suppose independent residual error has known variance $0.50$ squared concentration units. The total variance under the second-order prediction is: $$ 0.50+2.125=2.625 $$ squared concentration units.

The FO model has total variance $0.50+4\omega_{\mathrm{FO}}^2$. Matching it to $2.625$ gives: $$ 0.50+4\omega_{\mathrm{FO}}^2=2.625,\qquad \omega_{\mathrm{FO}}^2=\frac{2.625-0.50}{4}=0.53125. $$

The moment-matched FO variance is $0.53125$, compared with the true $0.25$, a difference of $0.28125$ in the dimensionless variance parameter. The FO fit inflates the variance here because it has no quadratic term with which to represent the extra variability. This direction follows from the example's fixed slope, known residual variance, and variance-matching setup. If those conditions change, or if other parameters can compensate, the direction of the estimated variance change need not be the same.

When FO Is Exact, and What to Compare It With

If the prediction is affine in $\eta$, the random effects are Gaussian, and the residual model is additive Gaussian with constant covariance independent of $\eta$, integrating the linearized model gives a Gaussian marginal distribution with mean $f(\theta,0)$ and covariance $R+J\Omega J^\mathsf{T}$, where $R$ is the residual covariance. Under those conditions the linearized model is exact for the specified model: there are no omitted curvature terms. This is the linear-Gaussian benchmark. The approximation becomes a concern when the original prediction is materially nonlinear over the random-effect distribution, or when other features of the conditional distribution also vary with $\eta$.

The main comparison is between fitting the FO-approximated likelihood and evaluating the marginal likelihood implied by the nonlinear model more directly, for example by numerical integration when feasible. Direct integration can represent the curvature that FO omits, but it may be computationally demanding, especially with multiple random effects. FO is simpler, but its computational convenience does not remove the need to check whether its approximation is adequate for the model and data. A linear-Gaussian benchmark is also useful: agreement with that benchmark is expected when the model really satisfies its assumptions, not merely because a model has been linearized.

Diagnostics, Failure Modes, and Limits

FO bias cannot be diagnosed from a single estimate or from a large standard error alone. Look for signs that the approximation is consequential and test sensitivity to the approximation where possible.

Assess curvature over plausible random effects. Compare the nonlinear prediction with its linear expansion across values of $\eta$ that have appreciable probability under the fitted $\Omega$. Large differences indicate that the local expansion may not describe much of the random-effect distribution.
Check residual patterns in context. Systematic residual structure across time, predicted concentrations, or other relevant covariates can signal inadequacy, although such patterns do not identify FO error uniquely; structural or residual-model misspecification can produce similar evidence.
Examine stability to random-effect variability. If estimates or the objective change substantially as plausible values of $\Omega$ are explored, the terms ignored by FO may be important. Stability is reassuring but does not prove that the approximate likelihood is correct.
Use simulation recovery when possible. Simulate data from a specified nonlinear model with known parameters and refit using FO. Repeated systematic differences between estimates and generating values can reveal bias for that design and model. This does not establish the bias in every application.
Compare with a more direct likelihood evaluation when feasible. Material changes in estimates can indicate approximation sensitivity. The comparison is informative only if both fits are numerically reliable and use otherwise comparable models.

A common mistake is to argue that FO must be accurate because $E(\eta)=0$, so the random effects “average out.” That reasoning ignores $E(\eta^2)=\omega^2$: the quadratic term has a nonzero expectation even when the random effect is centered. A related mistake is to claim that FO must always underestimate random-effect variance because it drops variability. In a fit, the variance parameter may increase to absorb omitted variability, as in the moment-matching example, or the fit may compensate in other ways. The correct conclusion is that FO can distort estimates, with the direction determined by the full model and estimation problem.

The Taylor calculations are local, and their usefulness depends on how well the expansion represents predictions over the relevant range of random effects. They do not quantify every source of error in an approximate likelihood. In particular, curvature in the mean prediction is only one mechanism; changes in residual variance or other features of the conditional distribution can also matter. The practical question is whether the approximation is adequate for the intended model, design, and inferential target.

Key takeaway: With nonlinear predictions, centered random effects can still shift the population mean and add variability. The leading mean correction is proportional to curvature times $\Omega$, while the scalar quadratic variance contribution is proportional to curvature squared times $\omega^4$. FO bias can therefore grow as random-effect variability or curvature increases, but its direction for any fitted parameter is model-dependent.

Check Your Understanding

Answer each question, then compare with the concise solutions.

For $\eta\sim N(0,\omega^2)$ and $f(\eta)=a+g\eta+\tfrac{1}{2}h\eta^2$, what is the second-order mean correction relative to the FO mean?
If $\omega^2$ doubles while $h$ stays fixed, by what factor does the quadratic variance contribution $\tfrac{1}{2}h^2\omega^4$ change?
Can a centered random effect shift the average prediction? State the condition in the quadratic example.
In a variance-matching FO fit, why might the fitted random-effect variance be larger than the generating variance?
Does omitted curvature imply that every FO estimate is biased upward? Explain briefly.

Solutions. (1) The correction is $\tfrac{1}{2}h\omega^2$. (2) It increases by a factor of four because the contribution is proportional to $(\omega^2)^2$. (3) Yes; the shift is $\tfrac{1}{2}h\omega^2$, which is nonzero when $h\ne0$ and $\omega^2>0$. (4) The linearized model may use a larger variance to represent variability contributed by the omitted quadratic term. (5) No. The direction depends on the model, data, and available parameter compensations.