Why Nest the Simulation?
A simulated population prediction can vary for two fundamentally different reasons. The population model parameters are uncertain because they were estimated from finite data; separately, individuals differ because their random effects vary around the population model. A simulation that changes both sources at once produces a spread of outputs, but that spread is difficult to interpret unless the two sources are deliberately organized and summarized.
The practical design is nested. First draw a plausible set of population parameters from the uncertainty distribution. Then, conditional on that draw, simulate individual random effects and calculate the target output for each individual. Repeating the inner simulation for each outer parameter draw creates groups of outputs: variation within a group reflects individual random effects, while changes among the groups' conditional summaries reflect parameter uncertainty.
This tutorial holds the scenario and covariates fixed unless otherwise stated. It concerns outputs generated by the model without residual observation error; add residual error only when the simulation target calls for observations rather than the underlying model output. The parameter draws must represent the intended uncertainty distribution, such as valid bootstrap or posterior draws, and the individual effects must be drawn from the model's specified random-effects distribution conditional on each parameter draw.
For example, in a clearance model, the typical clearance may be a component of \(\theta_m\), while an individual's deviation from typical clearance is represented by \(\eta_{mi}\). If clearance is modelled as \(CL_{mi}=CL_{\mathrm{typ},m}\exp(\eta_{mi})\), parameter draws change the typical clearance and may also change the random-effects variance. Inner draws then represent individual clearance deviations conditional on those population parameters. This is a simulation design, not a claim that the two sources have the same scientific meaning.
Separating the Two Components of Spread
For fixed population parameters \(\theta\), define the conditional mean and conditional variance of the output over individual random effects as
The overall output mean averages the conditional means over the parameter-uncertainty distribution. The total variance has two components: the average conditional variance, which is individual variability, and the variance of the conditional means, which is parameter uncertainty in the predicted mean.
The identity follows by adding and subtracting the conditional mean inside the squared deviation. Since \(Y-E(Y)=[Y-\mu(\theta)]+[\mu(\theta)-E(Y)]\), squaring and taking expectations gives a within term, a between term, and a cross-product. The cross-product is zero because \(E[Y-\mu(\theta)\mid\theta]=0\). Thus the remaining terms are \(E_{\theta}[\operatorname{Var}(Y\mid\theta)]\) and \(\operatorname{Var}_{\theta}[\mu(\theta)]\), as stated.
With \(M\) equally weighted parameter draws, let \(\mu_m\) and \(v_m\) denote the conditional mean and variance at draw \(m\), and let \(\bar{\mu}=M^{-1}\sum_m\mu_m\). The variance under the empirical distribution of those draws is calculated as follows. Use the divisor \(M\) here because the goal is to describe that empirical mixture; estimating a variance for a broader distribution from a random sample of draws may require a different estimator.
This decomposition is about the chosen output and target. Parameter uncertainty in an individual prediction can affect both its conditional mean and its conditional variance. If the random-effects variance is itself uncertain across parameter draws, that change contributes to the average within-draw variance as well as potentially changing the conditional means. The two terms still have the same definitions: average conditional variance and variance of conditional means.
Worked Example: Decomposing a Small Nested Simulation
Suppose three equally weighted parameter draws each produce three hypothetical individual clearances, in L/h. The fixed covariates, dose scenario, and model output definition are the same for all draws.
| Parameter draw | Individual outputs (L/h) | Conditional mean (L/h) | Conditional variance \((\mathrm{L/h})^2\) |
|---|---|---|---|
| 1 | 9, 10, 11 | 10 | 2/3 |
| 2 | 13, 14, 15 | 14 | 2/3 |
| 3 | 11, 12, 13 | 12 | 2/3 |
The within-draw variance for the first draw, for example, is the average squared deviation from 10: \(\{(9-10)^2+(10-10)^2+(11-10)^2\}/3=2/3\;(\mathrm{L/h})^2\). The other two draws have the same within-draw variance. Therefore \(V_{\mathrm{individual}}=(2/3+2/3+2/3)/3=2/3\;(\mathrm{L/h})^2\).
The mean of the three conditional means is \(\bar{\mu}=(10+14+12)/3=12\;\mathrm{L/h}\). Their variance under the empirical parameter-draw distribution is \(V_{\mathrm{parameter}}=\{(10-12)^2+(14-12)^2+(12-12)^2\}/3=8/3\;(\mathrm{L/h})^2\). Thus \(V_{\mathrm{total}}=2/3+8/3=10/3=3.333\;(\mathrm{L/h})^2\), and the corresponding standard deviation is \(\sqrt{10/3}=1.826\;\mathrm{L/h}\). The larger component here is parameter uncertainty in the conditional mean, not individual variability.
As a check, the mean of all nine outputs is 12 L/h. Their squared deviations from 12 sum to 30, so their variance using divisor 9 is \(30/9=3.333\;(\mathrm{L/h})^2\), matching the decomposition.
Lognormal Random Effects: A Clearance Calculation
The distinction is especially useful when the model output is nonlinear in a random effect. Suppose a hypothetical clearance model is \(CL=CL_{\mathrm{typ}}\exp(\eta)\), where \(\eta\sim N(0,\omega^2)\). Conditional on a fixed \(CL_{\mathrm{typ}}\) and \(\omega^2\), the expected clearance and its variance are
These expressions use the lognormal moments implied by a normally distributed \(\eta\). In particular, \(CL_{\mathrm{typ}}\) is the median of this lognormal clearance distribution, not its arithmetic mean.
Worked Example: Separating Typical-Clearance Uncertainty From Individual Effects
Consider three equally weighted hypothetical parameter draws with \(CL_{\mathrm{typ}}\) values of 8, 10, and 12 L/h. Let \(\omega^2=0.04\) for all three. Assume the lognormal random-effects model above, with no residual error and the same scenario for each draw.
For each parameter draw, the conditional mean is \(CL_{\mathrm{typ}}\exp(0.02)\). The average conditional mean is therefore \(10\exp(0.02)=10.202\;\mathrm{L/h}\), since the average typical clearance is \((8+10+12)/3=10\;\mathrm{L/h}\).
The parameter component is the variance of the conditional means. Their common multiplier is \(\exp(0.02)\), so its square is \(\exp(0.04)\). The variance of 8, 10, and 12 using divisor 3 is \(8/3\;(\mathrm{L/h})^2\). Keeping full precision, \(V_{\mathrm{parameter}}=(8/3)\exp(0.04)\approx2.7755\;(\mathrm{L/h})^2\).
For the individual component, average the conditional variance across the three parameter draws. The mean squared typical clearance is \((8^2+10^2+12^2)/3=308/3\;(\mathrm{L/h})^2\). Substitution gives \(V_{\mathrm{individual}}=(308/3)\exp(0.04)[\exp(0.04)-1]\approx4.3609\;(\mathrm{L/h})^2\). The total variance is \(2.7755+4.3609=7.1364\;(\mathrm{L/h})^2\), with standard deviation \(\sqrt{7.1364}\approx2.671\;\mathrm{L/h}\).
The result describes the spread in clearance across the empirical mixture of parameter draws and individual effects. It does not mean there is a 95% probability that a fixed population parameter lies within any interval, nor does it isolate uncertainty in the typical clearance alone if other parameters also varied. Here \(\omega^2\) was held fixed specifically to make the two sources transparent.
Weighted Draws and a Practical Nested Workflow
Not every simulation uses equally weighted parameter draws. If draw \(m\) has weight \(w_m\), where \(w_m\geq0\) and \(\sum_m w_m=1\), then \(\bar{\mu}=\sum_m w_m\mu_m\), \(V_{\mathrm{parameter}}=\sum_m w_m(\mu_m-\bar{\mu})^2\), and \(V_{\mathrm{individual}}=\sum_m w_m v_m\). The weights must reflect the intended uncertainty distribution. Do not treat an arbitrary collection of parameter sets as equally representative unless that is the simulation design.
Worked Example: Combining Unequally Weighted Draw Summaries
Suppose an analysis has three hypothetical parameter draws for a target AUC, measured in mg·h/L. Their weights are 0.2, 0.3, and 0.5. Inner simulations have estimated conditional means of 10, 12, and 16 mg·h/L and conditional variances of 1, 4, and 9 \((\mathrm{mg\cdot h/L})^2\), respectively.
The weighted overall mean is \(0.2(10)+0.3(12)+0.5(16)=13.6\;\mathrm{mg\cdot h/L}\). The parameter component is \(0.2(10-13.6)^2+0.3(12-13.6)^2+0.5(16-13.6)^2=2.592+0.768+2.880=6.240\;(\mathrm{mg\cdot h/L})^2\).
The weighted individual component is \(0.2(1)+0.3(4)+0.5(9)=0.2+1.2+4.5=5.900\;(\mathrm{mg\cdot h/L})^2\). The total variance is \(6.240+5.900=12.140\;(\mathrm{mg\cdot h/L})^2\), giving a standard deviation of \(\sqrt{12.140}=3.484\;\mathrm{mg\cdot h/L}\). Reporting only the pooled standard deviation would conceal that 6.240 of the variance is between parameter-draw means and 5.900 is within-draw individual variability.
Specify the model output, covariates, dosing conditions, and whether residual observation error belongs in the target.
Use draws from the intended uncertainty distribution and retain only sets that meet the model's valid-parameter requirements.
For each \(\theta_m\), sample the individual random effects from the specified distribution conditional on that parameter set, then calculate \(Y_{mi}\).
Calculate or estimate \(\mu_m\) and \(v_m\) from the inner outputs, keeping the parameter-draw identifier attached to every result.
Report the average conditional variance and the variance of conditional means separately, then verify that their sum agrees with the total mixture variance.
Diagnostics, Alternatives, and Limitations
A useful diagnostic is to plot or tabulate the conditional means and variances against the outer draw identifier. Large changes in conditional means indicate sensitivity of the prediction mean to parameter uncertainty. Large within-draw variances indicate substantial individual variability under the model. If the estimated components change markedly when the number of inner individuals increases, the inner Monte Carlo error is still influential; if the outer component changes markedly as parameter draws are added, the outer uncertainty distribution is not yet represented stably.
The key alternative is a single pooled simulation in which parameter draws and random effects are sampled together and only one overall spread is reported. That approach can estimate a marginal output distribution when the draws are correctly generated, but it does not by itself explain how much spread comes from each source. A nested design retains the grouping needed for that interpretation. Conversely, if the target is only the overall marginal distribution and no source attribution is needed, pooling may be sufficient, provided the simulation preserves the intended joint sampling structure.
The decomposition depends on the target and assumptions. It attributes variation only to the two sources represented in the construction; changing covariates, residual error, dosing histories, or other random components would add further sources. Finite inner simulations make \(\mu_m\) and \(v_m\) estimates rather than exact conditional moments, so their Monte Carlo error can affect both reported components. A simulated individual count does not repair a misspecified random-effects distribution, and valid outer draws do not compensate for an incorrect individual-level model.
- Pooling away the draw identifier. If the outer parameter index is discarded before calculating conditional summaries, the two components cannot be recovered from the pooled variance alone.
- Calling all spread “parameter uncertainty.” Within-draw random-effects variation is individual variability, even though it appears in the same pooled output distribution.
- Using typical values as arithmetic means. In a lognormal model, the typical value is not generally the arithmetic mean; use the conditional distribution's moments or simulate the effects.
- Ignoring weighting or unequal inner sizes. Outer weights must represent the target distribution. When inner sample sizes differ, calculate conditional summaries first and then apply the intended outer weights; do not let a larger inner sample accidentally give its parameter draw more target weight.
- Confusing Monte Carlo noise with model variation. More simulated individuals can reduce numerical noise in conditional summaries, but cannot reduce the model-implied individual variability itself.
A deliberate design may use matched random-effect draws across parameter sets to reduce Monte Carlo noise in comparisons. If so, document that choice and preserve the nested identifiers. Reusing one individual effect for every parameter draw without a deliberate pairing rationale can make the simulation difficult to interpret and may fail to represent the intended conditional individual distributions.
Check Your Understanding
Answer each question, then compare with the solutions below.
- In a nested simulation, what does variation among individuals within one fixed parameter draw represent?
- Three equally weighted parameter draws have conditional means 5, 7, and 9 mg/L and conditional variances 1, 2, and 3 \((\mathrm{mg/L})^2\). Calculate the parameter component, the individual component, and total variance.
- Why is the median clearance, rather than the arithmetic mean clearance, represented by \(CL_{\mathrm{typ}}\) in \(CL=CL_{\mathrm{typ}}\exp(\eta)\) when \(\eta\) is normal with mean zero?
- A pooled simulation reports only one output variance. Can that number alone identify the parameter and individual components? Explain.
- If outer parameter draws have unequal weights, should each conditional mean automatically receive equal weight? State the correct rule.
Solutions. (1) It represents individual random-effects variability conditional on that fixed population-parameter draw. (2) The mean of the conditional means is \((5+7+9)/3=7\;\mathrm{mg/L}\). The parameter component is \(\{(5-7)^2+(7-7)^2+(9-7)^2\}/3=8/3=2.667\;(\mathrm{mg/L})^2\). The individual component is \((1+2+3)/3=2.000\;(\mathrm{mg/L})^2\). Total variance is \(2.667+2.000=4.667\;(\mathrm{mg/L})^2\). (3) Since \(\exp(\eta)\) has median 1 when \(\eta\) is normal with mean zero, the median of \(CL\) is \(CL_{\mathrm{typ}}\); its arithmetic mean is \(CL_{\mathrm{typ}}\exp(\omega^2/2)\). (4) No. The pooled variance combines both components, so conditional summaries or equivalent retained grouping information are needed to separate them. (5) No. Use the weights that represent the intended parameter-uncertainty distribution, normalized to sum to one.