Tutorials › Nonlinear Mixed-Effects Models (NLME) › Parameter Uncertainty From Bayesian Posterior Draws

Simulation From Fitted Models · Tutorial 689 of 797

Parameter Uncertainty From Bayesian Posterior Draws

Use each joint posterior draw of an NLME parameter vector to calculate outcomes, retaining dependence among fixed effects and variance components.

Intermediate 10 min read

What You'll Learn

  • Define posterior-draw propagation for an NLME outcome as a transformation of joint parameter draws.
  • Explain why fixed effects and variance components must remain paired within each posterior draw.
  • Calculate model-based outcomes and summaries from a small set of posterior draws.
  • Interpret posterior intervals for simulated outcomes without confusing them with frequentist confidence intervals.
  • Diagnose common problems involving parameter constraints, MCMC mixing, and independent resampling.
  • Distinguish parameter uncertainty from additional random-effect and residual variability.

From Posterior Parameters to Outcomes

An NLME model can produce different predictions not only because subjects vary, but also because the population parameters are uncertain. Bayesian fitting represents this latter uncertainty with draws from a joint posterior distribution. To propagate it, evaluate the outcome of interest using each draw, then summarize the resulting set of outcomes. This is useful when the outcome is a nonlinear function of several parameters, or when the posterior distribution is asymmetric or constrained.

The central practical rule is to keep each posterior draw intact. A draw containing fixed effects and variance components is one joint parameter set, not a collection of unrelated marginal values. Pairing a fixed-effect draw with a variance-component draw from a different iteration can change the distribution of predicted outcomes.

Definition: Let \(D\) denote the observed data and let \(\boldsymbol{\theta}=(\boldsymbol{\beta},\boldsymbol{\Omega},\sigma^2,\ldots)\) denote the NLME parameters, including fixed effects \(\boldsymbol{\beta}\), random-effect covariance parameters \(\boldsymbol{\Omega}\), and, when applicable, residual variance \(\sigma^2\). Posterior-draw propagation uses draws \(\boldsymbol{\theta}^{(m)}\) from the joint posterior distribution \(p(\boldsymbol{\theta}\mid D)\) to calculate an outcome \(Q^{(m)}=g(\boldsymbol{\theta}^{(m)})\), or a model-based outcome conditional on that draw. The collection of \(Q^{(m)}\) approximates the posterior distribution of the outcome.

The function \(g\) might return a typical-subject concentration, a population-level probability, or another quantity defined by the model. The calculation does not require that \(g\) be linear. For a nonlinear model, in general, \(g\) evaluated at the posterior mean parameters is not the same as the posterior mean of \(g(\boldsymbol{\theta})\). Thus, a plug-in prediction can miss both the centre and the spread of the outcome distribution.

What the Posterior Draws Represent

Suppose the fitted model supplies \(M\) joint posterior draws. For every iteration \(m\), use all relevant parameter values from that same draw to calculate the outcome. For example, \(\boldsymbol{\beta}^{(m)}\) must stay paired with \(\boldsymbol{\Omega}^{(m)}\). This preserves posterior dependence, including dependence between fixed effects and variance components. It also preserves relationships among multiple fixed effects, such as clearance and volume.

$$ Q^{(m)} = g\!\left(\boldsymbol{\theta}^{(m)}\right), \qquad m=1,\ldots,M $$

The empirical distribution of the \(Q^{(m)}\) values is a Monte Carlo approximation to the posterior distribution of \(Q\). For instance, the posterior mean is estimated by \(\bar{Q}=M^{-1}\sum_{m=1}^M Q^{(m)}\), and posterior quantiles can be estimated from the ordered outcome draws. A central 95% posterior credible interval is commonly estimated by the 2.5th and 97.5th percentiles. Conditional on the model and data, that interval contains 95% of the posterior probability for \(Q\); it is not a statement that a fixed unknown value has a 95% probability of lying in the interval under a frequentist interpretation.

The justification is direct. For any threshold \(q\), the posterior probability of an outcome no greater than \(q\) is

$$ P(Q\leq q\mid D) = \int \mathbb{I}\!\left\{g(\boldsymbol{\theta})\leq q\right\} p(\boldsymbol{\theta}\mid D)\,d\boldsymbol{\theta} \approx \frac{1}{M}\sum_{m=1}^M \mathbb{I}\!\left\{g(\boldsymbol{\theta}^{(m)})\leq q\right\}. $$

Here \(\mathbb{I}\{\cdot\}\) equals 1 when its condition is true and 0 otherwise. The approximation follows by replacing the posterior integral with an average over posterior draws. As the number of effectively independent draws increases, this average estimates the posterior probability more precisely. Nonlinear transformations are handled by applying \(g\) to every draw, rather than by assuming the transformed outcome has a convenient distribution.

Worked Example: Typical Concentration From Joint Fixed-Effect Draws

