Tutorials › Nonlinear Mixed-Effects Models (NLME) › Worked Example of Parameter-Uncertainty Propagation

Simulation From Fitted Models · Tutorial 711 of 797

Worked Example of Parameter-Uncertainty Propagation

Learn to carry parameter and interindividual uncertainty through a fitted oral concentration model, and distinguish a typical plug-in profile from population prediction.

Intermediate 12 min read

What You'll Learn

  • Specify the one-compartment oral model and the log-parameter scale used for uncertainty propagation.
  • Use a delta-method gradient and fixed-effect covariance matrix to quantify uncertainty in concentration.
  • Distinguish uncertainty in fitted population parameters from variability between individuals.
  • Interpret plug-in profiles, confidence intervals, and prediction intervals in context.
  • Set up a simulation that carries parameter dependence and interindividual variability through an entire profile.
  • Recognize when a local approximation is inadequate and nonlinear simulation is preferable.

Why Carry Uncertainty Through a Concentration Model?

A fitted concentration-time curve is often displayed as though its estimated parameters were known exactly. They are not: repeated data sets would generally produce different fixed-effect estimates, and individuals also differ through their random effects. Propagating these two sources of uncertainty through the model shows how they affect predicted concentrations, rather than leaving them as uncertainty on a parameter scale.

The distinction matters. A curve evaluated at the fitted typical parameters is a plug-in typical profile; it is not a prediction interval for individual concentrations and it does not show the sampling uncertainty in the fitted parameters. The introductory material on confidence intervals and prediction intervals in simulation distinguishes these targets. Here the focus is a concrete propagation calculation for a one-compartment oral model.

Model, Parameters, and Uncertainty Targets

Assume a single oral dose \(D\), first-order absorption, first-order elimination, no lag time, and a one-compartment disposition model. Let \(CL/F\) and \(V/F\) denote apparent clearance and apparent volume, and let \(k_a\) be the absorption rate constant. For \(t>0\), set \(k=(CL/F)/(V/F)\). The concentration is

$$ C(t)=\frac{D}{V/F}\frac{k_a}{k_a-k} \left(e^{-kt}-e^{-k_a t}\right). $$

The expression has a continuous limiting form when \(k_a=k\), although the displayed quotient should not be evaluated by direct subtraction in that case. The examples below have \(k_a\ne k\). For fixed-effect uncertainty, use the log parameter vector \(\boldsymbol{\theta}=(\log(CL/F),\log(V/F),\log k_a)^\mathsf{T}\), with estimated covariance matrix \(\widehat{\Sigma}\). Log parameters are useful because rate and scale parameters must remain positive.

Definition: Fixed-effect uncertainty is uncertainty in the estimated population parameters, represented here by \(\widehat{\Sigma}\). Interindividual variability is the distribution of individual random effects around those population parameters, represented by \(\boldsymbol{\eta}_i\sim N(\mathbf{0},\Omega)\). These are distinct sources: a parameter draw represents uncertainty about the fitted population, while a random-effect draw represents a different individual within that population.

For a local calculation, approximate \(\boldsymbol{\theta}\) by a multivariate normal distribution centered at its estimate with covariance \(\widehat{\Sigma}\). This approximation and the delta method are most defensible when the estimate is reasonably regular, its covariance is well estimated, and uncertainty is not so large that the nonlinear model changes sharply over the plausible parameter region. The worked propagation assumes no covariance between the fixed-effect estimation error and the new individual’s random effect. If that independence is not appropriate for the intended prediction problem, the joint dependence must instead be represented.

Delta-Method Propagation for One Concentration

For a scalar model output \(C(t;\boldsymbol{\theta})\), let \(\mathbf{g}(t)\) be its gradient with respect to the log parameters, evaluated at the fitted values. A first-order Taylor expansion gives

$$ C(t;\boldsymbol{\theta})\ \approx\ C(t;\widehat{\boldsymbol{\theta}}) +\mathbf{g}(t)^\mathsf{T}(\boldsymbol{\theta}-\widehat{\boldsymbol{\theta}}), \qquad \operatorname{Var}_{\mathrm{fixed}}\{C(t)\}\ \approx\ \mathbf{g}(t)^\mathsf{T}\widehat{\Sigma}\mathbf{g}(t). $$

The quadratic form uses the full covariance matrix, not just its diagonal. Its off-diagonal terms matter whenever parameter estimates are correlated. This is a local approximation to uncertainty in the model-predicted concentration due to fixed-effect estimation; it is not, by itself, an interval for a new individual’s concentration.

Worked Example: Fixed-Effect Uncertainty at Two Hours

Suppose a hypothetical fitted model has \(D=100\) mg, \(CL/F=5.00\) L/h, \(V/F=50.0\) L, and \(k_a=1.00\) h\(^{-1}\). These are illustrative values, not results from a real study. Thus \(k=5/50=0.100\) h\(^{-1}\). At \(t=2.00\) h,

