Why Approximation Choice Depends on the Likelihood
In a nonlinear mixed-effects model, each individual contributes a likelihood averaged over a distribution of random effects. That average is generally an integral with no convenient closed form. FO, FOCE, and Laplace approximate it in different ways, so the choice should reflect the shape of the individual likelihood—not simply a preference for one method or the speed of a software default.
The key questions are whether the predictions change nearly linearly over plausible random effects, whether the residual-error distribution changes with those effects, and whether the conditional likelihood has a single, well-behaved mode. These considerations inform a starting choice. They do not replace fitting diagnostics or sensitivity checks.
FO uses a first-order Taylor expansion of the model prediction around the population value \(\eta_i=0\). It is most credible when prediction curvature over the relevant random-effect range is small. FOCE expands around an individual conditional estimate, commonly the conditional mode, so it can better represent individuals whose data support effects away from zero. Laplace approximates the log of the integrand quadratically around its conditional mode, retaining local curvature in the likelihood integral.
These are not interchangeable labels for the same calculation. In particular, FOCE’s expansion point can help with displaced individual effects, while Laplace uses the local curvature to approximate the integral itself. The earlier material on nonlinearity and Laplace curvature establishes why curvature matters; here it is used as part of a selection framework.
A Compact Justification for the Laplace Criterion
Write the log integrand for an individual as \(h_i(\eta)=\log f(y_i\mid\eta;\theta)+\log\phi(\eta;0,\Omega)\). Suppose it has a unique interior maximum at \(\hat\eta_i\), and define \(H_i=-\nabla^2h_i(\hat\eta_i)\). If \(H_i\) is positive definite and the integrand is concentrated near that mode, a second-order Taylor expansion gives \(h_i(\eta)\approx h_i(\hat\eta_i)-\tfrac12(\eta-\hat\eta_i)^\mathsf{T}H_i(\eta-\hat\eta_i)\). Integrating this local Gaussian approximation yields the Laplace log-likelihood approximation.
where \(q\) is the dimension of \(\eta_i\). The determinant term reflects the volume around the mode: a sharply concentrated integrand has larger curvature and smaller local volume. This derivation also identifies important conditions. A mode must be found, its negative Hessian must be positive definite, and the local quadratic shape must represent the important mass of the integrand. A numerical optimizer returning a mode does not, by itself, establish these conditions.
Worked Examples: Using Curvature and Error Structure
Worked Example: When FO Is a Reasonable Starting Point
Consider a hypothetical concentration model with prediction \(m(\eta)=100e^\eta\) mg/L, random-effect standard deviation \(0.05\), and additive residual standard deviation \(10\) mg/L. Assess the linearization at zero over a range of two random-effect standard deviations, \(-0.10\leq\eta\leq0.10\).
The FO approximation replaces \(e^\eta\) with \(1+\eta\). At the upper endpoint, the omitted second-order remainder in concentration units is $$ 100\{e^{0.10}-(1+0.10)\} =100(1.10517-1.10000) =0.517\text{ mg/L}. $$ At the lower endpoint it is $$ 100\{e^{-0.10}-(1-0.10)\} =100(0.90484-0.90000) =0.484\text{ mg/L}. $$ Both are small relative to the additive residual standard deviation of \(10\) mg/L. The linearized prediction changes by approximately \(10\) mg/L at either endpoint, while the curvature remainder is about half a milligram per litre. For this hypothetical scale of random effects and error, FO is a defensible starting approximation.
This is not a proof that FO will be adequate in a whole data set. The calculation checks prediction curvature over a chosen range; it does not check whether the residual distribution is correctly represented, whether other individuals have more informative data, or whether parameter estimates are sensitive to the approximation. Those checks remain necessary.
Worked Example: A Displaced Mode and Proportional Error
For one hypothetical observation, let \(y=130\) mg/L, \(m(\eta)=100e^\eta\) mg/L, the proportional residual standard deviation be \(0.20m(\eta)\), and \(\eta\sim N(0,0.30^2)\). Ignoring constants independent of \(\eta\), the log integrand is $$ h(\eta)=-\frac{1}{2}\left(\frac{y-m(\eta)}{0.20m(\eta)}\right)^2-\log\{0.20m(\eta)\}-\frac{\eta^2}{2(0.30)^2}. $$
Let \(r=y/m(\eta)\). Differentiating gives $$ h'(\eta)=\frac{r(r-1)}{0.04}-1-\frac{\eta}{0.09}. $$ Solving \(h'(\eta)=0\) numerically gives a conditional mode near \(\hat\eta=0.1645\). At that point, \(m=100e^{0.1645}\approx117.9\) mg/L, \(r\approx1.103\), and substitution into the derivative gives approximately \(2.84-1-1.83=0.01\), close to zero given rounding. The mode is displaced by \(0.1645/0.30\approx0.548\) random-effect standard deviations from zero. That displacement makes an expansion only at zero less attractive than one centered on the individual data.
The curvature at the mode can also be checked. Differentiating again gives $$ h''(\eta)=-\frac{r(2r-1)}{0.04}-\frac{1}{0.09}. $$ At \(r\approx1.103\), the negative second derivative is $$ H=-h''(\hat\eta)\approx\frac{1.103(1.206)}{0.04}+11.11\approx44.4. $$ This is positive, as required for a local maximum. FOCE is a more plausible starting choice than FO because the conditional mode is away from zero. Laplace is also a candidate if the profile around the mode is approximately symmetric and no competing mode is present.
The proportional-error term matters: its standard deviation is \(0.20m(\eta)\), so both the standardized residual and the log-scale contribution \(-\log\{0.20m(\eta)\}\) vary with \(\eta\). A calculation that treats the residual standard deviation as fixed would describe a different likelihood. The mode and curvature checks support considering Laplace, but one local Hessian alone cannot establish that its quadratic approximation captures the whole integral.
Worked Example: Choosing Laplace When the Local Shape Is Stable
Suppose a hypothetical rich-sampling data set is fitted with two standardized random effects \(z\), where each component is measured in its own random-effect standard-deviation units. The optimizer and profile checks report one conditional mode at \(\hat z=(0.20,-0.20)\). In the principal curvature directions, the negative Hessian has eigenvalues 25 and 100; the checked profiles are nearly symmetric near the mode and show no competing maximum. These are stipulated hypothetical diagnostics from the fit.
The displacement from zero is modest: each coordinate is \(0.20\) standard deviations in absolute value. More importantly for Laplace, both curvature eigenvalues are positive, so the local Hessian is positive definite. Its determinant is their product, \(25(100)=2500\). For \(q=2\), the curvature contribution to the Laplace approximation is $$ \frac{2}{2}\log(2\pi)-\frac12\log(2500) =1.8379-\frac12(7.8240) =-2.074 $$ to three decimals. Thus the approximate log marginal likelihood is \(h(\hat z)-2.074\), with \(h\) evaluated using a consistently normalized density.
The positive curvatures, a single mode, and nearly symmetric local profiles make Laplace defensible for this hypothetical fit. Rich within-subject information can help define a concentrated individual likelihood, but it does not guarantee a Gaussian shape. The profile check is therefore part of the justification, not an optional cosmetic step.
A Practical Decision Framework
A useful decision is conditional: choose a candidate approximation from the model and diagnostics, fit it, then test whether the evidence supports it. The following sequence keeps the relevant considerations distinct.
Compare the nonlinear prediction with its first-order expansion over a realistic range of random effects. Small remainders support considering FO; substantial curvature argues against relying on an expansion at zero without further checks.
Large random-effect variability can expose more of the nonlinear prediction curve. Sparse observations may leave individual effects weakly identified; rich data may sharpen the conditional likelihood but can also reveal substantial nonlinearity. Neither sampling pattern alone selects a method.
With additive error, the residual standard deviation need not change with the prediction. With proportional error, it typically does. Check whether changes in prediction alter the residual scale or standardized residual across the relevant random-effect range; a fixed-scale linearization may then be inadequate.
Record the mode’s displacement from zero, the eigenvalues of its negative Hessian, and whether profiles are asymmetric or multimodal. A displaced mode can favor FOCE over FO; a stable, approximately quadratic, single-mode integrand supports considering Laplace.
Check convergence, parameter estimates, standard errors, and variance components across plausible approximations. Investigate materially different results rather than selecting whichever fit has the preferred objective value.
Comparison, Diagnostics, and Limitations
| Approximation | Useful when | Warning signs |
|---|---|---|
| FO | Prediction changes nearly linearly around zero across plausible effects, with limited residual-scale interaction. | Large curvature, appreciable conditional-mode displacement, or strong dependence of residual scale on effects. |
| FOCE | Individual modes are meaningfully displaced from zero and a local expansion around those modes is more representative. | Unstable mode estimates, marked asymmetry, or multiple modes that make the chosen expansion point consequential. |
| Laplace | A dominant, interior mode has positive-definite curvature and the integrand is close to quadratic in the region that contributes materially. | Near-zero or negative curvature eigenvalues, strong skewness, heavy tails, or competing modes. |
Diagnostics should test the assumptions behind the choice. For Laplace, verify that the optimizer converges from multiple starting values, that the negative Hessian is positive definite, and that profiles around the mode do not reveal important skewness or another maximum. For FO or FOCE, compare the linearized predictions with the nonlinear predictions over plausible random effects and inspect whether residual-error scales are being represented consistently. In every case, nonconvergence, sensitivity to starting values, unstable variance estimates, or large shifts in fixed effects across methods warrant investigation.
No approximation is universally best. FO can be efficient when its local linearization is adequate, but may misrepresent a strongly curved likelihood. FOCE can accommodate a displaced expansion point, but still relies on a local approximation. Laplace uses curvature information, but may be poor when the integrand is asymmetric, diffuse, or multimodal. If diagnostics are unfavorable or estimates depend materially on the method, a more accurate integration-based check can help assess the approximation; it does not make a poor model or an unidentified parameter well behaved.
Common Pitfall
Wrong reasoning: “The random-effect standard deviation is small, so FO must be correct,” or “the Hessian is positive, so Laplace must be accurate.” Neither conclusion follows. A small standard deviation does not rule out an important residual-scale interaction or influential tail behavior. A positive Hessian verifies local maximum curvature, not that the integrand is globally dominated by a symmetric quadratic neighborhood. Correction: treat each diagnostic as evidence about a specific condition, then check the other conditions and assess sensitivity across plausible approximations.
Likewise, objective values from different approximations are not a stand-alone contest for choosing the truth. Each method changes the likelihood calculation. Interpret differences alongside parameter stability, convergence, curvature, and model adequacy.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- Why does a conditional mode far from zero make FO less attractive, and which approximation may better accommodate that displacement?
- A Laplace fit has a positive-definite Hessian but a strongly asymmetric profile. What has been checked, and what remains uncertain?
- Why can proportional residual error complicate an expansion around zero?
- A model gives similar estimates under FO and FOCE but a markedly different variance estimate under Laplace. What is a sound next step?
Solutions. (1) FO expands at zero, so it may approximate the model poorly near the data-supported effect; FOCE expands around the conditional mode. (2) Positive definiteness supports a local maximum, but asymmetry questions whether the quadratic approximation represents enough of the integrand; inspect the profile and consider sensitivity to another likelihood evaluation. (3) The residual standard deviation depends on the prediction, which changes with the random effect, so the residual scale and standardized residual can also change. (4) Check convergence, mode and curvature diagnostics, and the stability of estimates; investigate why the variance component changes rather than choosing a method solely by its objective value.