Why the Estimation Method Matters
In a nonlinear mixed-effects pharmacokinetic model, each subject’s likelihood depends on unobserved random effects. Estimation methods differ in how they approximate the integral over those effects. Choosing among first-order estimation (FO), first-order conditional estimation (FOCE), and the Laplace approximation is therefore not just a software setting: the approximation can affect parameter estimates, uncertainty, and the apparent fit.
The earlier material on first-order linearization, FOCE conditional likelihood, and the Laplace approximation introduced each method. Here the focus is comparison: what each approximation retains, which models are favorable or unfavorable to it, and how to check whether the choice is consequential. These methods are not interchangeable guarantees of accuracy. Their performance depends on the likelihood’s local shape, the amount of nonlinearity, the random-effects distribution, and the observation design.
Shared Model and Approximation Targets
For subject \(i\), let \(y_i\) be the vector of observations, \(\eta_i\) a \(q\)-dimensional vector of random effects, and \(f_i(\theta,\eta_i)\) the model-predicted observations. A standard starting model is \(y_i=f_i(\theta,\eta_i)+\epsilon_i\), with \(\eta_i\sim N(0,\Omega)\). In the simplest case, \(\epsilon_i\sim N(0,R_i)\), independently of \(\eta_i\), and subjects are independent conditional on the population parameters. The subject likelihood is \(p(y_i\mid\theta,\Omega)=\int p(y_i\mid\eta_i,\theta) p(\eta_i\mid\Omega)\,d\eta_i\).
The integral is generally not available in closed form when the predictions or residual variance are nonlinear functions of the random effects. FO, FOCE, and Laplace approximate this integral in different ways. Their usual Gaussian versions presume that the random-effects distribution and residual model are correctly specified, and that the model is sufficiently smooth for local approximations to be meaningful. If the residual variance depends on \(\eta_i\), or residual errors are non-Gaussian, the corresponding likelihood and approximation must account for that structure; the fixed-\(R_i\) formulas below do not automatically apply.
What the Methods Approximate
Write \(J_i\) for the derivative matrix of \(f_i\) with respect to \(\eta_i\). FO replaces the nonlinear prediction by \(f_i(\theta,\eta_i)\approx f_i(\theta,0)+J_i\eta_i\). Under the independent Gaussian assumptions above, this gives the approximate marginal distribution
This approximation is exact when the prediction is affine in \(\eta_i\) and the residual covariance does not depend on \(\eta_i\). Otherwise, it treats the prediction’s local slope at zero as representative of the random-effects range. In particular, it does not shift the expansion point to accommodate a subject whose data strongly favor a nonzero random effect.
FOCE instead finds a conditional estimate for each subject, typically the mode of the posterior density of \(\eta_i\) given \(y_i\), and expands the prediction locally there. This can represent subject-specific deviations and local slopes better than an expansion at zero. In a local-linear implementation, the approximation integrates a Gaussian model formed from that conditional expansion. A common conceptual distinction from Laplace is that FOCE’s local approximation is based on the linearized observation model, whereas Laplace uses the curvature of the full log integrand. Software implementations and options differ, so the actual objective function should be checked rather than inferred from a method label alone.
For Laplace, define \(g_i(\eta)=\log p(y_i\mid\eta,\theta)+\log p(\eta\mid\Omega)\), including the density normalizing constants, and let \(\hat\eta_i\) maximize \(g_i\). If \(A_i=-\nabla^2g_i(\hat\eta_i)\) is positive definite, the approximation is
The determinant term represents the width of the integrand around its mode. Laplace is exact when the log integrand is quadratic, as in the affine Gaussian case. With a non-quadratic integrand it remains an approximation: a good mode alone is not enough if the local curvature is unrepresentative, the distribution is strongly skewed, or important mass lies away from that mode.
Worked Example: An Affine Model Where FO Is Exact
Worked Example: An Affine Model Where FO Is Exact
Suppose one hypothetical concentration observation follows \(y=10+2\eta+\epsilon\), where \(\eta\sim N(0,0.25)\), \(\epsilon\sim N(0,1)\), and the errors are independent. The prediction is affine in \(\eta\), the residual variance is constant, and both distributions are Gaussian, so the marginal calculation is exact. For an observed concentration of 11 units, the marginal mean is 10 units and the variance is \(2^2(0.25)+1=2\) concentration-units squared.
Thus the log likelihood contribution is
The conditional mean of \(\eta\) given \(y=11\) is \(\operatorname{Cov}(\eta,y)\operatorname{Var}(y)^{-1}(11-10)\). Here \(\operatorname{Cov}(\eta,y)=2(0.25)=0.5\), so the conditional estimate is \(0.5/2=0.25\). In this affine Gaussian case, FO’s population-centered approximation is exact, and an exact Gaussian integration, including a Laplace calculation, gives the same marginal likelihood. This is a useful benchmark, not evidence that the methods agree in nonlinear models.
Worked Example: Different Local Approximations in a Nonlinear Model
Worked Example: Different Local Approximations in a Nonlinear Model
Consider a hypothetical scalar observation \(y=11\) with prediction \(f(\eta)=10+\eta+\eta^2\), residual error \(\epsilon\sim N(0,1)\), and \(\eta\sim N(0,1)\), independently. All quantities are on a consistent measurement scale. The random-effects variance is one, the residual variance is fixed, and the smooth prediction permits local expansions. The quadratic term makes the prediction nonlinear.
At \(\eta=0\), FO has prediction 10 and slope \(f'(0)=1\). Its approximate marginal variance is \(1^2(1)+1=2\), giving
For the conditional expansion, the log integrand excluding normalizing constants is \(h(\eta)=-\tfrac12(1-\eta-\eta^2)^2-\tfrac12\eta^2\). Its derivative is \(h'(\eta)=1-3\eta^2-2\eta^3\), which is zero at \(\eta=0.5\). The other stationary point, \(\eta=-1\), is not a maximum; evaluating \(h\) gives \(h(0.5)=-0.15625\) and \(h(-1)=-1\). The tails tend to minus infinity, so the mode is \(\hat\eta=0.5\).
At this mode, \(f(\hat\eta)=10.75\), \(f'(\hat\eta)=2\), and the residual is \(11-10.75=0.25\). A local-linear conditional approximation has curvature \(1+(2)^2=5\): the first term is the residual precision times the squared slope, and the second is the prior precision. Its approximate log likelihood is
Laplace uses the full curvature. In this example, \(h''(\eta)=-(f'(\eta))^2+(y-f(\eta))f''(\eta)-1\), so at \(\eta=0.5\) the negative second derivative is \(A=4-0.25(2)+1=4.5\). Therefore
The difference between the local-linear and Laplace values arises because the latter retains the residual-times-second-derivative contribution to curvature. These are approximations, not known exact likelihood values; their disagreement signals that the approximation choice matters for this example. The calculation does not establish which estimate is closest to the exact integral.
Worked Example: High Variability and a Poor FO Expansion
Worked Example: High Variability and a Poor FO Expansion
Suppose an individual elimination-rate parameter is \(k_i=k_{\mathrm{pop}}\exp(\eta_i)\), with \(k_{\mathrm{pop}}=1\ \mathrm{h}^{-1}\) and \(\eta_i\sim N(0,0.64)\). This is a hypothetical lognormal random-effects model with a smooth prediction parameterization and substantial variability. The population parameter \(k_{\mathrm{pop}}\) is the median of \(k_i\), not its arithmetic mean.
For a normal random variable, \(E[\exp(\eta_i)]=\exp(0.64/2)\), so \(E[k_i]=1\cdot\exp(0.32)\approx1.377\ \mathrm{h}^{-1}\). Also, \(\operatorname{Var}(\exp(\eta_i))=\exp(0.64)(\exp(0.64)-1)\), giving
FO’s local expansion at zero is \(k_i\approx k_{\mathrm{pop}}(1+\eta_i)\). It therefore gives an approximate mean of \(1\ \mathrm{h}^{-1}\) and variance \(0.64\ (\mathrm{h}^{-1})^2\), compared with the exact lognormal moments above. The discrepancy follows from curvature in \(\exp(\eta_i)\), which the linear approximation discards. This example does not imply that FO always fails at this variability level: the impact on a data likelihood also depends on how the parameter enters the observation model and on the sampling design.
Choosing and Checking an Approximation
The following comparison is a starting point, not a universal ranking. Computational cost depends on the implementation, number of random effects, number of observations per subject, and numerical settings.
| Feature | FO | FOCE | Laplace |
|---|---|---|---|
| Expansion point | Population mean, \(\eta_i=0\) | Subject-specific conditional estimate | Mode of the full log integrand |
| Local information retained | Prediction slope at zero | Prediction slope at the conditional estimate | Mode and Hessian of the full log integrand |
| Favorable setting | Approximately linear predictions and modest variability | Nonlinear predictions where subject-specific modes improve the local expansion | A smooth, concentrated integrand well represented by a quadratic near its mode |
| Common concern | Bias when local behavior at zero is unrepresentative | Local linearization may still miss curvature or dependence in the residual model | Curvature or a single mode may not represent the full integral |
| Relative computational demand | Often lowest | Usually higher because conditional estimates are needed | Requires modes and curvature calculations; cost can rise with random-effect dimension |
For a nearly linear model with low or moderate variability and a well-behaved residual model, FO can be a reasonable and computationally efficient choice. When predictions are nonlinear over the plausible random-effects range, FOCE’s subject-specific expansion can improve the approximation. Laplace can be attractive when the conditional integrand is smooth, concentrated, and close to quadratic, particularly when its curvature is informative. High variability, sparse observations, strong residual-error dependence on random effects, or multiple modes warrant extra caution with all three.
Diagnostics should focus on whether the approximation is stable and consequential, not on a single convergence message. Useful checks include:
- Fit plausible alternatives from comparable starting values and compare parameter estimates, standard errors, objective values, and whether the same substantive conclusions hold.
- Inspect conditional estimates and gradients or curvature for extreme values, boundary behavior, or subjects whose likelihood is poorly represented by a local expansion.
- Check whether a Laplace Hessian is positive definite at the reported mode. A non-positive-definite Hessian prevents the stated local Gaussian approximation.
- Assess predictions and residual behavior across time and covariate ranges. A stable objective does not show that the structural or residual model is correctly specified.
- Check sensitivity to initial estimates and numerical settings. A method comparison is informative only if the fits have reached comparable solutions.
FOCE and Laplace may differ little when the local integrand is close to quadratic, but agreement is not proof of model correctness. Conversely, a difference does not automatically identify the more accurate method. When an approximation choice changes estimates or conclusions, report the method and investigate the model and design rather than selecting a result solely because it is numerically convenient. A local approximation can be unreliable when the conditional distribution is broad, skewed, multimodal, or constrained by parameterization. Method choice cannot repair misspecified random-effects distributions, residual errors, or structural predictions.
Common Pitfall
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- Under what conditions is the FO marginal Gaussian calculation exact?
- Why can a subject-specific expansion point help when FO is inaccurate?
- In the nonlinear worked example, why is the Laplace curvature 4.5 rather than 5?
- For the lognormal rate example, what does \(k_{\mathrm{pop}}\) represent, and why does the FO mean differ from the exact arithmetic mean?
- Name two checks that are useful when alternative approximations produce different estimates.
Solutions. (1) The prediction must be affine in the random effects, with Gaussian random effects, independent Gaussian residuals, and residual covariance not depending on those effects. (2) A conditional mode can reflect the subject’s data-supported deviation and local prediction slope, rather than assuming the slope at zero represents that subject. (3) Laplace includes the second derivative of the nonlinear prediction multiplied by the residual; here \(A=1+(2)^2-0.25(2)=4.5\). The linearized conditional approximation omits that curvature contribution and uses 5. (4) \(k_{\mathrm{pop}}\) is the median because \(\eta_i\) has median zero; the FO expansion discards the exponential function’s curvature, so its mean approximation is \(k_{\mathrm{pop}}\), not \(k_{\mathrm{pop}}\exp(0.64/2)\). (5) Examples include comparing estimates and standard errors from fits using comparable starting values, checking the Laplace Hessian for positive definiteness, and examining conditional estimates and residual behavior.