From a Model Specification to an Estimation Plan
Choosing estimation settings is not a search for one universally best optimizer configuration. It is a set of linked decisions: what values give the model a plausible starting point, which method can approximate the likelihood for this design, and how the numerical algorithm will judge that it has made adequate progress. These choices should follow from the dosing and sampling scheme, the expected variability, and the parameters the data can inform.
Consider a hypothetical study of an intravenously administered drug with known dose and sampling times. The structural model has a central and a peripheral compartment, linear elimination from the central compartment, and distribution between compartments. For an individual, let \(A_c(t)\) and \(A_p(t)\) be the amounts of drug in the compartments. The model is
Here \(CL\) is clearance, \(Q\) is intercompartmental clearance, \(V_c\) is central volume, and \(V_p\) is peripheral volume. Each is positive. For subject \(i\), a common population model is \(P_i=P_{\mathrm{typ}}\exp(\eta_{Pi})\) for a parameter \(P\), where \(\eta_{Pi}\) is a normally distributed individual random effect with mean zero and variance \(\omega_P^2\). Conditional on an individual’s parameters, suppose concentration observations have independent normal residual errors with mean \(C_{ij}\) and variance \(\sigma_{\mathrm{add}}^2+(\sigma_{\mathrm{prop}}C_{ij})^2\). The marginal likelihood integrates the conditional observation likelihood over the random-effects distribution.
Use the Design to Choose Starting Values
Starting values are hypotheses about a plausible region of parameter space, not estimates that the data have already established. In a two-compartment model, early observations help characterize the rapid distribution phase and central volume; later observations inform the terminal phase and clearance. A design with few early samples may not support a precise central volume, while one with little late follow-up may give weak information about the terminal slope. Choose starts that are physiologically plausible and consistent with the observations, but do not treat a visual fit or a single summary as proof of identifiability.
Worked Example: Constructing Structural Starts From a Concentration Curve
Suppose a hypothetical IV bolus study has dose 100 mg. A preliminary summary of the concentration-time profile is \(C(t)=8e^{-1.0t}+2e^{-0.10t}\) mg/L, with time in hours. Use this biexponential summary to construct approximate starts for \(V_c\), \(CL\), \(V_p\), and \(Q\). Treat the summary as an approximation for initialization, not as a substitute for fitting the population model.
At time zero, \(C(0)=8+2=10\) mg/L. For an IV bolus, the initial concentration is dose divided by central volume, so
The area under the curve for this biexponential profile is
Thus, the dose-to-AUC clearance approximation is \(CL\approx100/28=3.571\) L/h. The mean residence time from the same profile is the first moment divided by AUC:
For this linear IV model, \(V_{ss}=CL\times MRT=3.571\times7.429=26.53\) L, so \(V_p\approx V_{ss}-V_c=16.53\) L. To obtain a compatible \(Q\), use the relationship between the sum of the two disposition rate constants and the microconstants:
Substitution gives \(Q=(1.10-3.571/10.0)/(1/10.0+1/16.53)=4.63\) L/h. The resulting starts are therefore \(CL=3.57\) L/h, \(V_c=10.0\) L, \(Q=4.63\) L/h, and \(V_p=16.5\) L. They are mutually consistent with the approximate curve: the product of its two rate constants is \(1.0\times0.10=0.10\ \mathrm{h}^{-2}\), and \(CLQ/(V_cV_p)\approx0.100\ \mathrm{h}^{-2}\). The data design must still provide enough early and late information for the fitted model to refine these values.
Population variability starts need the same care. For a lognormal parameter with coefficient of variation \(CV\), the corresponding log-scale variance is \(\omega^2=\ln(1+CV^2)\). This follows because a lognormal variable with log-scale variance \(\omega^2\) has squared coefficient of variation \(e^{\omega^2}-1\). The conversion is useful when prior knowledge or a pilot summary is expressed as a coefficient of variation rather than a variance.
Worked Example: Converting Variability Summaries to Random-Effects Starts
Suppose hypothetical pilot information suggests a 35% coefficient of variation for clearance and a 30% coefficient of variation for central volume. For lognormal random effects, calculate the variance and standard deviation starts on the eta scale.
For clearance, \(CV=0.35\), so \(\omega_{CL}^2=\ln(1+0.35^2)=\ln(1.1225)=0.1155\), and \(\omega_{CL}=\sqrt{0.1155}=0.3399\). For central volume, \(CV=0.30\), so \(\omega_{V_c}^2=\ln(1.09)=0.08618\), and \(\omega_{V_c}=\sqrt{0.08618}=0.2936\). In a model that estimates variances, the starting values are approximately 0.116 and 0.0862; in a model parameterized by eta standard deviations, they are approximately 0.340 and 0.294. These pilot summaries justify starts, not a decision that both variances must be estimable in the current study. If the design is sparse, begin with only variability components that the data can reasonably support.
Match the Estimation Method to the Likelihood
The population likelihood is generally not a simple normal likelihood in the fixed effects: individual random effects are unobserved and must be integrated out. A method that approximates this integration must be selected with the model’s nonlinearity and the amount and timing of information per subject in mind. For a reasonably informative concentration-time design and a smooth two-compartment model, first-order conditional estimation with interaction (FOCEI) is a defensible starting choice. It approximates the individual random-effects contribution around a conditional mode and accounts for the dependence of residual variance on individual predictions when that interaction is included.
The choice is conditional, not automatic. A first-order method may be an option for a model whose response is nearly linear in random effects, but can be less reliable when that approximation is poor. Laplace approximation also uses local information around a conditional mode and may be useful for suitable likelihoods. Stochastic approximation expectation-maximization (SAEM) can be attractive when the likelihood is strongly nonlinear or local approximations are unsatisfactory, but it introduces simulation settings and Monte Carlo variability that must be controlled. These methods are alternatives for approximating the same model likelihood; they are not interchangeable guarantees of a correct answer.
| Method or approach | When it may be reasonable | Consideration for this model |
|---|---|---|
| FOCEI | Nonlinear mixed-effects models with informative individual data and a smooth conditional likelihood | A practical initial choice for a two-compartment model with adequate sampling and proportional-error interaction |
| Laplace approximation | Likelihoods well represented locally around each individual’s conditional mode | Compare if supported by the software and model; local approximation quality remains relevant |
| SAEM | Models where integration is difficult or nonlinearities make local approximations concerning | Can provide a useful alternative, but requires adequate simulation and reproducible Monte Carlo settings |
| First-order approximation | Settings where response changes are sufficiently close to linear in random effects | May be less dependable when distribution-phase curvature or prediction-dependent residual variance is important |
A compact justification is therefore: with rich early and late samples, moderate expected variability, and a smooth two-compartment concentration model, start with FOCEI and the stated residual interaction. If the design is very sparse, or the likelihood is especially curved or difficult to integrate, plan a method comparison rather than assuming the first choice is adequate. Keep the model, data, and parameterization fixed during that comparison; otherwise changes in the objective function cannot be attributed to the estimation method alone.
Scale Parameters and Set Joint Stopping Criteria
The optimizer’s coordinates should make a numerical step meaningful across parameters with different units and magnitudes. Positive parameters can be represented on a log scale, as in \(z_j=\ln(\theta_j/\theta_{j,\mathrm{ref}})\), where \(\theta_{j,\mathrm{ref}}\) is a fixed positive reference value. A unit change in \(z_j\) then corresponds to a multiplicative change in the original parameter, rather than the same absolute change in liters, liters per hour, and variance. For parameters not transformed, use explicit scaling so that typical coordinate changes have comparable numerical size. The scale should reflect plausible variation, not be chosen merely to make the printed values look similar.
Stopping should require evidence from more than one measure of progress. A small objective-function change alone may occur on a flat surface while parameters continue moving; a small parameter step alone may occur because the step has been damped. Define criteria in advance and require them jointly, subject to the optimizer’s valid termination rules. For example, with \(N\) observations and objective function \(OFV\), an implementation-independent set of targets could require, for two consecutive iterations,
These are illustrative targets, not universal software defaults or proof that a model is correct. The gradient condition is stated per observation to make its scale less dependent on sample size; software may define and report gradients differently, so map the intended criterion to its documented convention. If the method uses stochastic simulation, also specify settings that make the approximation sufficiently stable for the intended comparison. Tightening tolerances is informative only if the model and objective remain unchanged and the resulting estimates and objective are compared consistently.
Worked Example: Applying Scaled Coordinates and Tolerances
Suppose the starts from the concentration summary are \(CL_{\mathrm{ref}}=3.57\) L/h, \(V_{c,\mathrm{ref}}=10.0\) L, \(Q_{\mathrm{ref}}=4.63\) L/h, and \(V_{p,\mathrm{ref}}=16.5\) L. After an optimization iteration, \(CL=3.60\) L/h and \(V_c=10.1\) L. Using \(z=\ln(\theta/\theta_{\mathrm{ref}})\), calculate the scaled coordinate changes from the starts.
For clearance, \(z_{CL}=\ln(3.60/3.57)=0.00837\); for central volume, \(z_{V_c}=\ln(10.1/10.0)=0.00995\). The changes are thus roughly 0.84% and 1.00% in multiplicative terms. They are both much larger than a \(10^{-4}\) relative step target, so the optimizer has not met the illustrative parameter-step condition at this iteration, regardless of a small objective change. At a later iteration, if the maximum scaled step is \(6.0\times10^{-5}\), the objective changes by 0.004 for \(N=200\), and the maximum absolute gradient component is 0.12, then the three checks are: step \(6.0\times10^{-5}<10^{-4}\); objective ratio \(0.004/(200+|OFV|)\leq0.004/200=2.0\times10^{-5}\), which is not below \(10^{-6}\); and per-observation gradient \(0.12/200=6.0\times10^{-4}<10^{-3}\). The joint stopping rule is not met because the objective criterion fails. The example shows why checking one criterion cannot substitute for checking the others.
Preflight Checks, Limitations, and a Common Pitfall
Before launching a long estimation, verify that the model returns finite, positive predictions at the observed times for the proposed starts; that dose, time, and concentration units agree; and that the parameter transformations and bounds are consistent with the scientific model. Simulate or calculate representative individual concentration profiles to check that the starts place observations in a plausible range. Confirm that the sampling design contains information relevant to each intended random-effects component. If multiple initial values are plausible, use the earlier material on sensitivity to starting values and layered estimation strategies to plan a controlled comparison. Record the method, initial values, scaling, tolerances, and any integration settings so the run can be reproduced.
A common mistake is to respond to an optimizer that stops quickly by tightening every tolerance or increasing every iteration limit. That reasoning assumes the difficulty is insufficient numerical effort. The correction is to separate the questions: are the starts plausible and internally consistent; is the chosen likelihood approximation suitable; are the coordinates well scaled; and are the stopping criteria being assessed jointly? More iterations cannot compensate for incompatible units, unsupported variability components, or an unsuitable likelihood approximation. Conversely, a converged run is not evidence by itself that the structural model or its assumptions are adequate.
The estimates and numerical behavior also have limits that settings cannot remove. A two-compartment model may be weakly informed by a design lacking early distribution samples or late terminal samples. Random-effects variances can be difficult to estimate when there are few observations per subject or limited between-subject information. Starting values derived from a summary curve ignore the population distribution and residual error, so they are initialization aids only. Method comparisons and tighter tolerances can assess numerical sensitivity, but they do not establish scientific validity or make unsupported parameters identifiable.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- A design has several early samples but only one sample after the distribution phase. Which aspect of the two-compartment model deserves particular caution when choosing starts, and why?
- A lognormal random effect has a coefficient of variation of 50%. Calculate its eta-scale variance.
- Why is a small change in the objective function not, by itself, sufficient evidence that an optimizer has converged?
- For an adequately sampled, smooth two-compartment model with proportional residual interaction, name a reasonable initial estimation method and one circumstance that could motivate comparing an alternative.
Solutions. (1) The terminal slope and parameters informed by late concentrations, including clearance-related behavior, deserve caution because one late observation gives limited information about the terminal phase. (2) \(\omega^2=\ln(1+0.50^2)=\ln(1.25)=0.2231\). (3) The objective can change little while parameters still move along a flat direction, or while the gradient remains appreciable; parameter-step and gradient checks should also be considered. (4) FOCEI is a reasonable initial choice under these conditions. A strongly nonlinear likelihood or concern that a local approximation is inadequate could motivate a controlled comparison with an alternative such as SAEM or Laplace approximation.