Consider hypothetical posterior draws for a one-compartment, intravenous bolus model. For a typical subject, the concentration at time \(t\) is \(C(t)=\frac{\text{Dose}}{V}\exp(-CLt/V)\). Let the hypothetical dose be 100 mg and \(t=5\) h. The paired posterior draws below give clearance \(CL\) in L/h and volume \(V\) in L.

Draw \(m\)\(CL^{(m)}\) (L/h)\(V^{(m)}\) (L)
1420
2525
3620
4840

For draw 1, \(C^{(1)}(5)=\frac{100}{20}\exp[-(4/20)(5)]=5\exp(-1)=1.839\) mg/L. For draw 2, \(C^{(2)}(5)=\frac{100}{25}\exp[-(5/25)(5)]=4\exp(-1)=1.472\) mg/L. For draw 3, \(C^{(3)}(5)=\frac{100}{20}\exp[-(6/20)(5)]=5\exp(-1.5)=1.116\) mg/L. For draw 4, \(C^{(4)}(5)=\frac{100}{40}\exp[-(8/40)(5)]=2.5\exp(-1)=0.9197\) mg/L.

The mean of these four illustrative outcome draws is \((1.839+1.472+1.116+0.9197)/4=1.337\) mg/L. This is an estimate of the posterior mean concentration based on only four draws, not a reliable real-world posterior summary. It shows the calculation: each prediction uses the paired \(CL\) and \(V\) from one draw. The predicted concentration is for a typical subject; no subject-level random effect or residual error was sampled here.

Variance Components Also Affect Outcomes

A variance component can influence an outcome even when the target is a population probability or fraction, rather than an individual simulated observation. Consider a lognormal random-effects model for clearance:

$$ \log(CL_i)=\log(CL_{\mathrm{typ}})+\eta_i,\qquad \eta_i\sim N(0,\omega^2). $$

Here \(CL_{\mathrm{typ}}\) is typical clearance, and \(\omega^2\) is the variance of the subject random effect on the log-clearance scale. For a fixed threshold \(c>0\), the model-implied fraction of subjects with \(CL_i<c\) is \(P(CL_i<c\mid CL_{\mathrm{typ}},\omega^2)=\Phi((\log c-\log CL_{\mathrm{typ}})/\omega)\), where \(\omega=\sqrt{\omega^2}\). This expression uses both the fixed effect and the variance component. If either parameter is uncertain, the fraction is uncertain too.

Worked Example: Propagating Uncertainty Into a Population Fraction

Suppose four hypothetical joint posterior draws for \(CL_{\mathrm{typ}}\) and \(\omega^2\) are \((4,0.04)\), \((5,0.09)\), \((6,0.16)\), and \((8,0.25)\), with clearance in L/h. Calculate the model-implied fraction below \(c=6\) L/h for each draw. The standard deviations on the log scale are \(\omega=0.2, 0.3, 0.4,\) and \(0.5\), respectively.

For draw 1, the standardized threshold is \((\log 6-\log 4)/0.2=2.027\), giving \(\Phi(2.027)\approx0.9787\). For draw 2, it is \((\log 6-\log 5)/0.3=0.6076\), giving approximately 0.7283. For draw 3, it is \((\log 6-\log 6)/0.4=0\), giving 0.5000. For draw 4, it is \((\log 6-\log 8)/0.5=-0.5754\), giving approximately 0.2825.

The mean of these four model-implied fractions is \((0.9787+0.7283+0.5000+0.2825)/4=0.6224\). The interpretation is an estimated posterior mean of the population fraction below 6 L/h, conditional on the model and data, based on these four hypothetical parameter draws. The calculation includes uncertainty in both the typical clearance and the random-effect variance. It does not mean that 62.24% of a particular simulated sample must fall below the threshold.

Why Pairing Matters

The joint posterior can associate large values of one parameter with large or small values of another. Preserving each draw retains this dependence in the outcome distribution. Independently sampling from separate marginal distributions discards that information and generally creates parameter combinations that are not representative of the joint posterior. The resulting mean might happen to be similar for a particular outcome, but the spread, tail probabilities, and quantiles can change.

Worked Example: Independent Resampling Changes a Predicted Fraction

Use the four paired draws from the clearance-threshold example. Their mean fraction is 0.6224. Now keep the same marginal values but, purely to illustrate the consequence of breaking pairing, assign the variances in the order \(0.25, 0.04, 0.09, 0.16\) to typical clearances \(4, 5, 6, 8\) L/h. The marginal values have not changed, but the joint pairings have.

The recalculated fractions are \(\Phi(\log(6/4)/0.5)=\Phi(0.8109)\approx0.7913\), \(\Phi(\log(6/5)/0.2)=\Phi(0.9116)\approx0.8189\), \(\Phi(0)=0.5000\), and \(\Phi(\log(6/8)/0.4)=\Phi(-0.7192)\approx0.2360\). Their mean is \((0.7913+0.8189+0.5000+0.2360)/4=0.5866\), rather than 0.6224.

