Separating Between-Subject Spread from Parameter Uncertainty

What the Two Sources of Spread Mean

A PBPK simulation can produce a wide distribution of exposures for two distinct reasons. Individuals can differ even when the population parameters are held fixed, for example because their body size or enzyme activity differs. In addition, the population parameters themselves may be uncertain, so predictions can shift when those parameters change. A nested simulation represents these levels separately: an outer draw selects a parameter set, and inner draws simulate individuals under that set.

This distinction matters when a wide prediction distribution is used to inform a dose or to describe exposure in a target population. A broad distribution caused mainly by individual differences has a different interpretation from one that is broad because plausible parameter values produce substantially different population predictions. The total predictive spread alone does not say which explanation dominates.

Here, \(\boldsymbol{\theta}\) is drawn from a distribution representing parameter uncertainty. That distribution might come from a Bayesian analysis, a bootstrap, or another specified uncertainty analysis; its interpretation depends on how it was constructed. Conditional on each draw, subject-level variability is generated from the population model. The individual exposure \(Y\) may be an AUC, a concentration at a specified time, or another defined endpoint.

The Variance Decomposition and Its Meaning

The law of total variance gives the exact decomposition for the predictive distribution obtained by drawing both \(\boldsymbol{\theta}\) and an individual:

The first term is the average conditional subject-level variance across uncertain parameter values. The second is the variance of the conditional mean across parameter values. Both terms have units equal to the square of the exposure units. Their sum is the variance of the predictive distribution that includes both levels of randomness.

A short derivation shows why these terms add. Write \(m(\boldsymbol{\theta})=E[Y\mid\boldsymbol{\theta}]\) and let \(\mu=E[Y]=E_{\boldsymbol{\theta}}[m(\boldsymbol{\theta})]\). Add and subtract \(m(\boldsymbol{\theta})\) from \(Y-\mu\), then expand the square:

The cross term is zero: conditional on \(\boldsymbol{\theta}\), \(m(\boldsymbol{\theta})-\mu\) is fixed and \(E[Y-m(\boldsymbol{\theta})\mid\boldsymbol{\theta}]=0\). The first remaining expectation is \(E_{\boldsymbol{\theta}}[\operatorname{Var}(Y\mid\boldsymbol{\theta})]\); the second is \(\operatorname{Var}_{\boldsymbol{\theta}}(m(\boldsymbol{\theta}))\). This establishes the decomposition.

The first term is not necessarily the variance at one “typical” parameter set. It averages conditional variances over the parameter-uncertainty distribution. Also, parameter uncertainty can change the shape or spread of the conditional distribution, not just its mean. Such changes are reflected in how \(\operatorname{Var}(Y\mid\boldsymbol{\theta})\) varies across parameter sets; the first term averages those conditional variances. The second term specifically measures how conditional means vary.

Estimating the Components from Nested Outputs

Let \(j=1,\ldots,J\) index outer parameter draws and \(i=1,\ldots,n\) index inner subject simulations at each draw. Write the simulated exposure as \(Y_{ij}\). For outer draw \(j\), calculate the inner sample mean \(\bar{Y}_{j}\) and sample variance \(s_j^2\). A direct estimate of the average conditional variance is the average of the inner sample variances:

The variation among the \(\bar{Y}_{j}\) values reflects variation in conditional means, but it also contains Monte Carlo noise because each conditional mean was estimated from only \(n\) individuals. If all inner runs have the same size, the expected contribution of this noise to the sample variance of \(\bar{Y}_{j}\) is approximately \(\widehat{W}/n\). A corresponding estimate of the between-parameter variance is

This correction is an approximation based on independent outer draws, independent inner simulations conditional on each draw, and a common inner sample size. In a finite run, \(\widehat{B}\) can be negative because of sampling noise; that does not imply a negative true variance. It signals that the between-parameter component is not reliably resolved by that simulation. Increasing the number of outer draws and, where needed, inner simulations can improve precision. With unequal inner sample sizes, the noise contribution depends on each \(n_j\), so the equal-size correction should not be applied unchanged.

