A Simulation Study of FO, FOCE, and Laplace

Why a Simulation Comparison Needs a Design

FO, FOCE, and Laplace are alternative approximations to the marginal likelihood in a nonlinear mixed-effects model. A simulation study can show how their behavior changes under controlled conditions: the number of subjects, the amount of between-subject variability, and the curvature of the individual prediction function. The goal is not to declare one approximation universally best. It is to determine which differences arise under the specified model and design, and how large those differences are relative to simulation uncertainty.

The comparison is most informative when each approximation is applied to the same simulated datasets. This paired design holds the generated data fixed while changing the estimation method. It also makes method differences within a replicate interpretable, rather than mixing approximation effects with differences between independently generated samples.

Definition: A simulation replicate is one dataset generated from a specified model and parameter vector, followed by the planned fitting procedures. For a parameter $\theta_k$, its Monte Carlo bias over $R$ replicates for method $m$ is the average fitted-minus-true difference: $$ \widehat{\operatorname{Bias}}_{m,k} = \frac{1}{R}\sum_{r=1}^{R} \left(\hat{\theta}_{m,k,r}-\theta_{0,k}\right). $$ Here $\theta_{0,k}$ is the generating value and $\hat{\theta}_{m,k,r}$ is the estimate from replicate $r$. Positive bias indicates average overestimation on the parameter's original scale.

Specify the Model and Vary One Factor at a Time

A useful design begins with a fully specified generating model. For illustration, suppose subject $i$ has two observations at scaled design points $x_j=1$ and $x_j=2$:

y_{ij}=\beta\exp(\kappa x_j\eta_i)+\epsilon_{ij}, \qquad \eta_i\sim N(0,\omega^2), \qquad \epsilon_{ij}\sim N(0,\sigma^2).

The random effects $\eta_i$ are independent across subjects, the residual errors are independent across observations and subjects, and the two error sources are independent. The parameter $\beta$ sets the response scale, $\omega^2$ is the random-effect variance, and $\sigma^2$ is the residual variance. The dimensionless multiplier $\kappa$ changes how strongly the prediction curves with $\eta_i$. This is a hypothetical model for designing a comparison, not a description of a real study.

For example, fix $\beta=1.20$ response units and $\sigma=0.15$ response units, and cross three factors: $N$, the number of subjects; $\omega^2$, the random-effect variance; and $\kappa$, the curvature multiplier. A $2\times2\times2$ factorial design with $N\in\{50,200\}$, $\omega^2\in\{0.04,0.36\}$, and $\kappa\in\{0.5,1.5\}$ has eight conditions. Keeping the residual model and observation schedule fixed helps attribute changes to the factors being studied. If the scientific question concerns sparse versus rich sampling, that should be introduced as an additional planned factor rather than changed informally between conditions.

FO linearizes the individual prediction around the population value of the random effect. FOCE uses conditional information about each subject's random effect in its approximation, while Laplace uses the conditional mode and local curvature of the individual contribution. Their approximations can therefore respond differently to nonlinear prediction and random-effect variability. The earlier material on the Laplace approximation and FOCE objective differences establishes these method definitions; here the concern is how to compare their behavior in repeated simulated data.

Worked Example: Laying Out a Factorial Simulation

Suppose the planned factors are $N=50$ or $200$, $\omega^2=0.04$ or $0.36$, and $\kappa=0.5$ or $1.5$, using the model above. Choose $R=500$ replicates per condition.

There are $2\times2\times2=8$ conditions, so the study generates $8\times500=4{,}000$ datasets. Fit FO, FOCE, and Laplace to every dataset: this produces $4{,}000\times3=12{,}000$ planned fits before accounting for any retries or reference calculations.

Each replicate uses the same generated dataset for all three methods. Thus, for example, the FO and Laplace estimates in replicate 27 under $N=50,\omega^2=0.36,\kappa=1.5$ can be compared as a paired difference. The high-variance, high-curvature condition is not presumed to have a particular outcome; it is included to test whether the methods behave differently under that combination.

Summarize Parameter Performance

Bias alone is not enough. A method can have small average bias but substantial replicate-to-replicate variability. For each parameter and condition, report at least the average estimate, bias, and root mean squared error (RMSE). For positive parameters such as $\omega^2$, relative bias can supplement bias when comparisons across different generating scales are useful.

\widehat{\operatorname{RMSE}}_{m,k} = \sqrt{ \frac{1}{R}\sum_{r=1}^{R} (\hat{\theta}_{m,k,r}-\theta_{0,k})^2 }, \qquad \widehat{\operatorname{Relative\ Bias}}_{m,k} = \frac{\widehat{\operatorname{Bias}}_{m,k}}{\theta_{0,k}}.