This rearrangement is not a valid substitute for the original posterior draws: it changes their joint distribution. It demonstrates why separately drawing parameter columns can alter a nonlinear outcome. In actual computation, use the full parameter vector from each posterior iteration, not a hand-chosen rearrangement.

Implementation, Checks, and Limitations

A reliable workflow is simple, but it depends on having usable posterior draws and applying the model consistently.

1
Identify the outcome.
Write down the function or model-based quantity to calculate, its units, and which parameters it requires.
2
Take one joint draw.
Read all required fixed effects and variance components from the same posterior iteration. If parameters were stored on transformed scales, convert them draw by draw to the scale required by the model.
3
Calculate and validate the outcome.
Evaluate the model using that draw. Check units, parameter constraints, and whether the resulting value is plausible under the stated model.
4
Summarize the outcome draws.
Calculate the requested posterior mean, quantiles, or probability. Report the outcome scale and the uncertainty summary used.

Before interpreting the summaries, check that the posterior sampling has mixed adequately and that the draws cover the target posterior. Trace plots, multiple-chain comparisons, and effective sample size help reveal poor exploration or strong serial autocorrelation. A large raw number of saved iterations does not guarantee a large effective sample size. Also check that variance parameters remain in their allowed range and covariance matrices remain valid; a calculation that produces impossible values may indicate a scale, parameterization, or coding error.

When parameters are represented on transformed scales, apply the inverse transformation to each draw before using a model equation defined on the natural scale. For example, if a log-variance is sampled, exponentiate each draw to obtain the variance. Transforming a posterior mean instead is generally not equivalent to averaging the transformed draws.

This procedure propagates parameter uncertainty, not every source of simulation variability. The typical-subject calculation above contains no random effects or residual error. The threshold example computes a population fraction from the random-effects distribution but does not generate a finite sample of subjects. If a task calls for individual outcomes, additional random effects and possibly residual errors must be generated conditional on parameter draws; those introduce further variability and should be distinguished from the spread caused by the posterior parameters alone.

Posterior-draw propagation is especially useful when the outcome is nonlinear, parameters are constrained, or the joint posterior is not well represented by a multivariate normal approximation. A covariance-matrix approximation can be convenient when posterior uncertainty is nearly symmetric and locally well behaved, but it may poorly represent skewness or boundary behavior. A nonparametric bootstrap instead approximates sampling variability through repeated resampling and refitting; posterior draws represent uncertainty conditional on the specified Bayesian model, prior, and observed data. Neither source of uncertainty should be described as interchangeable with the other.

Common Pitfalls

  • Pairing parameters independently. Wrong reasoning: each marginal draw is valid, so any pairing is valid. Correction: the target is the joint posterior; retain the parameter vector as sampled at each iteration.
  • Predicting only at posterior mean parameters. Wrong reasoning: the mean parameter vector represents the mean prediction. Correction: for nonlinear \(g\), generally \(g(E[\boldsymbol{\theta}\mid D])\ne E[g(\boldsymbol{\theta})\mid D]\); calculate the outcome draw by draw.
  • Calling all outcome spread parameter uncertainty. Wrong reasoning: simulated outcomes vary, therefore the parameters are uncertain. Correction: separate posterior variation in parameters from random effects, residual error, and finite-sample variation introduced by additional simulation.
  • Reporting a credible interval as a confidence interval. Correction: describe it as a posterior interval conditional on the model and data, and state the outcome and units.
Key takeaway: To propagate Bayesian parameter uncertainty, calculate the outcome for every joint posterior draw and summarize those outcome draws. Keeping fixed effects and variance components paired preserves their posterior dependence; independent marginal resampling does not.

Check Your Understanding

Answer each question, then compare with the short solutions below.

  1. What must be kept together when calculating outcomes from posterior draws, and why?
  2. If \(Q^{(m)}=g(\boldsymbol{\theta}^{(m)})\), what does the empirical distribution of the \(Q^{(m)}\) values approximate?
  3. In the lognormal-clearance example, what is the standard deviation on the log scale when \(\omega^2=0.09\)?
  4. Why can the prediction at posterior mean parameters differ from the posterior mean prediction?
  5. Does a typical-subject prediction that uses only fixed effects and variance-component draws include residual error or a newly sampled subject random effect?

Solutions. (1) Keep the entire parameter vector from the same posterior iteration together, because its components have a joint posterior distribution and may be dependent. (2) It approximates the posterior distribution of the model-based outcome, conditional on the model and data. (3) \(\omega=\sqrt{0.09}=0.3\) on the log-clearance scale. (4) Nonlinearity means, in general, \(g(E[\boldsymbol{\theta}\mid D])\ne E[g(\boldsymbol{\theta})\mid D]\). (5) No. Those are additional sources of variability that must be generated separately if individual simulated outcomes are required.