Choosing the Number of Predictive Simulations

Why the Number of Simulations Matters

A predictive check compares observed data with data generated under a fitted model. The simulated datasets provide a numerical approximation to the model’s predictive distribution; that approximation is itself subject to random variation. With too few simulations, a plotted median or quantile may move noticeably if the simulation is repeated, even when the model and observed data have not changed. Choosing the replicate count is therefore part of making the check stable enough to interpret.

The relevant count is the number of independent simulated replicates contributing to the quantity being estimated. A replicate is typically a simulated dataset generated under the intended design and model. Measurements within one dataset can be dependent, for example because they come from the same simulated subject. They should not automatically be counted as independent replicates when assessing the precision of a statistic calculated across datasets.

MCSE is not the uncertainty in the model parameters, nor is it the sampling uncertainty that would arise from collecting a different study dataset. It quantifies only the numerical variability introduced by finite simulation. A small MCSE makes the predictive reference more reproducible; it does not establish that the model is adequate.

Monte Carlo Precision for Quantiles

Suppose \(R\) independent simulation draws are used to estimate the \(p\)-quantile \(q_p\) of a continuous predictive distribution with cumulative distribution function \(F\) and density \(f\). Assume \(0<p<1\), \(f(q_p)>0\), and that the density is reasonably smooth near \(q_p\). For sufficiently large \(R\), the empirical quantile has approximate variance

Here \(\hat q_p\) is the quantile estimated from the simulations. The formula shows two important features. First, MCSE decreases in proportion to \(1/\sqrt{R}\): multiplying the replicate count by four approximately halves the MCSE. Second, quantiles are less precise where the predictive density is low. A small change in cumulative probability near such a quantile can correspond to a large change in the outcome value.

This approximation concerns independent draws from the predictive distribution, or independent replicate-level values when the target is a statistic calculated once per simulated dataset. If observations are clustered within simulated subjects or datasets, the dependence structure matters. The formula should not be applied as though every measurement were independent when the estimate is based on clustered measurements.

Why the Quantile Formula Works

A short approximation explains the formula. At the true quantile, \(F(q_p)=p\). The empirical cumulative distribution function at \(q_p\), denoted \(\hat F_R(q_p)\), is an average of \(R\) indicator variables: each is one if its draw is at or below \(q_p\), and zero otherwise. Each indicator has success probability \(p\), so

For a small displacement \(\delta\) from \(q_p\), smoothness of the density gives \(F(q_p+\delta)\approx p+f(q_p)\delta\). Thus, to first order, an error in cumulative probability is converted to an error in the quantile by dividing by \(f(q_p)\). Applying that local conversion to the variance above gives

This is a large-sample approximation, not an exact guarantee for every finite simulation. It can be unreliable for very small \(R\), very sparse tails, a density close to zero near the target quantile, or discrete predictive outcomes with many ties.

Simulation Counts and NPDE

Normalized prediction distribution errors (NPDE) use simulated predictive distributions to express observed deviations on a standardized, decorrelated scale. Their detailed transformation is treated in the later material on NPDE. For choosing the number of simulations, the important point is that the predictive distribution’s finite-simulation estimates—such as its location, covariance, or probability assigned to an observed value—can vary from one simulation run to another. That variation can carry through to the calculated NPDE.

The familiar quantile MCSE formula is useful for a probability or quantile estimated from independent simulation draws, but it does not by itself give the MCSE of every NPDE calculation. In particular, covariance estimates and transformations involving several correlated observations can have more complicated finite-simulation behavior. A simple illustration is the estimated fraction of simulated values below an observation. If the true predictive probability is \(p\), its Monte Carlo standard error is approximately \(\sqrt{p(1-p)/R}\). A normal-scale score based on that estimated probability will also vary; the amount depends on the probability’s location and on the transformation. Tail probabilities deserve care because a small probability error can produce a substantial score change.

Checking Whether a Pattern Is Simulation Noise

A predictive check can show irregular wiggles even when its underlying target is stable. Conversely, a persistent discrepancy may remain after simulation noise has been reduced. A practical assessment asks whether the feature is reproducible across independent runs and whether its size is plausible relative to the MCSE.

Look especially carefully at extreme quantiles, where few simulated values determine the estimate, and at intervals or time regions with limited information. Quantile estimates that jump between batches, disappear as the replicate count rises, or change direction repeatedly are warning signs of simulation noise. Systematic departures that persist across independent batches are less easily explained by Monte Carlo variability alone, though they still require scientific interpretation.

Choosing a Count and Avoiding Common Errors

A useful planning rule follows directly from the quantile formula. If an approximate MCSE no larger than \(\delta\) is desired, then

The target \(\delta\) should be chosen in the units and context of the check: for example, small relative to a clinically meaningful concentration difference or to the visual resolution needed to distinguish a trend. If \(f(q_p)\) is unknown, a pilot simulation can provide a provisional density estimate, followed by a larger run and a stability check. The final count should be large enough for the key summaries, not merely convenient for one central statistic.

Increasing the number of simulations reduces Monte Carlo error, but it does not repair model misspecification, biased parameter estimates, an inappropriate predictive target, or a poorly chosen analysis. It also costs computation. Alternatives or complements include comparing independent batches, estimating MCSE directly for the summary of interest, or using a justified analytic calculation where available. These approaches diagnose precision; they do not make finite simulation error vanish.

• Do not confuse simulation replicates with observations. Measurements from the same simulated subject or dataset may be dependent; treating them as independent can understate uncertainty.
• Do not assume all quantiles have equal precision. The density near the quantile matters, so a tail quantile can be substantially less stable than the median.
• Do not call Monte Carlo error model uncertainty. More simulations stabilize the numerical approximation under the specified model; they do not account for uncertainty about the model or its parameters.

The central practical lesson is to make simulation precision commensurate with the question. A replicate count that is adequate for a median may be inadequate for a tail quantile or a simulation-based normal score. Precision should be evaluated for the specific summaries used in the check, while keeping Monte Carlo variability distinct from uncertainty and model misfit.