Relative bias is undefined when the generating value is zero, so report the absolute bias in that case. For a variance parameter, state whether summaries are on the variance scale or the standard-deviation scale; those are different parameters and their biases need not match. Report estimates on the scale used in the model, with units where applicable.

Monte Carlo error describes the numerical uncertainty in a simulation summary due to using a finite number of replicates. If $d_r=\hat{\theta}_{m,k,r}-\theta_{0,k}$, let $s_d$ be the sample standard deviation of the $R$ values $d_r$. The estimated Monte Carlo standard error (MCSE) of the estimated bias is $s_d/\sqrt{R}$. It is not the standard error of a parameter estimate from one fitted dataset.

Worked Example: Calculating Bias and RMSE

In one hypothetical condition, let the true value be $\beta=1.20$ response units. Across five pilot replicates, suppose FO estimates are $1.08, 1.16, 1.21, 1.12,$ and $1.18$. Laplace estimates on the same datasets are $1.16, 1.19, 1.23, 1.17,$ and $1.20$. These five replicates illustrate the calculations; they are not enough to establish stable performance.

For FO, the mean estimate is $(1.08+1.16+1.21+1.12+1.18)/5=1.15$ response units, so its estimated bias is $1.15-1.20=-0.050$ response units. The squared errors are $0.0144, 0.0016, 0.0001, 0.0064,$ and $0.0004$. Therefore,

\widehat{\operatorname{RMSE}}_{\mathrm{FO}} = \sqrt{\frac{0.0144+0.0016+0.0001+0.0064+0.0004}{5}} = \sqrt{0.00458} \approx 0.0677

response units. For Laplace, the mean estimate is $1.19$ response units, giving bias $-0.010$ response units. Its squared errors are $0.0016, 0.0001, 0.0009, 0.0009,$ and $0$, so its RMSE is $\sqrt{0.0035/5}=\sqrt{0.0007}\approx0.0265$ response units. In this small illustrative set, Laplace has smaller absolute bias and RMSE; the calculation does not establish that it will do so in other conditions.

Compare Objective Functions on a Common Basis

Objective-function comparisons require more care than parameter summaries. Define the objective as $-2$ times the log likelihood, using the same likelihood constants and scaling for every method. For a given dataset, evaluate the methods at the same parameter vector before comparing their objective values. One possible reference is the marginal likelihood evaluated accurately enough that its numerical error is negligible for the intended comparison.

Let $\operatorname{OBJ}_{m,r}(\theta_0)$ be method $m$'s objective in replicate $r$, evaluated at the generating parameter vector, and let $\operatorname{OBJ}_{\mathrm{ref},r}(\theta_0)$ be the reference objective on that same dataset and at that same vector. Define

D_{m,r} = \operatorname{OBJ}_{m,r}(\theta_0) - \operatorname{OBJ}_{\mathrm{ref},r}(\theta_0).

A positive $D_{m,r}$ means the method's objective is larger than the reference objective under this sign convention. Summarize the replicate distribution of $D_{m,r}$, for example with its mean and spread. The mean can reveal systematic displacement, while a spread indicates how much that displacement varies across generated datasets. This is a comparison of approximations at the generating values, not a test of a scientific null hypothesis.

Worked Example: Summarizing Objective Differences

Suppose that, in a hypothetical pilot with four datasets, FO objective differences relative to the reference at $\theta_0$ are $3.2, 2.8, 3.6,$ and $2.4$. The corresponding FOCE differences are $0.8, 1.1, 0.5,$ and $0.8$, and Laplace differences are $0.1,-0.2,0.2,$ and $-0.1$. All values use the same $-2$ log-likelihood convention.

The mean FO difference is $(3.2+2.8+3.6+2.4)/4=3.0$ objective units. The mean FOCE difference is $(0.8+1.1+0.5+0.8)/4=0.80$ objective units. The mean Laplace difference is $(0.1-0.2+0.2-0.1)/4=0.00$ objective units. In these four illustrative datasets, FO is on average 3.0 units above the reference and FOCE 0.80 units above it, while the Laplace differences average to zero. The Laplace differences are not all zero, and four replicates cannot characterize their distribution reliably.

This calculation compares each method to the reference at common parameter values. Subtracting the minimized FO objective from the minimized Laplace objective is a different comparison: the methods may have been optimized at different estimates and may use different approximations. Such a difference cannot by itself be interpreted as approximation error at a fixed parameter vector.

