Why the Model and Design Should Guide the Method
An estimation method is not a neutral switch that changes only computation time. Methods differ in how they approximate the likelihood contributed by unobserved individual random effects. Their accuracy can therefore depend on the shape of the model, the amount and timing of data within each individual, the size of the random-effects variability, and the residual-error model. A method that is adequate for a nearly linear model with modest variability may be less reliable when predictions are strongly nonlinear in the random effects.
The practical task is to make a defensible choice before fitting, then check whether that choice behaves adequately for the data at hand. A large total sample does not by itself guarantee that individual random effects are well informed: many individuals with only one or two observations can still provide little within-individual information. Conversely, dense sampling does not make every approximation accurate if the likelihood is sharply curved, asymmetric, or otherwise poorly represented by a local approximation.
A common model statement is \(y_{ij}=f(t_{ij},\theta,\eta_i)+\epsilon_{ij}\), where \(y_{ij}\) is the observation for individual \(i\) at time \(t_{ij}\), \(f\) is the structural model, and \(\epsilon_{ij}\) is residual error. A standard random-effects assumption is \(\eta_i\sim N(0,\Omega)\), with \(\Omega\) the random-effects covariance matrix. The precise conditional distribution of \(\epsilon_{ij}\) depends on the chosen residual-error model. Under these assumptions, the individual likelihood has the form
where \(y_i\) denotes all observations for individual \(i\), and \(\psi\) denotes residual-error parameters. Estimation methods differ in how they handle this integral. The earlier material on FO, FOCE, and Laplace establishes the individual approximations; the selection decision here is how those approximations fit the model and design.
What the Main Methods Approximate
First-order estimation (FO) linearizes the model in the random effects around their population mean, usually \(\eta_i=0\). This can be computationally economical, but the approximation may be poor when the model changes substantially over the range of plausible random effects. FO does not condition its expansion on each individual’s data in the way FOCE does.
First-order conditional estimation (FOCE) expands around an individual-specific conditional estimate of the random effects. It can better represent nonlinear relationships when the individual data locate that conditional estimate reasonably well. FOCE with interaction (FOCEI) additionally accounts for dependence between the random effects and residual variance when the residual model makes that variance depend on the individual prediction. The interaction option is relevant only when the residual-error structure supports it; it is not a general-purpose accuracy setting.
The Laplace approximation uses a local quadratic approximation to the log integrand around its conditional mode. It uses local curvature as well as the mode location, and can be useful for nonlinear models when the conditional distribution around that mode is well represented by a Gaussian shape. It remains a local approximation: a poor mode, strong asymmetry, or multiple important regions can undermine it.
The Stochastic Approximation Expectation-Maximization method (SAEM) uses simulation and stochastic approximation to estimate model parameters while accounting for latent random effects. It avoids relying on the same local likelihood expansion as FOCE or Laplace, which can make it attractive for difficult nonlinear models. Its practical cost includes stochastic variability and the need for adequate simulation and run settings; Monte Carlo precision and stability must be assessed rather than assumed.
Why Local Approximation Quality Depends on Design
Sparse sampling matters because a conditional estimate of \(\eta_i\) is determined by an individual’s observations relative to the population model. With few observations, the data may not locate that estimate sharply. In that case, a method centered on the conditional estimate can still depend heavily on the local approximation around a weakly informed location. Richer sampling can sharpen individual information, but only if the observations occur at times that distinguish relevant model features. Several measurements in an uninformative region do not necessarily identify clearance, absorption, or a turnover rate.
Random-effects variability also changes the range over which a local approximation must represent the model. For a lognormally distributed individual parameter with coefficient of variation \(CV\), the corresponding normal random-effect variance is \(\omega^2=\ln(1+CV^2)\). Greater variability means that more of the random-effects distribution lies farther from zero, where a Taylor approximation around zero may be less representative for a nonlinear structural model.
Residual error is another design-and-model consideration. In a proportional-error model, the residual variance depends on the prediction. If the prediction itself changes with \(\eta_i\), then residual variance changes with \(\eta_i\) as well. FOCEI addresses this interaction within its approximation. If residual variance is constant and independent of the random effects, interaction is not the reason to choose FOCEI.
A Short Derivation of the Approximation Trade-Off
Write the individual log integrand as \(g_i(\eta_i)=\log p(y_i\mid\eta_i,\theta,\psi)+\log p(\eta_i\mid\Omega)\). The individual likelihood integrates \(\exp(g_i(\eta_i))\) over \(\eta_i\). FO uses a first-order expansion of the model around zero random effects. For a scalar random effect, for example,
This replaces the nonlinear dependence on \(\eta_i\) with a local linear dependence around zero. If the second derivative is large over plausible values of \(\eta_i\), the omitted curvature can matter. FOCE instead centers its local approximation on a conditional estimate \(\hat{\eta}_i\), which is informed by that individual’s data. A quadratic approximation of the log integrand at a conditional mode \(\hat{\eta}_i\) gives the Laplace form
where \(q\) is the dimension of \(\eta_i\), and the approximation requires a well-defined local maximum with negative-definite \(H_i\). The determinant term reflects local curvature: a wider peak contributes more integrated likelihood than a narrower peak of the same height. This derivation explains the trade-off, not a universal ranking. If the integrand is poorly represented by a quadratic near its mode, Laplace can still be inaccurate; if the model is nearly linear around zero, FO may be adequate.
Worked Examples of Method Selection
Worked Example: Sparse Concentrations and High Clearance Variability
Consider a hypothetical oral pharmacokinetic analysis with 120 individuals, each contributing two concentrations: one at 1 hour and one at 8 hours. The structural model is nonlinear in clearance and absorption parameters, and the expected between-individual coefficient of variation in clearance is 80%. Choose a primary method and explain the trade-off.
Design calculation. The total number of observations is \(120\times2=240\), but each individual has only 2 observations. The relevant limitation is not the total 240 alone: with two measurements, an individual’s data may not distinguish several individual parameters or locate the conditional random-effects estimate precisely.
Variability calculation. For lognormal clearance variability, \(CV=0.80\), so
The substantial spread and nonlinear dependence make a first-order expansion around \(\eta_i=0\) a potentially consequential approximation. FO is therefore not the default choice without evidence that its results are adequate. A simulation-based method such as SAEM is a defensible primary candidate because it does not rely on the same local expansion, provided the stochastic estimation is run and checked adequately. FOCEI or Laplace can be useful comparisons, but sparse data may leave conditional locations weakly informed, and Laplace additionally depends on a good local quadratic representation.
Decision and checks. I would begin with SAEM, then examine stability across suitable runs and compare parameter estimates and predictive diagnostics with a likelihood-approximation fit if feasible. The choice is not a claim that SAEM must win: convergence behavior, model identifiability, and fit diagnostics still matter. Sparse observations also limit what any method can recover from the design.
Worked Example: Prediction-Dependent Residual Variance
In a hypothetical rich-sampling study, an individual has 10 observations spanning a wide concentration range. The model uses proportional residual error, with residual standard deviation \(\sigma_{\mathrm{prop}}f_{ij}\), where \(f_{ij}\) is the individual prediction. Let \(\sigma_{\mathrm{prop}}=0.15\). Compare residual variance at predictions of 5 and 20 concentration units, and choose between FOCE and FOCEI.
At prediction 5, the residual standard deviation is \(0.15\times5=0.75\) concentration units, giving variance \(0.75^2=0.5625\) concentration-units squared. At prediction 20, it is \(0.15\times20=3.00\) concentration units, giving variance \(3.00^2=9.00\) concentration-units squared. Thus the variance ratio is \(9.00/0.5625=16.0\).
Because the individual prediction depends on the random effects, this residual variance also depends on the random effects. FOCEI is the more directly matched candidate: it accounts for the interaction between random effects and prediction-dependent residual variance. Ordinary FOCE would not represent that interaction in the same way. The 10 observations per individual may help locate conditional estimates, although their timing and information content still need review. A Laplace fit could be compared if the conditional log integrand is suitably well behaved, but it does not replace the need to specify the residual model correctly.
The arithmetic supports choosing FOCEI for the stated error structure; it does not establish that FOCEI will fit well. Residual plots and predictive checks should still be inspected for systematic error, and an interaction term should not be added merely because it is available.
Worked Example: A Rich, Nearly Linear Model With Modest Variability
Suppose a hypothetical study has 60 individuals with 12 well-spaced observations per individual. Over the fitted parameter range, the structural model is approximately linear in a single random effect, residual variance is constant, and the random-effects coefficient of variation is 20%. Compare FO with the more elaborate candidates.
There are \(60\times12=720\) observations, or 12 per individual. The design is relatively rich within individuals. For the lognormal random effect,
This corresponds to a comparatively modest random-effects spread. Together with approximate linearity, it makes the FO expansion around zero more plausible than in a strongly nonlinear, high-variability setting. Since residual variance is constant and independent of the prediction, the interaction feature of FOCEI is not specifically needed. FO is therefore a reasonable economical candidate, with FOCE as a useful comparison if computationally practical.
This is a provisional selection, not proof that FO is accurate. If FO and FOCE give materially different parameter estimates or predictions, or if residual and prediction diagnostics show a systematic pattern, the approximation may not be adequate despite the favorable design summary. A more computationally demanding method is justified when the difference affects the scientific conclusions, not simply because it is available.
A Practical Selection and Review Workflow
Method selection is best documented as a reasoned sequence. The aim is to make the assumptions visible, choose a method consistent with them, and specify what evidence would prompt reconsideration.
Identify important nonlinearities, the random-effects structure, and whether residual variance depends on individual predictions.
Record the number and timing of observations per individual. Ask whether they inform the individual parameters and random effects of interest, rather than relying only on the total sample size.
Consider the plausible random-effects spread, curvature, and whether the conditional likelihood is likely to be well represented by a local quadratic approximation.
Choose a method whose assumptions fit the model and design. When feasible, compare with a materially different approximation or simulation-based approach.
Check estimation behavior, parameter plausibility, residual and predictive diagnostics, and whether method comparisons change conclusions. Treat unexplained disagreement as a reason to investigate.
Diagnostics, Alternatives, and Limits
Warning signs include markedly different parameter estimates across plausible methods, unstable results across stochastic runs, implausible estimates, poor predictive checks, or systematic residual patterns. These signs do not uniquely diagnose an approximation problem: misspecified structural or residual models, weak identifiability, and inadequate data can produce similar symptoms. Diagnose the model and design alongside the method rather than changing algorithms until one produces a preferred answer.
| Method | Potentially favorable setting | Main approximation or practical concern |
|---|---|---|
| FO | Approximately linear random-effects relationship and modest variability | Expansion around zero can be poor with substantial curvature or broad variability |
| FOCE | Nonlinear model with useful individual information for conditional estimates | Conditional approximation may be weak when individual data are sparse or uninformative |
| FOCEI | FOCE-type setting with prediction-dependent residual variance | Interaction is not a remedy for misspecification or a universal improvement |
| Laplace | Conditional integrand well represented by a local quadratic near its mode | Can be inaccurate with poor modes, strong asymmetry, or multiple important regions |
| SAEM | Nonlinear or otherwise challenging random-effects likelihood where simulation is suitable | Stochastic variability and computational demands require adequate run settings and stability checks |
No method can compensate for a design that does not identify the parameters of interest, a misspecified model, or insufficient information. Nor does a method comparison prove which estimate is correct: it reveals sensitivity to approximation choices. In practice, the preferred method is the one justified by the model and design, with acceptable stability and diagnostics, and with limitations reported clearly.
Common Pitfall
Wrong reasoning: “There are 1,000 observations in the study, so FOCE will have well-informed conditional estimates for every individual.” This confuses total sample size with information per individual. For example, 1,000 observations from 500 individuals still means only 2 observations per individual on average, and those observations may not identify the relevant individual effects.
Correction: Examine the number, timing, and informativeness of observations within individuals, together with model nonlinearity and variability. Then treat the method choice as a hypothesis to check through stability, comparisons, and diagnostics—not as a guarantee based on sample size alone.
Check Your Understanding
Answer each question, then compare with the concise solutions below.
- A nonlinear model has few observations per individual and substantial random-effects variability. Why might FO be a questionable default, and what alternative could be considered?
- A model has constant residual variance independent of individual predictions. Does this fact alone favor FOCEI over FOCE? Explain.
- For a lognormal random effect with \(CV=0.50\), calculate \(\omega^2\) and \(\omega\) to three significant figures.
- State one condition that makes a Laplace approximation more credible and one failure mode that can weaken it.
- Two methods give different estimates but both runs appear numerically stable. What should the analyst do before choosing one?
Solutions. (1) High variability and nonlinearity can make a linear expansion around zero inaccurate; sparse data may also give weak conditional estimates. SAEM is a candidate that avoids that same local expansion, subject to stochastic-stability checks. (2) No. Constant residual variance gives no specific reason for the interaction adjustment; method choice still depends on nonlinearity, design, variability, and diagnostics. (3) \(\omega^2=\ln(1+0.50^2)=\ln(1.25)\approx0.223\), so \(\omega\approx\sqrt{0.223}=0.472\). (4) A well-defined conditional mode with a locally near-quadratic log integrand supports Laplace; strong asymmetry or multiple important regions can weaken it. (5) Compare the model assumptions and approximation trade-offs, inspect parameter plausibility and predictive diagnostics, and investigate whether the discrepancy changes scientific conclusions. Stability alone does not establish accuracy.