A simpler descriptive option is to report the outer conditional means and variances separately, together with their Monte Carlo precision. This makes the hierarchy visible without treating a noisy decomposition as exact. The general Monte Carlo diagnostics for means and variances remain important: more simulated rows do not guarantee that the outer uncertainty distribution or its tails are adequately represented.

Why Pooled Outputs Can Mislead

Pooling every \(Y_{ij}\) gives a sample from the combined predictive distribution, provided the outer draws and inner sampling correctly represent the intended distributions. Its spread includes both conditional subject-level variation and variation in conditional means. Therefore, the pooled standard deviation is useful for describing total predictive spread, but it cannot by itself attribute that spread to one source.

The following hypothetical calculation shows why conditional results also matter. Suppose two equally likely parameter sets give normally distributed AUCs, each with a conditional standard deviation of 2 mg·h/L. Their conditional means are 10 and 20 mg·h/L. The average conditional variance is \(4\) (mg·h/L)\(^2\), and the between-parameter variance of means is \(25\) (mg·h/L)\(^2\), since the overall mean is 15 mg·h/L. Total variance is 29 (mg·h/L)\(^2\), giving a pooled standard deviation of \(\sqrt{29}=5.385\) mg·h/L.

At the threshold of 15 mg·h/L, the first conditional distribution has probability \(P(Y>15\mid\boldsymbol{\theta}_1)=1-\Phi(2.5)\approx0.0062\); the second has probability \(P(Y>15\mid\boldsymbol{\theta}_2)=\Phi(2.5)\approx0.9938\), where \(\Phi\) is the standard normal cumulative distribution function. With equal weights, the combined probability is approximately \((0.0062+0.9938)/2=0.5000\). Reporting only the combined probability hides the fact that predictions at the two parameter sets behave very differently.

This illustrates a practical comparison. A fixed-parameter simulation describes subject variability conditional on the chosen parameters, but omits parameter uncertainty. A pooled nested simulation describes the combined prediction, but does not explain its sources. A nested summary reports conditional distributions or moments across parameter draws, allowing both questions to be addressed. None of these summaries alone establishes that the parameter uncertainty distribution is appropriate; that depends on the underlying evidence and modeling choices.

Assumptions, Diagnostics, and Common Pitfalls

Interpret the decomposition as intended when the outer draws represent the stated uncertainty distribution for population parameters, and the inner simulations represent the target population conditional on each draw. The target exposure and time point must be consistent across runs. Inner subjects should be generated according to the specified population model; if simulated individuals are dependent, or if outer and inner draws are coupled in an unaccounted way, the simple variance estimates may not apply. The variance identity itself is general, but estimating its terms requires the simulation to represent this hierarchy.

• Check the hierarchy. Confirm that each inner run uses its associated outer parameter draw and that subject-level variation is regenerated as intended.
• Check outer coverage. Inspect whether the conditional means and variances are stable as additional parameter draws are added. A large inner sample cannot compensate for too few outer draws.
• Check inner precision. If conditional means fluctuate noticeably when inner sample size is increased, their apparent between-parameter spread may be substantially simulation noise.
• Check units and scales. Conditional variances and between-parameter variances must use squared exposure units before they are added.
• Do not infer source from a pooled interval alone. A wide interval is evidence of total predictive spread, not a diagnosis of which component produced it.

A common mistake is to subtract the conditional standard deviation from the pooled standard deviation and call the result “parameter uncertainty.” Standard deviations do not add in this decomposition. Under the stated hierarchy, subtract variances: the between-parameter variance is the total variance minus the average conditional variance. Even that subtraction is an estimate affected by Monte Carlo error, and it specifically concerns variation in conditional means—not every way parameter uncertainty can alter conditional distributions.