Why compare approximations on the same concentration model?
FO, FOCE, and Laplace fitting address the same basic task: estimating population parameters when each subject has unobserved random effects. Their estimates can differ because they approximate the integration over those effects differently. In a nonlinear concentration model, the curvature of the individual prediction as the random effect changes can make that distinction consequential. A comparison is informative only when the model, data, estimation target, and objective convention are held fixed.
This tutorial uses a one-compartment intravenous bolus model with a random effect on the elimination rate. It focuses on carrying out and interpreting a method-to-method comparison, not on treating any one approximation as universally best.
Model and assumptions for the comparison
The model assumes a single exponential decline after the bolus, the same fixed initial concentration for every subject, normally distributed random effects on the log-rate scale, and independent normal residual errors with constant variance. For the simple comparisons below, subjects are independent and measurements are conditionally independent given the subject's random effect. In an application, these assumptions need to be justified by the design and checked against the observed data; the three approximation methods do not remove the assumptions.
The marginal likelihood for subject \(i\) integrates the conditional data density over the distribution of \(\eta_i\). FO approximates the prediction by a first-order Taylor expansion around \(\eta_i=0\). FOCE instead uses a local expansion around a subject-specific conditional mode; this comparison uses the interaction form, FOCE-I, which retains the dependence of the conditional residual variation on the random effect. Laplace approximates the integral locally around that mode using its curvature. Implementations can differ in details, so the specific approximation and objective convention should be documented when reporting results.
What is being approximated?
Write the conditional residual sum of squares as $$ RSS_i(\eta_i)=\sum_{j=1}^{n_i}\{y_{ij}-f_{ij}(\eta_i)\}^2. $$ Ignoring constants that do not depend on \(\eta_i\), the negative twice log of the joint contribution from subject \(i\)'s data and random effect is $$ q_i(\eta_i)=\frac{RSS_i(\eta_i)}{\sigma^2}+\frac{\eta_i^2}{\omega^2}. $$ FO replaces each \(f_{ij}(\eta_i)\) with \(f_{ij}(0)+f'_{ij}(0)\eta_i\). Under this linearized model, the subject's vector of observations has a normal marginal distribution with mean vector \(f_i(0)\) and covariance \(\sigma^2 I+\omega^2 g_i g_i^\mathsf{T}\), where \(g_i\) contains the derivatives \(f'_{ij}(0)\). The FO objective is formed by evaluating this approximate marginal model across subjects.
For FOCE-I and Laplace, define \(h_i(\eta_i)=q_i(\eta_i)/2\), with any parameter-independent density constants handled consistently. Let \(\hat\eta_i\) be a conditional mode, a value minimizing \(h_i\). Laplace uses the local curvature \(H_i=h_i''(\hat\eta_i)\). The conditional mode and curvature change as the population parameters change, so they are part of the numerical fitting problem, not fixed data summaries.
A local justification for the Laplace contribution
At a well-behaved interior minimum, \(h_i'(\hat\eta_i)=0\) and \(H_i>0\). A second-order Taylor expansion near the mode gives
Substituting this quadratic form into the integral gives a Gaussian integral. Specifically, \(\int \exp\{-H_i(\eta_i-\hat\eta_i)^2/2\}\,d\eta_i=\sqrt{2\pi/H_i}\). Thus, up to the stated density constants, the subject's approximate negative twice log marginal contribution is
The curvature term is essential: the Laplace contribution is not just the conditional objective evaluated at the mode. The approximation relies on the log integrand being well represented by a quadratic over the region that contributes appreciably to the integral. Strong asymmetry, a nearly flat mode, or multiple relevant modes can undermine that representation.
Worked Example: Comparing Three Fits to Hypothetical Concentration Data
Suppose four hypothetical subjects are measured at 1, 2, and 4 hours after a bolus. The recorded concentrations, in mg/L, are as follows. The initial concentration is fixed at 5 mg/L as specified in the model.
| Subject | 1 h | 2 h | 4 h |
|---|---|---|---|
| A | 4.200 | 3.565 | 2.533 |
| B | 4.050 | 3.310 | 2.135 |
| C | 4.340 | 3.770 | 2.850 |
| D | 4.160 | 3.420 | 2.300 |
For illustration, assume numerical minimization of each method's specified approximate objective converged for these data. The following rounded results are hypothetical optimizer output, included to show how a same-data comparison is read; they are not a general prediction of how the methods must differ.
| Method | \(k\) (h\(^{-1}\)) | \(\omega^2\) (log-rate scale) | \(\sigma\) (mg/L) | Approximate objective |
|---|---|---|---|---|
| FO | 0.206 | 0.025 | 0.105 | 12.84 |
| FOCE-I | 0.201 | 0.030 | 0.091 | 11.27 |
| Laplace | 0.199 | 0.034 | 0.087 | 10.96 |
The FO estimate corresponds to a typical rate of \(0.206\ h^{-1}\), while the Laplace fit gives \(0.199\ h^{-1}\). The FOCE-I and Laplace random-effect variances are larger than the FO estimate in this example. Because \(\omega^2\) is a variance on the log-rate scale, it is dimensionless; it is not a variance in concentration units. The residual standard deviation decreases from 0.105 to 0.087 mg/L across the displayed fits.
The approximate objective is lower for the displayed FOCE-I fit than for FO, and lower again for Laplace. The differences are \(12.84-11.27=1.57\) and \(11.27-10.96=0.31\). These are differences between method-specific approximations evaluated at their fitted parameter values. They are not exact likelihood-ratio statistics and, on their own, do not establish that one method has recovered the true parameters or provides a universally preferable fit.
Worked Example: Calculating a Conditional Objective for One Subject
Use subject A's hypothetical observations at 1, 2, and 4 hours and consider candidate values \(k=0.200\ h^{-1}\), \(\eta_A=-0.137\), \(\sigma=0.050\) mg/L, and \(\omega^2=0.040\). These are candidate values for this calculation, not a claim that they are the fitted values from the preceding table.
At 1 hour, the predicted concentration is \(5\exp\{-0.200\exp(-0.137)\}\) mg/L, so the residual based on the unrounded observation is $$ 4.20-5\exp\{-0.200\exp(-0.137)\}\approx 0.000171\ \text{mg/L}. $$ At 2 and 4 hours, the residuals are approximately 0.037 and 0.044 mg/L, respectively. The small first residual rounds to 0.000 mg/L at three decimal places; retaining its unrounded value does not materially affect the following sum.
Using the defined conditional objective, $$ q_A=\frac{0.000171^2+0.037^2+0.044^2}{0.050^2}+\frac{(-0.137)^2}{0.040}. $$ The residual contribution is approximately \((0.001369+0.001936)/0.0025=1.322\). The random-effect penalty is \(0.018769/0.040=0.4692\). Therefore, \(q_A\approx1.322+0.4692=1.791\). This is a conditional joint objective at the stated candidate effect; it is not by itself the marginal likelihood contribution, and the candidate \(\eta_A\) has not been shown to be the conditional mode.
Worked Example: Seeing the Curvature Term in a Laplace Calculation
For a one-subject illustrative calculation, suppose a conditional mode has \(h(\hat\eta)=0.500\), the curvature there is \(H=6.00\) per squared random-effect unit, and \(\omega^2=0.040\). The curvature-and-random-effect part of the Laplace contribution, before adding the residual-density normalization for this one-observation illustration, is $$ 2h(\hat\eta)+\log(\omega^2H)=2(0.500)+\log(0.040\times6.00). $$ Since \(0.040\times6.00=0.240\) and \(\log(0.240)\approx-1.427\), this contribution is approximately \(1.000-1.427=-0.427\). This is only the curvature-and-random-effect part: because \(n=1\), the full contribution additionally includes \(\log(2\pi\sigma^2)\), and cannot be evaluated numerically here because \(\sigma\) is unspecified. The calculation shows why the mode value alone, \(2h(\hat\eta)=1.000\), does not give the Laplace approximation: the local curvature changes the approximate integral.
How to read a method comparison
Use the same subjects, observations, structural model, residual model, parameter constraints, and treatment of constants for each fit.
Confirm that each method has reached a stable solution and that the reported estimates are not merely intermediate values from different optimization states.
Describe changes in \(k\), \(\omega^2\), and \(\sigma\), then report method-specific objective values and their differences with the objective convention stated.
A lower approximate objective is meaningful as a comparison only with care. FO, FOCE-I, and Laplace use different local representations of the same marginalization problem. In particular, a difference between approximations is not automatically distributed as a standard likelihood-ratio statistic. Objective differences are most useful here as a signal to inspect the approximation, alongside parameter shifts and diagnostic evidence.
Diagnostics, alternatives, and limitations
The earlier material on diagnostic patterns of approximation bias emphasizes that estimates alone do not reveal whether an approximation is behaving well. For this concentration model, inspect conditional residuals against time and predicted concentration, subject-level fitted curves, estimated random effects, and the stability of the conditional modes and curvature. A systematic time trend in residuals may indicate a structural or residual-model problem rather than an integration-approximation problem. Large changes in \(k\) or \(\omega^2\) across methods, unstable modes, or very small curvature \(H_i\) warrant particular attention.
The main alternatives differ in where they represent the likelihood:
| Method | Local representation | What to watch |
|---|---|---|
| FO | Linearization of predictions around zero random effects | May be inaccurate when predictions are strongly nonlinear over the random-effect distribution. |
| FOCE-I | Conditional expansion around subject-specific modes, retaining interaction | Depends on locating stable conditional modes; mode-based local behavior may not represent distant regions well. |
| Laplace | Quadratic approximation to the log integrand around its mode | Can be unreliable with strong asymmetry, multiple important modes, or poor local quadratic fit. |
None of these checks proves that the structural model is correct. Nor does a close match among the three methods guarantee accuracy: all can share misspecification, such as an unsuitable one-compartment structure, an incorrect residual-error assumption, or an inadequate random-effect distribution. The approximations also do not substitute for assessing identifiability and whether the sampling times contain information about the parameters.
Common mistake
Check Your Understanding
Try each question first, then compare with the solutions below.
- In the specified model, what does \(\omega^2\) describe, and what are its units?
- If \(H_i\) is close to zero and positive at a conditional mode, why should the Laplace result be examined carefully?
- For approximate objectives 9.80 under FO and 9.10 under FOCE-I, calculate the difference and state one conclusion that is not warranted.
- A residual pattern shows concentrations increasingly underpredicted at late times for most subjects under every approximation. Is changing from FO to Laplace necessarily the solution? Explain.
Solutions. (1) \(\omega^2\) is the variance of the subject-specific log-rate effect \(\eta_i\); it is dimensionless because the log-rate deviation is dimensionless. (2) The Laplace contribution includes \(\log(\omega^2 H_i)\). Near-zero curvature indicates a broad, relatively flat region, so a quadratic local approximation may not adequately represent the integral. (3) The difference is \(9.80-9.10=0.70\), with FOCE-I lower by 0.70. That does not prove FOCE-I is more accurate or that the difference is a likelihood-ratio statistic. (4) No. A pattern shared by all methods points first to checking the structural or residual model and the data; changing the integration approximation alone need not fix systematic model misspecification.