Monte Carlo Precision and Predictive Uncertainty Diagnostics

Predictive Uncertainty and Monte Carlo Error

A QSP simulation can produce a distribution of predicted clinical outcomes by propagating uncertain inputs and patient variation through the model. That distribution may be wide because the predicted outcomes genuinely vary across the target population or because model inputs are uncertain. Separately, a finite number of simulations means its summaries are estimated with limited numerical precision. This second component is Monte Carlo error.

The distinction matters when reporting a predictive mean, a quantile, or the probability of an event. A wide predictive distribution does not necessarily imply that its mean is imprecise: a sufficiently large simulation can estimate the mean precisely even when individual predictions vary substantially. Conversely, a narrow distribution can still have an imprecisely estimated tail quantile if few simulations represent that tail.

The usual standard error formulas below assume that the \(N\) outcomes are independent draws from the target predictive distribution and have equal weight. They describe simulation precision, not whether the QSP model is clinically credible or whether the target distribution adequately represents the intended population.

Monte Carlo Precision for a Predictive Mean

For outcomes with finite variance, the simulation estimates the predictive mean \(\mu_Y=E(Y)\) by the sample mean \(\bar{Y}\). The sample standard deviation \(s_Y\) describes the spread of simulated outcomes; the standard error of the estimated mean describes its Monte Carlo precision.

This estimate follows because the variance of the mean of independent draws is \(\operatorname{Var}(\bar{Y})=\operatorname{Var}(Y)/N\), and \(s_Y^2\) estimates \(\operatorname{Var}(Y)\). For a sufficiently large \(N\), a rough 95% Monte Carlo interval for the predictive mean is \(\bar{Y}\pm1.96SE(\bar{Y})\). It describes uncertainty in the numerical estimate of the mean, not the range containing 95% of individual predicted outcomes.

Precision of Event Probabilities

For a prespecified event \(A\), such as an outcome exceeding a threshold, define \(I_i=1\) when simulated outcome \(Y_i\) satisfies \(A\), and \(I_i=0\) otherwise. The estimated predictive event probability is the fraction of simulations meeting the event. Under independent draws, the event indicators are Bernoulli with probability \(p=P(A)\), so their sample proportion has an estimated standard error.

where \(K\) is the number of simulated outcomes satisfying \(A\). A large-sample approximate 95% Monte Carlo interval is \(\hat{p}\pm1.96SE(\hat{p})\). This approximation is less reliable when the event is rare or the number of event simulations is small; in that setting, a normal approximation may be poor. The issue is not cured simply by reporting more decimal places.

Precision of Predictive Quantiles

A predictive quantile \(q_p\) is defined by \(P(Y\leq q_p)=p\) for a continuous outcome distribution. The simulated quantile is obtained from the ordered outcomes. Unlike a mean or an event probability, its standard error depends on the density of outcomes near the target quantile: estimates are less stable where the distribution is sparse.

For independent draws from a continuous distribution with density \(f\) that is positive and sufficiently smooth near \(q_p\), a large-sample approximation is

In practice, \(f(q_p)\) must itself be estimated, and the approximation can be sensitive to that estimate, especially in a tail. A useful diagnostic is to divide the simulation into \(B\) independent batches of equal size, calculate the same quantile in each batch, and examine how much those batch quantiles vary. If their sample standard deviation is \(s_{\text{batch}}\), then \(s_{\text{batch}}/\sqrt{B}\) estimates the Monte Carlo standard error of their average. This is an approximation to the precision of the full-simulation quantile; the average of batch quantiles need not equal the quantile from all pooled outcomes.

Derivation and Practical Diagnostics

The mean and probability formulas follow directly from sampling variation across independent draws. For the mean, each draw has variance \(\sigma_Y^2\), so independence gives \(\operatorname{Var}(\bar{Y})=\sigma_Y^2/N\). Replacing the unknown \(\sigma_Y\) with \(s_Y\) and taking the square root yields \(SE(\bar{Y})=s_Y/\sqrt{N}\). For the event probability, each indicator has variance \(p(1-p)\); the mean of \(N\) indicators has variance \(p(1-p)/N\), estimated by replacing \(p\) with \(\hat{p}\).

The quantile approximation can be motivated by the empirical cumulative distribution function \(F_N\). Near \(q_p\), a small change \(\Delta q\) changes the true cumulative probability by approximately \(f(q_p)\Delta q\). The empirical cumulative probability at \(q_p\) has approximate variance \(p(1-p)/N\). Consequently, the quantile error has approximate variance \(p(1-p)/(Nf(q_p)^2)\), giving the stated standard error. This argument explains why low density near a quantile makes that quantile harder to estimate.

Useful diagnostics focus on the reported summary, not just the number of simulated outcomes. Compare estimates across independent batches, check whether the estimate changes materially as \(N\) increases, and report Monte Carlo standard errors alongside summaries when precision matters. For quantiles, inspect the spread of batch quantile estimates; for event probabilities, also report the event count so readers can see how much information supports the estimate.

Batch comparisons can reveal unstable estimates, but a reassuring result does not prove that a simulation is correct. All batches can share a coding error, omit an important source of uncertainty, or represent the wrong population. A stable summary demonstrates repeatability under the simulation procedure, not clinical validity.

Assumptions, Alternatives, and Common Pitfalls

The standard errors above rely on independent, equally weighted draws from the target predictive distribution. If simulated outcomes are dependent, the effective information can be less than \(N\), and the independent-draw formulas can understate Monte Carlo error. If outcomes have unequal weights, an unweighted sample standard error is generally not appropriate. In those cases, use independent replications of the full simulation procedure or a precision method that respects its dependence and weighting structure.

For means and event probabilities, analytic standard errors are fast and interpretable when their assumptions apply. Repeated batches provide a flexible empirical check and are especially useful for quantiles, but require enough independent batches and can give a noisy precision estimate when \(B\) is small. The density-based quantile formula is another option, but depends on estimating the density near the target quantile. These methods estimate simulation error; none replaces assessment of model assumptions or predictive uncertainty.

A frequent mistake is to report \(\bar{Y}\pm1.96s_Y\) as a 95% interval for the predictive mean. That uses the spread of individual outcomes instead of the standard error of their mean. The correction is to divide \(s_Y\) by \(\sqrt{N}\) when estimating the simulation precision of the mean. Conversely, a predictive interval for individual outcomes describes the predictive distribution and should not be interpreted as a Monte Carlo confidence interval for a summary.