Why Estimate a Complex Model in Stages?
A complex pharmacometric model may contain typical structural parameters, interindividual variability, correlations among random effects, residual-error parameters, and covariate effects. Estimating all of these from weak or arbitrary initial values can make optimization difficult: one parameter block may compensate for another before either is well informed by the data. A staged strategy gives the optimizer a more organized route to a usable solution.
Staging is a computational workflow, not a different statistical model or a substitute for checking whether the model is appropriate. Intermediate estimates are conditional on the parameters held fixed in that stage. They are useful as stepping stones, but the final estimates should generally come from a fit in which all intended parameters are estimated together.
For a model with typical structural parameters \(\theta\), random-effects covariance matrix \(\Omega\), residual parameters \(\sigma\), and covariate coefficients \(\beta\), write the estimation objective as \(Q(\theta,\Omega,\sigma,\beta)\). For a likelihood-based fit, \(Q\) is commonly the negative twice log likelihood, with any constants treated consistently. The workflow below assumes the data, structural model, and likelihood approximation are specified; it does not make a misspecified model correct.
A Practical Sequence of Parameter Blocks
A useful default is to establish structural behavior first, then add flexibility for unexplained between-subject differences, then correlations, and finally covariate effects. The order is a starting point, not an inflexible rule: identifiability, scientific knowledge, and the structure of the data can justify adjustments.
Use scientifically reasonable initial values, and, if useful, a simpler fit to obtain starting values. Estimate the typical structural parameters and basic residual terms while holding random-effects terms at provisional values or using a simpler model. Do not interpret this fit as the final hierarchical result.
Introduce or free the intended random-effects variance components, initially using a diagonal \(\Omega\) when that is a reasonable simplification. Re-estimate the structural parameters and residual terms as well, so they can adapt to the newly estimated variability.
Free selected off-diagonal elements of \(\Omega\), or their correlation parameterization, after the marginal variances have plausible estimates. Re-optimize the structural parameters, variances, and other already active terms. Add correlations only when the data and model can support them.
Introduce covariate coefficients in a controlled sequence, starting from zero or another scientifically justified value. At the final stage, estimate all intended structural, variability, correlation, residual, and covariate parameters together, using the selected final estimation settings.
The key is to distinguish the block being introduced from the parameters being re-optimized. When a new variance is freed, for example, previously estimated typical values should usually be free to adjust. Freezing every earlier block can leave the new block compensating for estimates that are no longer appropriate.
Keep the numerical estimation method and relevant accuracy settings consistent when comparing objective values across stages. If the integration method, quadrature order, or objective-function scaling changes, the objective may change for numerical reasons as well as because the model changed. A change in optimizer controls may be reasonable to obtain a stable fit, but record it and avoid treating the resulting objective difference as a clean comparison of model structure. Earlier material on estimation settings and objective-function changes provides the relevant background.
Why Conditional Stages Can Help
A staged fit can be understood as a succession of conditional optimization problems. If \(A\) is the parameter block currently being estimated and \(B\) is held fixed, a stage finds a value of \(A\) that minimizes \(Q(A,B)\). This can be easier than beginning with every component free, particularly when the starting values of one block would otherwise distort the estimates of another.
But a conditional optimum is not necessarily a joint optimum. The joint problem minimizes \(Q(A,B)\) over both blocks; after \(B\) is freed, the best value of \(A\) may change. Thus, staging is most defensible when each step provides plausible starting values and the final fit releases all intended parameters.
A simple derivation illustrates how structural and variability estimates can be separated temporarily. Suppose independent subject-level estimates \(y_i\) follow \(Y_i\sim N(\theta,\omega^2+s^2)\), where \(\theta\) is a typical value, \(\omega^2\) is between-subject variance, and \(s^2\) is known within-subject estimation variance. With \(\omega^2\) fixed, the likelihood is maximized at the mean of the \(y_i\), because all observations have the same variance. Once \(\theta=\bar{y}\), the log likelihood as a function of total variance \(v=\omega^2+s^2\), ignoring constants, is
Differentiating and setting the derivative to zero gives \(-n/(2v)+\sum_i(y_i-\bar{y})^2/(2v^2)=0\), so the maximizing total variance is \(\sum_i(y_i-\bar{y})^2/n\). Therefore, the unconstrained estimate of \(\omega^2\) is this quantity minus \(s^2\); if it is negative, the nonnegative variance constraint places the maximum at the boundary \(\omega^2=0\). This simplified result motivates a variance stage, but a nonlinear mixed-effects likelihood generally requires numerical optimization and does not reduce to this calculation.
Worked Example: Establishing a Typical Value Before Variability
Consider a hypothetical simplified model with three subject-level estimates, \(y_1=8\), \(y_2=10\), and \(y_3=12\), in concentration units. Suppose the known within-subject variance is \(s^2=1\) concentration-units squared. First hold the between-subject variance at a provisional value of \(1\), and estimate the typical value \(\theta\).
Because all three observations have the same variance, the estimate is their mean:
Next free \(\omega^2\). The squared deviations sum to \((8-10)^2+(10-10)^2+(12-10)^2=4+0+4=8\). The estimated total variance is \(8/3=2.667\) concentration-units squared, so the estimated between-subject variance is
In this equal-variance example, freeing \(\omega^2\) leaves \(\widehat{\theta}=10\), because the mean remains the likelihood estimate of the typical value. In a nonlinear mixed-effects model, structural and variability estimates commonly move together. The staged estimates here are teaching values for a hypothetical model, not results from a clinical study.
Correlations and Covariates: Add Flexibility Deliberately
A diagonal random-effects covariance matrix assumes distinct random effects are uncorrelated. Allowing correlations can represent additional structure, but it also adds parameters and can make estimation more sensitive to the available information. A correlation stage should therefore follow reasonable estimates of the marginal variances, and the full covariance matrix must remain positive definite.
Worked Example: Freeing a Random-Effects Correlation
Suppose a hypothetical two-random-effect model has fitted marginal variances \(0.16\) and \(0.09\), in squared units of their respective random effects. With those variances held fixed for this illustrative profile, the same estimation method gives these objective-function values at candidate correlations:
| Correlation \(\rho\) | Objective value \(Q\) |
|---|---|
| 0.0 | 102.4 |
| 0.3 | 100.8 |
| 0.6 | 100.2 |
| 0.9 | 101.1 |
Among these evaluated candidates, \(\rho=0.6\) has the lowest objective. The corresponding covariance is \(0.6\sqrt{0.16}\sqrt{0.09}=0.6(0.4)(0.3)=0.072\), so the candidate covariance matrix is
Its determinant is \(0.16(0.09)-(0.072)^2=0.0144-0.005184=0.009216>0\), and its diagonal entries are positive, so this two-dimensional matrix is positive definite. The objective is \(2.2\) lower at \(\rho=0.6\) than at \(\rho=0\) in the stated profile. This is a hypothetical illustration of choosing a candidate, not a significance test or evidence that the correlation is scientifically important. In a full stage, other active parameters should generally be re-optimized rather than permanently held at their profile values.
Covariates are often easier to assess after the basic structural and variability components have plausible estimates. Centering a continuous covariate makes the typical parameter interpretable at the centering value; scaling can make numerical optimization easier. Add covariates in a scientifically defensible form, and re-estimate the affected typical parameters and variability terms. A fitted covariate association is not, by itself, evidence of a causal effect.
Worked Example: Adding a Covariate Effect
In a hypothetical clearance model, let \(x=0\) denote the reference group and \(x=1\) the comparison group. Use \(CL_i=\theta_{CL}\exp(\beta x_i+\eta_{CL,i})\), where \(\theta_{CL}\) is clearance for the reference group, \(\beta\) is the covariate coefficient, and \(\eta_{CL,i}\) is an individual random effect. A fit without the covariate gives \(\theta_{CL}=5.0\) L/h and \(Q=250\). After adding the covariate and jointly re-optimizing the active parameters, suppose the hypothetical fit gives \(\theta_{CL}=5.2\) L/h, \(\beta=0.35\), and \(Q=243\).
For the comparison group, the estimated typical clearance is \(5.2\exp(0.35)\) L/h. Since \(\exp(0.35)\approx1.419\), this is approximately \(7.38\) L/h, compared with \(5.2\) L/h for the reference group. The multiplicative difference is about \(41.9\%\), calculated as \((1.419-1)\times100\%\). It is an estimated model association for these hypothetical groups, not a causal conclusion. The objective is 7 units lower in the covariate model under the stated setup; the example does not use that difference to make an inferential decision.
Settings, Diagnostics, and Limits
A staged workflow is useful only if each stage is interpretable. Keep the likelihood approximation and numerical accuracy settings appropriate to the model, and use consistent settings when the purpose is to compare objective values. Revisit convergence tolerances and evaluation limits when necessary, but do not mistake a larger iteration budget or a more permissive tolerance for evidence that the model is well identified.
Practical checks during staging include confirming which parameters are fixed and which are free, checking that the new estimates are finite and within their intended bounds, and verifying that the intended model actually changes at each stage. Record starting values, parameter blocks, estimation settings, and objective values. If a new block produces implausible estimates or the optimizer repeatedly fails to settle, do not simply continue adding blocks: reconsider the starting values, model complexity, or whether the available data support that parameterization. Detailed failure signatures require their own diagnostics.
There are alternatives to this staged sequence. A one-step fit estimates all parameters from the outset and can be efficient when starting values are good and the problem is well conditioned. Simplified models can provide starting values, but their estimates may not transfer directly to the richer model. Staging offers an organized compromise, yet it cannot guarantee the global optimum: a nonlinear objective may have multiple local minima, and the path taken can affect the result. Assessing sensitivity to starting values remains relevant.
Finally, an intermediate conditional optimum is not a final uncertainty estimate for the full model. Once the final parameterization is fitted, evaluate it as a joint model and obtain uncertainty estimates using appropriate settings. If a block cannot be estimated reliably, omitting it or simplifying the model may be more defensible than forcing a staged sequence to produce an answer.
Common Pitfall
Correction: Stage 1 estimates are conditional on the values fixed there. When variability or covariates are added, re-optimize the relevant earlier parameters and finish with a joint fit of all intended free parameters. Otherwise the reported values need not minimize the objective for the stated full model.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- In a staged fit, why should structural parameters usually be re-optimized after random-effects variances are freed?
- A two-random-effect model has variances \(0.25\) and \(0.04\), with correlation \(0.5\). What covariance does this imply?
- A covariate model estimates \(\beta=-0.20\) in an exponential parameter model. What is the multiplicative comparison for a one-unit increase in the covariate?
- Why should objective values be interpreted cautiously if the integration accuracy settings changed between stages?
Solutions. (1) The new variability estimates can change the likelihood’s preferred structural values; keeping them fixed leaves a conditional, not joint, optimum. (2) The covariance is \(0.5\sqrt{0.25}\sqrt{0.04}=0.5(0.5)(0.2)=0.05\), in the product units of the random effects. (3) The multiplier is \(\exp(-0.20)\approx0.819\), or about an 18.1% decrease per one-unit increase under the model. (4) The objective can change because its numerical approximation changed, so the difference cannot be attributed solely to the parameter or model changes.