$$ C(2)=\frac{100}{50}\frac{1}{1-0.1} \left(e^{-0.1(2)}-e^{-1(2)}\right) =2.2222(0.81873-0.13534) \approx 1.5187\ \mathrm{mg/L}. $$

Suppose the estimated covariance matrix of the three log parameters, in the stated order, is

$$ \widehat{\Sigma}= \begin{pmatrix} 0.01000&0.00480&-0.00100\\ 0.00480&0.00464&-0.00020\\ -0.00100&-0.00020&0.02000 \end{pmatrix}. $$

Differentiating the concentration expression with respect to the log parameters and evaluating at these fitted values gives the local gradient \(\mathbf{g}(2)\approx(-0.1951,-1.3235,0.4328)^\mathsf{T}\) mg/L per unit change in the corresponding log parameter. The covariance propagation is the quadratic form. Written as diagonal and cross-product contributions, its arithmetic is

$$ \begin{aligned} \mathbf{g}^\mathsf{T}\widehat{\Sigma}\mathbf{g} &\approx (0.1951)^2(0.01000)+(1.3235)^2(0.00464)+(0.4328)^2(0.02000)\\ &\quad+2(-0.1951)(-1.3235)(0.00480) +2(-0.1951)(0.4328)(-0.00100)\\ &\quad+2(-1.3235)(0.4328)(-0.00020)\\ &\approx 0.01513\ (\mathrm{mg/L})^2. \end{aligned} $$

The approximate standard error due to fixed-effect uncertainty is therefore \(\sqrt{0.01513}=0.1230\) mg/L. A pointwise, normal-approximation 95% interval for the model-predicted concentration, holding the random effect at zero, is \(1.5187\pm1.96(0.1230)=[1.278,1.760]\) mg/L. It describes uncertainty in this fitted-model prediction under the approximation; it does not mean there is a 95% probability that a fixed true concentration lies in this realized interval.

Adding Interindividual Variability

For an individual, a common log-normal parameterization is \(\boldsymbol{\theta}_i=\boldsymbol{\theta}+\boldsymbol{\eta}_i\), with \(\boldsymbol{\eta}_i\sim N(\mathbf{0},\Omega)\). For example, the individual apparent clearance is \(CL_i/F=\exp(\log(CL/F)+\eta_{CL,i})\). Each random effect must be carried into the parameter it represents before recalculating \(k\) and the concentration. A population distribution of individual profiles is not obtained by adding random effects to concentrations after evaluating the typical profile.

When fixed-effect uncertainty and individual random effects are independent, a first-order approximation to the variance at a fixed time separates their contributions:

$$ \operatorname{Var}\{C_i(t)\}\ \approx\ \mathbf{g}(t)^\mathsf{T}\widehat{\Sigma}\mathbf{g}(t) +\mathbf{g}(t)^\mathsf{T}\Omega\mathbf{g}(t). $$

This local decomposition is useful for understanding sources of spread. It does not capture nonlinear shifts in the mean or skewness in the concentration distribution. In particular, a nonlinear function of normally distributed log parameters need not itself be normally distributed.

Worked Example: Separating Parameter Uncertainty From Individual Variability

Continue the hypothetical example at 2.00 h. Suppose, for illustration, the random-effect covariance is diagonal, \(\Omega=\operatorname{diag}(0.0400,0.0225,0.0400)\), on the same log-parameter scale. With the gradient from the first example, the local interindividual contribution is

$$ \begin{aligned} \mathbf{g}^\mathsf{T}\Omega\mathbf{g} &\approx (0.1951)^2(0.0400)+(1.3235)^2(0.0225)+(0.4328)^2(0.0400)\\ &\approx 0.04843\ (\mathrm{mg/L})^2, \end{aligned} $$

so its local standard deviation is \(\sqrt{0.04843}=0.2201\) mg/L. The fixed-effect contribution from the first example is \(0.01513\ (\mathrm{mg/L})^2\). Under the stated independence assumption, the total local variance is \(0.01513+0.04843=0.06356\ (\mathrm{mg/L})^2\), giving standard deviation \(\sqrt{0.06356}=0.2521\) mg/L.

This calculation does not say that the 0.2521 mg/L value is an exact prediction standard deviation. It is a first-order approximation around the typical parameter values and zero random effect. A full nonlinear simulation evaluates the concentration after drawing both parameter uncertainty and random effects. Its distribution can be asymmetric, and its center can differ from the typical plug-in concentration of 1.5187 mg/L.

Propagating Uncertainty Across a Profile

For a profile, use the same parameter draw at every time point. Otherwise the simulated curve can jump from one parameter combination to another over time and no longer represents a coherent individual or fitted population. A full simulation separates two useful targets:

  • Uncertainty in a population profile: for each fixed-effect draw, average the concentrations over random effects; compare these population-average profiles across fixed-effect draws.
  • Prediction for individuals: draw fixed effects and individual random effects, then calculate each individual’s complete concentration profile. The resulting distribution includes both sources of uncertainty.