Justification for Pairing and Monte Carlo Uncertainty

Pairing reduces noise in a direct method comparison when the estimates from the same dataset tend to move together. For methods $A$ and $B$, define $q_r=\hat{\theta}_{A,r}-\hat{\theta}_{B,r}$. The estimated average difference is $\bar{q}=R^{-1}\sum_r q_r$, with MCSE $s_q/\sqrt{R}$. The reason this can be more precise than comparing two independently simulated sets is the variance identity:

\operatorname{Var}(\hat{\theta}_{A}-\hat{\theta}_{B}) = \operatorname{Var}(\hat{\theta}_{A}) + \operatorname{Var}(\hat{\theta}_{B}) - 2\operatorname{Cov}(\hat{\theta}_{A},\hat{\theta}_{B}).

When estimates from the same dataset have positive covariance, the covariance term reduces the variability of their difference. The paired design does not guarantee a reduction in every setting; its advantage should be assessed from the replicate differences themselves. For objective differences, apply the same logic to paired $D_{m,r}$ values or paired differences between two methods' $D$ values.

Diagnostics, Alternatives, and Limits

A simulation summary is credible only if its generation, fitting, and reporting steps are checked. Record convergence status for every method and condition, inspect estimates at parameter boundaries, and report how many replicates failed or required a prespecified retry. Do not silently remove failed fits: if failures are more common in one method or scenario, summaries among successful fits alone can give a misleading impression of performance. Report failure rates and explain how they were handled.

Check the generated data as well as the fits. Confirm that the random effects have the intended variance, residuals follow the specified model, and the observation schedule is identical across methods. A coding error in the data-generating model can produce precise summaries of the wrong experiment. For each metric, report the number of replicates contributing to it and an MCSE or other appropriate indication of simulation precision. Increase $R$ if Monte Carlo uncertainty is too large to distinguish practically relevant differences.

A factorial design provides a direct comparison of conditions, while a one-factor-at-a-time design is simpler but can miss interactions. Here an interaction means, for example, that the effect of increasing $\omega^2$ differs between low and high $\kappa$. Replicate-level summaries by condition preserve that information. Averaging over all conditions first can conceal it.

Approach	Strength	Limitation for this comparison
Paired simulation across all methods	Isolates method differences on identical generated data.	Requires careful control of seeds, data, and fitting settings.
Separate datasets for each method	May be convenient when fitting pipelines cannot share data.	Method differences are mixed with random differences between generated datasets.
Comparison against a reference marginal likelihood	Measures objective displacement at a fixed parameter vector.	Reference accuracy must be verified; it does not alone summarize fitted-parameter bias.

The conclusions remain conditional on the chosen model, parameter values, design points, residual model, and fitting implementation. A result under a two-observation design does not automatically extend to rich sampling, different random-effect distributions, or correlated effects. Moreover, objective differences at the generating values do not directly determine how estimated parameters will differ after optimization. Parameter bias, objective differences, convergence, and Monte Carlo precision answer related but distinct questions.

Common pitfall: “Laplace has the smallest minimized objective, so it is closest to the true marginal likelihood.” This reasoning is not justified when each method minimizes a different approximate objective, potentially at a different parameter estimate. Instead, compare objectives at common parameter values against a consistently scaled reference, and separately summarize the fitted estimates against the known generating values.

Check Your Understanding

Use the definitions and comparison principles above. Brief solutions follow each question.

If a simulation has four sample-size and variance combinations and three curvature levels, how many conditions are there? With 300 replicates per condition, how many datasets are generated? Solution: There are $4\times3=12$ conditions and $12\times300=3{,}600$ datasets.
A parameter has generating value $0.40$ and average fitted value $0.46$. What are its estimated bias and relative bias? Solution: Bias is $0.46-0.40=0.06$; relative bias is $0.06/0.40=0.15$, or 15%.
Why fit each method to the same replicate dataset? Solution: It makes method differences paired and removes between-dataset variation from the direct comparison; the replicate differences still need to be summarized for their Monte Carlo uncertainty.
FO and Laplace minimized objectives differ by 2.5 units. Can this difference be treated as Laplace's objective error? Solution: No. The objectives may be evaluated at different fitted parameter values and are based on different approximations. Evaluate both at a common parameter vector against a suitable reference to assess objective displacement.
What should be reported if FO fails to converge in 8 of 100 replicates while Laplace fails in 1? Solution: Report both failure rates, the prespecified handling of failed fits, and the number of replicates used in each summary. Do not hide failures by reporting only successful-fit estimates.