Why Look for Patterns, Not Just One Difference?
FO, FOCE, and Laplace replace or approximate the integration over individual random effects in an NLME likelihood. Because the approximations can differ, they may produce different parameter estimates on the same data. A difference alone, however, does not establish that one approximation is biased or that it changes the scientific conclusion. It could also reflect ordinary sampling variability, a poorly optimized objective, or a mismatch between the fitted and data-generating models.
A useful diagnostic therefore combines three kinds of evidence: parameter shifts in simulations with known generating values, objective behavior within each approximation, and stability across optimizer starting values. The strongest concern arises when estimates shift consistently across simulated datasets, the shifts are large enough to matter in context, and repeated fits show that the differences are not simply caused by unstable optimization.
For an NLME model, the conditional likelihood contribution for subject \(i\) is integrated over that subject’s random effects. The marginal likelihood is the product of those contributions. An approximation replaces the generally difficult integral with an approximate contribution; fitting then maximizes the resulting approximate likelihood. In what follows, \(A\) denotes the chosen approximation and \(J_A\) denotes its objective, conventionally \(-2\) times its approximate log likelihood, so smaller values indicate a better fit within that same approximation.
Here \(y_i\) is subject \(i\)’s data, \(\eta_i\) is the vector of random effects, \(\theta\) represents fixed effects, and \(\Omega\) and \(\sigma\) represent random-effect and residual-error parameters. The notation suppresses the simultaneous estimation of all model parameters. The comparison is meaningful only when the approximations are applied to the same model specification, data, parameter constraints, and estimation target.
Simulation-Based Parameter Shifts
A simulation study can expose systematic estimation patterns because its generating parameters are known. To assess approximation behavior, generate datasets under a specified NLME model, fit each dataset using the candidate approximations, and summarize estimates across replications. If the goal is to compare methods rather than designs, use the same simulated dataset for each method in a replication. This paired comparison removes some dataset-to-dataset variation from the method contrast.
For \(R\) successful fits, the estimated bias is \(\bar{\hat{\theta}}_A-\theta_0\). If the sample standard deviation of the \(R\) estimates is \(s_A\), the Monte Carlo standard error of their mean is \(s_A/\sqrt{R}\). This describes simulation noise in the estimated mean, not the standard error that a single fitted model would report for its parameter estimate. For a positive parameter, relative bias is often summarized as \(100(\bar{\hat{\theta}}_A-\theta_0)/\theta_0\) percent; the signed absolute bias and the parameter’s units should also be retained.
When comparing two approximations on the same replications, calculate a paired difference for each dataset, \(d_r=\hat{\theta}_{A,r}-\hat{\theta}_{B,r}\). The mean of these differences describes the average shift, and its Monte Carlo standard error is \(s_d/\sqrt{R}\). The difference includes any relative finite-sample bias of the two procedures; it does not, by itself, isolate an approximation error from every other source of estimator behavior.
Worked Example: A Persistent Clearance Shift in Simulation
Hypothetical setup. Suppose 200 datasets are generated from the same correctly specified concentration model, with true typical clearance \(CL_0=5.00\) L/h. Each dataset is fitted by FO, FOCE, and Laplace using identical model structure and parameter constraints. Assume all 200 fits converge and that their mean estimates and across-replication standard deviations are as follows.
| Approximation | Mean estimate (L/h) | Across-replication SD (L/h) |
|---|---|---|
| FO | 5.62 | 1.40 |
| FOCE | 5.18 | 1.10 |
| Laplace | 5.09 | 1.08 |
For FO, the estimated bias is \(5.62-5.00=0.620\) L/h, or \(100(0.620/5.00)=12.4\%\). For FOCE, it is \(5.18-5.00=0.180\) L/h, or \(3.60\%\). For Laplace, it is \(5.09-5.00=0.0900\) L/h, or \(1.80\%\).
The Monte Carlo standard errors are \(1.40/\sqrt{200}=0.0990\) L/h for FO, \(1.10/\sqrt{200}=0.0778\) L/h for FOCE, and \(1.08/\sqrt{200}=0.0764\) L/h for Laplace. The FO mean’s upward shift is large relative to its Monte Carlo error in this simulation. FOCE and Laplace are closer to the generating value, although their estimated biases are not exactly zero.
Suppose the paired FO-minus-Laplace differences have a mean of \(5.62-5.09=0.530\) L/h and a standard deviation of \(0.720\) L/h. Their Monte Carlo standard error is \(0.720/\sqrt{200}=0.0509\) L/h. This supports a stable difference between the two procedures under this simulation design. It does not prove that the same difference will occur under another design or model. The inference depends on the assumed generating model, sampling schedule, error model, parameter values, and successful, comparable fitting.
The generating assumptions deserve explicit attention. For the example to support an approximation-focused interpretation, the simulated model should match the fitted structural and error models, and the intended population and sampling design should be represented. If the fitted model differs from the generating model, observed shifts combine approximation behavior with model misspecification. If some fits fail, reporting only successful fits can also distort the comparison; convergence rates and reasons for failure should be summarized by method.
Reading Objective Changes Carefully
Objective values help diagnose whether an optimizer has made progress and whether different starting values lead to a stable solution. Their interpretation must be kept within the method that produced them. FO, FOCE, and Laplace objectives approximate the likelihood differently, so a lower reported objective under one approximation is not automatically evidence that it is more accurate or that it fits the data better on a common likelihood scale. In particular, do not treat a difference across approximations as a likelihood-ratio statistic.
Worked Example: Comparing Objective Movement Within Each Method
Hypothetical setup. On one dataset, FOCE starts from \(CL=4.5\) L/h with objective \(J_{\mathrm{FOCE}}=428.0\) and ends at \(CL=5.6\) L/h with \(J_{\mathrm{FOCE}}=410.7\). Laplace starts from \(CL=4.5\) L/h with \(J_{\mathrm{Laplace}}=423.1\) and ends at \(CL=5.2\) L/h with \(J_{\mathrm{Laplace}}=408.9\). Assume each objective is evaluated and minimized according to its own approximation.
Within FOCE, the objective decreases by \(428.0-410.7=17.3\). Within Laplace, it decreases by \(423.1-408.9=14.2\). These reductions show improvement from the respective starting points, but they do not establish that either final fit is globally optimal.
The final estimates differ by \(5.6-5.2=0.4\) L/h, which may be important if it changes a dosing or exposure conclusion. The final objective values differ by \(410.7-408.9=1.8\), but that cross-approximation subtraction is not a valid test of whether Laplace is preferable: the objectives come from different approximate likelihoods. To interpret the parameter shift, inspect simulation behavior and start-value stability as well as the scientific size of the change.
A stable reduction in an objective is reassuring about progress, not a complete convergence diagnosis. Check the optimizer’s convergence status, parameter gradients or other available convergence measures, boundary estimates, and whether the final values reproduce from alternative starts. If objective changes are negligible but gradients remain substantial, the reported solution may not satisfy the intended stopping criterion. If estimates or objective values change markedly across starts, numerical or local-optimization problems remain plausible.
Starting-Value Sensitivity and Failure Modes
Repeated initialization is a practical diagnostic. Fit the same dataset with the same approximation, model, and constraints, but with dispersed, plausible starting values. A stable pattern consists of estimates that agree to a meaningful precision and objective values that are essentially the same, with acceptable convergence diagnostics. Disagreement can indicate local optima, a flat or poorly identified direction, parameter-boundary behavior, or an optimizer that has not converged. It is not, on its own, evidence of approximation bias.
Worked Example: A Competing Laplace Solution
Hypothetical setup. A Laplace fit is repeated three times on the same data and under the same parameter constraints. Two starts give \(CL=4.81\) L/h and \(J_{\mathrm{Laplace}}=408.7\), with acceptable convergence checks. A third gives \(CL=6.30\) L/h and \(J_{\mathrm{Laplace}}=415.9\), also marked as converged.
The third run’s objective exceeds the lower value by \(415.9-408.7=7.2\). Since Laplace is minimizing this same objective in all three runs, the results do not support treating the \(CL=6.30\) L/h solution as equivalent to the two \(CL=4.81\) L/h solutions. The difference warrants checking initialization, gradients, parameter boundaries, and whether the optimizer has found a local rather than a better solution. A “converged” label alone does not settle the issue.
If repeated starts and adequate optimization diagnostics consistently support the lower-objective solution, the starting-value sensitivity may be resolved as an optimization failure rather than an approximation-induced estimate shift. If the issue remains unresolved, interpretation of the Laplace estimate is uncertain. Only after establishing stable fitting should a difference between Laplace and FOCE, for example, be interpreted as approximation sensitivity—and then it should be checked across simulated datasets rather than inferred from this single dataset.
Hold the dataset, structural and error models, constraints, and estimation target fixed when comparing approximations. For simulation, also document the generating model and sampling design.
Record convergence, objective value, parameter estimates, boundary behavior, and relevant gradient or optimizer diagnostics. Investigate failed fits rather than silently excluding them.
Within each approximation, look for reproducible estimates and objective values. Resolve or report material instability before attributing a parameter shift to approximation choice.
Report mean estimates, bias relative to the generating values, Monte Carlo uncertainty, and paired method differences. Ask whether the shift changes a scientifically relevant interpretation.
Comparing Evidence and Its Limits
| Diagnostic pattern | What it can indicate | What it cannot establish alone |
|---|---|---|
| A consistent parameter shift across simulated datasets | Approximation-dependent estimation behavior under the specified simulation design | That the shift will generalize to other designs or models |
| Large objective reduction from an initial value | The optimizer improved its objective within that approximation | Global optimality or superiority to another approximation |
| Different estimates across starts within one approximation | Possible local optimization, weak identification, or boundary issues | Approximation bias as the sole cause |
| Similar estimates across approximations in one dataset | Little observed sensitivity for those estimates and that dataset | Absence of bias in other parameters, data designs, or repeated samples |
The main alternative to a simulation-based diagnostic is to inspect one fitted dataset and compare estimates and within-method objective behavior. That can flag a concern quickly, but it cannot show whether a shift is systematic across repeated samples. Simulation can estimate such patterns only for the generating conditions actually examined. A trusted reference likelihood calculation, when available, offers another comparison target, but agreement with a reference does not remove the need to check model adequacy and optimizer behavior.
Several limitations matter in practice. Simulation results depend on the chosen parameter values and design; behavior under sparse sampling or strong nonlinearity need not match behavior under richer data. A finite number of replications gives noisy bias estimates, especially when estimator variability is high. Large parameter shifts may be statistically clear yet clinically negligible, while smaller shifts may change a threshold-based conclusion. Report parameter units and connect the magnitude to the decision or scientific interpretation at issue.
Check Your Understanding
Answer each question, then compare with the short solutions below.
- In a simulation, the generating value is \(2.00\) mg/L and the mean estimate from 100 replications is \(2.16\) mg/L. What is the estimated bias and relative bias?
- Why is it generally inappropriate to subtract a FOCE objective from a Laplace objective and interpret the result as a likelihood-ratio statistic?
- One approximation produces different parameter estimates from three starting values, with materially different objective values. What should be checked before describing the pattern as approximation bias?
- What does the Monte Carlo standard error \(s_A/\sqrt{R}\) quantify in a simulation study?
Solutions. (1) The estimated bias is \(2.16-2.00=0.160\) mg/L; the relative bias is \(100(0.160/2.00)=8.00\%\). (2) FOCE and Laplace use different approximate likelihood objectives, so their objective difference is not generally a likelihood-ratio statistic. (3) Check convergence diagnostics, gradients, parameter boundaries, and repeated starts; determine whether a stable solution can be obtained within that approximation before interpreting between-method differences. (4) It estimates the simulation-to-simulation uncertainty in the mean estimate across the \(R\) replications, not the standard error for a parameter estimate from one dataset.