The first target concerns uncertainty about a population quantity; the second concerns variation in individual outcomes as well as uncertainty in the fitted parameters. They should not be given the same label or interpreted interchangeably.

Worked Example: Plug-In Profile Versus Full Uncertainty

For the fitted values in the first example, the plug-in typical profile at three times is obtained by substituting each time into the same model:

$$ \begin{aligned} C(0.5)&=2.2222(e^{-0.05}-e^{-0.5})\approx0.7660\ \mathrm{mg/L},\\ C(2.0)&=2.2222(e^{-0.2}-e^{-2})\approx1.5187\ \mathrm{mg/L},\\ C(8.0)&=2.2222(e^{-0.8}-e^{-8})\approx0.9978\ \mathrm{mg/L}. \end{aligned} $$

These three numbers are a single typical curve, not an uncertainty band. At 2.00 h, the fixed-effect-only local interval calculated earlier is approximately [1.278, 1.760] mg/L. Adding the random-effect contribution gives a local total standard deviation of 0.2521 mg/L; a purely normal local approximation would give \(1.5187\pm1.96(0.2521)=[1.025,2.013]\) mg/L. The latter is only an illustrative local approximation to individual predictive spread. It is not a substitute for simulating the nonlinear model when an actual prediction distribution is required.

A simulation carries the quantities through the model as follows: draw one fixed-effect vector from the fitted uncertainty distribution, draw individual random effects from the fitted random-effect distribution, transform to positive individual parameters, and evaluate the concentration formula at every required time. Across draws, summaries for the population average answer the population-profile question; quantiles across individual profiles answer a predictive question. The fixed-effect draw must be shared across the individuals in a simulated population replicate, since it represents one uncertain fitted model, not a new model for every individual.

Assumptions, Checks, and Alternatives

The delta method assumes a sufficiently smooth model near the fitted values and uses a local linear approximation. It can misrepresent uncertainty when parameter uncertainty is large, parameters are near constraints, the model is strongly curved over plausible draws, or the concentration is near a region where small parameter changes produce large changes in the profile. The multivariate-normal approximation for fixed-effect estimates should also be considered in light of the fit and the parameter covariance estimate; a nearly singular covariance matrix or unstable parameter estimates are warning signs.

For simulation, check that the fixed-effect covariance is symmetric and positive semidefinite, that simulated parameter values are valid for the model, and that the same covariance and parameter ordering are used consistently. Examine simulated profiles for non-finite values, implausible shapes, or extreme tails. If failed draws occur, investigate their cause rather than silently discarding them: failures may indicate a model-domain problem or a poorly behaved uncertainty approximation. Also check that the uncertainty summary matches the target: a population-average uncertainty band should not be constructed from individual-level spread alone.

ApproachWhat it givesMain limitation
Plug-in profileConcentration at fitted typical parametersDoes not quantify parameter uncertainty or individual variability
Delta methodFast local approximation using a gradient and covarianceMay be inaccurate for curved or strongly nonlinear propagation
Simulation from fitted modelNonlinear uncertainty or predictive distributions, depending on what is drawnRequires a correctly specified draw structure and enough valid simulations

Simulation is not automatically superior: it inherits the fitted model and covariance approximation, and its output can be misleading if uncertainty sources are combined twice. In particular, do not treat empirical Bayes estimates as population random-effect draws and then also add a fresh draw from \(\Omega\) for the same individual. Use the distribution appropriate to the target, and include each source of uncertainty once.

Key takeaway: A plug-in curve is one model evaluation. Fixed-effect uncertainty is propagated through the gradient-covariance quadratic form or by drawing fitted parameter vectors; interindividual variability is propagated by drawing random effects and transforming individual parameters. Keep the targets distinct and carry each draw through the complete nonlinear profile.

Check Your Understanding

Answer each question, then compare with the short solutions.

  1. What does the interval [1.278, 1.760] mg/L from the first example quantify?
  2. Why must off-diagonal entries of the fixed-effect covariance matrix be included in the delta-method variance?
  3. In a simulated population replicate, why should one fixed-effect draw be used across the replicate’s individuals?
  4. What is the difference between a population-profile uncertainty band and an individual prediction interval?

Solutions. (1) It is a pointwise, approximate 95% interval for the model-predicted typical concentration at 2 h due to fixed-effect estimation uncertainty, under the local normal approximation. It is not an individual prediction interval. (2) Correlated parameter estimates contribute cross-products to \(\mathbf{g}^\mathsf{T}\widehat{\Sigma}\mathbf{g}\); omitting them changes the propagated variance. (3) The draw represents one uncertain fitted population model shared by the replicate, rather than a different fitted model for each individual. (4) A population-profile band summarizes uncertainty in a population quantity across fitted-parameter draws; an individual prediction interval includes the distribution of individual random effects, and may also include fixed-effect uncertainty if that is the specified target.