Why Predictions Need Two Levels
A simulated outcome can vary for two fundamentally different reasons. The model parameters may be uncertain, so plausible parameter values imply different population predictions. In addition, individuals may differ even when the parameters are held fixed. A predictive distribution that combines these sources is useful, but its spread does not by itself say which source is responsible.
Nested simulation keeps the two levels visible. An outer draw selects a plausible parameter vector; an inner simulation then generates individuals conditional on that draw. The summaries can describe uncertainty in a population-level prediction across outer draws, subject-level spread within a draw, or the full predictive distribution. These are different targets and should not be reported as if they were interchangeable.
The distribution used for \(\Theta\) must match the prediction target. It may, for example, be represented by bootstrap draws or posterior draws as covered in the earlier material on parameter uncertainty. This tutorial treats the outer distribution as specified. The inner distribution must likewise match the intended individual-level simulation. If a prediction includes residual error, that error contributes to its conditional variance; if the target excludes residual error, it should not be added merely to make the simulation look more realistic.
The Variance Decomposition
For a single predictive outcome \(Y\), the total variance can be divided into a between-parameter component and an average within-parameter component. The decomposition requires finite second moments and a coherent joint model in which \(Y\) is generated conditional on \(\Theta\).
A short derivation makes the distinction precise. Write \(\mu=E(Y)=E_{\Theta}\{m(\Theta)\}\), and add and subtract \(m(\Theta)\) inside the centered outcome:
After squaring and taking expectations, the cross-product is zero because \(E\{Y-m(\Theta)\mid\Theta\}=0\). The first squared term has expectation \(E_{\Theta}\{v(\Theta)\}\), while the second has expectation \(\operatorname{Var}_{\Theta}\{m(\Theta)\}\). Therefore,
This is a decomposition of the variance of an individual predictive outcome. It is not the same as the variance of the mean outcome in a future sample of \(n\) individuals. Under a common fixed \(\Theta\) and conditionally independent individuals, the individual component of the variance of their sample mean is \(v(\Theta)/n\); uncertainty in \(m(\Theta)\) remains because the parameter draw is shared. Other sampling designs or dependence assumptions require their own calculation.
Worked Example: Separating Two Sources of Predictive Spread
Suppose a hypothetical outcome \(Y\), measured in mg/L, follows \(Y=\theta+\eta\). The parameter distribution assigns probability \(1/2\) to each of \(\theta=9\) and \(\theta=11\) mg/L. Conditional on either parameter value, individual effects are \(\eta=-1\) and \(\eta=1\) mg/L with equal probability. Assume the effects are generated independently of the parameter draw.
Conditional population means. For either value of \(\theta\), \(E(\eta\mid\theta)=0\), so \(m(\theta)=\theta\). The mean across parameter draws is \((9+11)/2=10\) mg/L. The between-parameter variance of conditional means is \[ \{(9-10)^2+(11-10)^2\}/2=1\;(\mathrm{mg/L})^2. \]
Conditional individual variance. At either parameter value, the two outcomes are one unit below and above \(\theta\). Thus \(v(\theta)=\{(-1)^2+1^2\}/2=1\;(\mathrm{mg/L})^2\), and the average conditional variance is also 1 \((\mathrm{mg/L})^2\).
Total predictive variance. The decomposition gives \(1+1=2\;(\mathrm{mg/L})^2\), so the predictive standard deviation is \(\sqrt{2}=1.414\) mg/L. The first variance component describes uncertainty in the population mean; the second describes individual spread at a fixed parameter value. Neither component alone describes the full individual predictive distribution.
Computing Nested Summaries
In a nested simulation, each outer draw \(m\) has \(R\) inner outcomes \(Y_{mr}\), where \(m=1,\ldots,M\) and \(r=1,\ldots,R\). The inner mean \(\bar{Y}_m\) estimates \(m(\theta^{(m)})\); the inner sample variance \(s_m^2\) estimates \(v(\theta^{(m)})\). If the conditional mean is available analytically, use it rather than adding inner Monte Carlo noise to its estimate.
A basic nested summary reports the distribution of the \(M\) conditional means separately from the average of the \(M\) within-draw variances. The sample variance of the inner means is not automatically a clean estimate of parameter uncertainty: each inner mean has Monte Carlo error. With equal \(R\), independent inner simulation, and finite variances, its expected value conditional on the outer parameter draws is the sample variance of their true conditional means plus the average \(v(\theta^{(m)})/R\). Thus it contains an extra component that decreases as \(R\) increases.
Worked Example: Correcting for Inner-Simulation Noise
Consider two fixed outer parameter values, 9 and 11 mg/L, and the additive individual effect \(\eta\) from the previous example, with variance 1 \((\mathrm{mg/L})^2\). Generate \(R=2\) independent inner outcomes at each outer value. The inner mean for draw \(m\) is \(\bar{Y}_m=\theta^{(m)}+\bar{\eta}_m\), where \(E(\bar{\eta}_m)=0\) and \(\operatorname{Var}(\bar{\eta}_m)=1/2\).
The sample variance of the two fixed conditional means 9 and 11 is \[ S_{\theta}^{2}=\frac{(9-10)^2+(11-10)^2}{2-1}=2\;(\mathrm{mg/L})^2. \] For the two simulated inner means, the sample variance is \(B=(\bar{Y}_1-\bar{Y}_2)^2/2\). The difference in their conditional means is \(-2\) mg/L, and the variance of the difference in their independent inner-mean errors is \(1/2+1/2=1\;(\mathrm{mg/L})^2\). Hence \[ E(B\mid\theta^{(1)}=9,\theta^{(2)}=11)=\frac{(-2)^2+1}{2}=2.5\;(\mathrm{mg/L})^2. \] The additional \(0.5\;(\mathrm{mg/L})^2\) is inner-mean noise, not extra parameter uncertainty. If instead the two inner effects at each outer value are deliberately paired as \(-1\) and \(1\), both inner means equal their conditional means exactly and \(B=2\;(\mathrm{mg/L})^2\). If the outer parameter draws are also newly sampled independently from a distribution with variance 1 \((\mathrm{mg/L})^2\), averaging over those outer draws as well gives an expected \(B=1+0.5=1.5\;(\mathrm{mg/L})^2\). These are different sampling conditions and should not be conflated.
An estimated between-parameter variance can be formed by subtracting the estimated inner-mean noise from the observed variance of inner means. With equal \(R\), the corresponding correction is the average \(s_m^2/R\). This is a Monte Carlo adjustment, not a license to report a negative variance estimate as meaningful parameter uncertainty. A negative adjusted value signals that the between-parameter component is small relative to simulation noise, or that the replicate counts are inadequate; increase the outer and/or inner replication or use analytically computed conditional means where possible.
Nonlinear Predictions and a Lognormal Example
For nonlinear transformations, conditional means and variances need not equal the transformation of a parameter estimate or the variance on the transformed scale. This is one reason to calculate summaries at each outer draw rather than transform only a single point estimate.
Worked Example: A Positive Outcome With Lognormal Individual Variation
In a hypothetical model, let \(Y=\theta\exp(\eta)\), where \(Y\) and \(\theta\) are measured in mg/L and \(\eta\mid\theta\sim N(0,0.2^2)\). Let the outer parameter draw be equally likely to be 8 or 12 mg/L. Assume the stated normal distribution is the complete individual-effect model for this example.
For a normal random variable with mean zero and variance \(0.2^2\), \(E\{\exp(k\eta)\}=\exp(k^2(0.2^2)/2)\). Therefore the conditional mean is \(m(\theta)=\theta\exp(0.02)\). The two conditional means are \(8\exp(0.02)=8.1616\) mg/L and \(12\exp(0.02)=12.2424\) mg/L. Their average is \(10\exp(0.02)=10.2020\) mg/L.
The conditional second moment is \(E(Y^2\mid\theta)=\theta^2\exp(0.08)\), while the squared conditional mean is \(\theta^2\exp(0.04)\). Thus \(v(\theta)=\theta^2\{\exp(0.08)-\exp(0.04)\}\). Averaging across the two parameter values gives \[ E_{\Theta}\{v(\Theta)\} =\frac{8^2+12^2}{2}\{\exp(0.08)-\exp(0.04)\} =104(0.0424763) =4.4175\;(\mathrm{mg/L})^2. \]
The variance of the two conditional means, with equal probabilities, is \(4\exp(0.04)=4.1632\;(\mathrm{mg/L})^2\). The total individual predictive variance is therefore \(4.4175+4.1632=8.5808\;(\mathrm{mg/L})^2\), corresponding to a predictive standard deviation of approximately 2.929 mg/L. The result distinguishes within-parameter individual spread from uncertainty across the two conditional means; it does not imply that every individual has the same variance or that the predictive distribution is normal.
Diagnostics, Alternatives, and Limitations
Nested summaries are most useful when their simulation design and reported target are explicit. Record the number of outer draws and inner replicates, the random-effects and residual-error components included in \(Y\), and whether reported intervals summarize conditional means or individual outcomes. Check that the inner outcomes at each outer draw are generated under that draw’s parameters, and that the intended dependencies among simulated quantities are preserved.
- Check inner precision. Compare inner means across repeated runs or increase \(R\). Large changes in the conditional means indicate that inner simulation noise may obscure the between-parameter pattern.
- Check outer precision. Plot or summarize the sequence of outer-draw means as \(M\) grows. A few extreme parameter draws can dominate variance or tail summaries, particularly with skewed predictions.
- Keep the estimand visible. A spread across conditional means is not an individual-level prediction interval. Conversely, the spread of all pooled individual outcomes includes both sources and cannot identify either component by itself.
- Inspect model plausibility. Implausible parameter draws or unstable conditional simulations can produce extreme predictions. A large predictive spread may indicate genuine uncertainty, a model problem, or both; variance decomposition alone cannot decide which.
A plug-in prediction evaluates the model at one parameter estimate. It is simpler and can be adequate for a conditional point prediction, but it suppresses parameter uncertainty and cannot represent a predictive distribution over uncertain parameters. Pooling every simulated outcome is a direct way to approximate the full predictive distribution, but it combines the two sources of spread. Nested summaries retain that pooled distribution while also separating its variance components. When conditional moments can be integrated analytically or computed accurately by numerical integration, those approaches reduce inner Monte Carlo error; they do not remove uncertainty in the outer parameter distribution.
The variance identity applies to finite-variance outcomes, but it does not by itself decompose quantiles or guarantee that a predictive interval has a particular shape. Quantiles should be calculated from the simulated distribution corresponding to the stated target. In addition, a parameter distribution represents uncertainty according to its source and assumptions; it is not automatically a distribution of biological differences among individuals. Careful labels prevent epistemic uncertainty about model parameters from being mistaken for subject-level heterogeneity.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- In a nested simulation, what does the variance across exact conditional means represent, and what does the average conditional variance represent?
- For fixed outer conditional means 9 and 11 with units mg/L, suppose each inner mean is based on two independent outcomes with conditional variance 1 \((\mathrm{mg/L})^2\). What is the expected sample variance of the two inner means?
- Why can the variance of pooled individual outcomes not be reported as uncertainty in the population mean alone?
- Under what assumption does the individual-variance contribution to the variance of a sample mean decrease by a factor of \(n\)?
Solutions. (1) The variance across exact conditional means estimates uncertainty in the population mean induced by the parameter distribution. The average conditional variance describes individual predictive spread at fixed parameters, averaged over that distribution. (2) The sample variance of the two fixed means is 2 \((\mathrm{mg/L})^2\). Each inner mean has variance \(1/2\), so the expected sample variance adds \(1/2\), giving \(2.5\;(\mathrm{mg/L})^2\). (3) Pooled outcomes vary both across parameter draws and among individuals conditional on a draw; the law of total variance includes both components. (4) Conditional independence of the \(n\) individuals under a common parameter value gives conditional variance \(v(\theta)/n\) for their sample mean. Dependence or a different sampling structure changes that result.