Diagnosing Monte Carlo Instability in Trial Simulations

Why Monte Carlo Stability Matters

A clinical trial simulation can produce a precise-looking mean, quantile, or probability even when the number of simulated trials is too small for that estimate to be numerically stable. Re-running the same simulation may then give a materially different result—not because the model or clinical assumptions changed, but because the random draws changed. This is Monte Carlo instability.

The distinction matters in nested simulations. An outer draw may represent uncertain model parameters or a virtual population, while inner draws generate patients, trial outcomes, or other random trial features conditional on that outer draw. The spread of the simulated trial outcomes is part of the predictive distribution. Monte Carlo error, by contrast, is the numerical error from approximating a target quantity with a finite number of simulation draws. A wide predictive distribution can be estimated stably; a narrow one can be estimated unstably.

The diagnostics below assume that the simulation code correctly represents the intended model, that runs use the same design and target quantity, and that the random draws used as independent units are independent (or that dependence is explicitly accounted for). A numerical stability check cannot establish that the model is clinically appropriate. It assesses whether finite simulation noise is small enough for the quantity being reported.

Replication and Batch Comparisons

A direct diagnostic is to repeat the complete simulation independently. Each replication should use new random draws, preserve the same model and design, and compute the same target using the same number of simulated units. If \(R\) independent estimates are \(\hat{\theta}_1,\ldots,\hat{\theta}_R\), their sample standard deviation describes variation between runs of that size. The estimated MCSE of their average is that standard deviation divided by \(\sqrt{R}\). For assessing one run, the between-run standard deviation itself is the relevant estimate of its run-to-run error.

When a single run is large enough, divide its independent simulation units into \(B\) batches. Compute the target separately in each batch. For a mean, the standard error across batch means estimates the MCSE of the overall mean, provided the batches are representative and the underlying units are independent. Batches are a diagnostic device, not a way to manufacture independence: if units are clustered by an outer parameter draw, batches should be formed at the independent outer-draw level, or the clustering must be reflected in the variance calculation.

A short derivation clarifies the formula. If the batch means are independent with common variance \(\tau^2\), then the variance of their average is

Replacing the unknown \(\tau\) with the observed sample standard deviation \(s_{\text{batch}}\) gives the batch-mean MCSE. If batch sizes differ, a simple average of batch means need not equal the overall mean; use size-weighted estimates and an appropriate variance calculation rather than applying this equal-batch formula mechanically.

Quantiles Need Their Own Stability Check

Quantiles are often more sensitive to finite simulation size than means, especially in the tails. A practical batch diagnostic is to calculate the target quantile separately in each of \(B\) equal batches, then calculate the sample standard deviation \(s_q\) of those batch quantiles. The quantity \(s_q/\sqrt{B}\) estimates the standard error of the average batch quantile. When batches are sufficiently large and representative, this serves as an approximation to the Monte Carlo error of the pooled quantile. It is a diagnostic approximation, not an identity: quantiles are nonlinear, and the average of batch quantiles need not equal the pooled quantile.

Independent full-run replications provide a useful cross-check: compute the same quantile in each complete run and compare the results. A density-based approximation is another option when the sampling distribution has a well-behaved, positive density near the target quantile, but it can be unreliable in sparse tails or near discontinuities. In practice, report the quantile’s units, the number and size of batches or replications, and the observed spread across them. A small estimated MCSE does not guarantee accuracy if all runs share a coding error or an unrepresentative sampling scheme.

Treatment-Effect Probabilities and Nested Draws

A treatment-effect probability is often estimated as the fraction of simulated trials meeting a specified criterion. For example, let \(\Delta\) be a simulated trial-level treatment effect and let \(c\) be a prespecified threshold. With \(N\) independent simulated trials, define \(I_i=1\) if trial \(i\) has \(\Delta_i>c\), and \(I_i=0\) otherwise. The estimate is \(\hat{p}=\sum_{i=1}^{N}I_i/N\). Since each indicator is Bernoulli with variance \(p(1-p)\), the estimated MCSE is

This formula assumes independent trials from the target simulation distribution. It can be misleading if many inner trials share the same outer parameter or virtual-population draw: outcomes within an outer group may be dependent after averaging over outer uncertainty. In that setting, calculate the event fraction within each independent outer unit or batch and estimate variability across those units. Do not count every inner trial as independent merely because it has a separate random number stream.

For an equal number \(K\) of inner simulations per outer draw, the variance of an estimate of an expectation can be understood using the law of total variance. If \(Y_{jk}\) is inner result \(k\) under outer draw \(j\), then the variance of the overall estimate has contributions from variation in the outer-conditional means and from inner simulation noise:

under independent outer draws, independent inner draws conditional on each outer draw, and common inner size \(K\). Here \(M\) is the number of outer draws. Increasing \(K\) reduces the second contribution but does not remove the first; more outer draws are needed to reduce uncertainty from sampling outer draws. This decomposition concerns numerical error in estimating an expectation, not the predictive spread of \(Y\).

Convergence Plots, Failure Modes, and Alternatives

A convergence plot displays an estimate as the number of simulated units increases. For a mean or event probability, plot the cumulative estimate against simulation count; for a quantile, do the same while noting that early tail estimates can be especially erratic. Add batch estimates or independently replicated estimates where possible. A stable-looking trace is supportive, not conclusive: plateaus can occur by chance, and a plot of the estimate alone does not quantify its Monte Carlo error. Pair it with batch or replication variability and the precision required for the intended use.

Useful warning signs include a cumulative estimate that continues to drift, substantial disagreement among batches or full runs, a quantile that jumps as a few extreme trials enter the sample, and a probability based on very few threshold exceedances. Also investigate failed or rejected simulation runs. Dropping failures can change the target distribution if failure is related to the simulated outcome; document the rule and establish that exclusions are appropriate before interpreting the result.

The main alternatives answer related but different questions. Independent full-run replication directly reveals run-to-run variability for the chosen simulation size and is easy to explain, but can be computationally expensive. Batch diagnostics reuse one large run and are economical, but depend on batches representing independent units and can be noisy when \(B\) is small. An analytic MCSE formula, such as the Bernoulli formula, is efficient when its independence and sampling assumptions hold. For dependent or clustered simulations, use cluster-level replication or a variance estimate that respects the dependence rather than an independent-unit formula. None of these methods substitutes for checking the simulation design, the target definition, or model assumptions.

A common mistake is to interpret a large predictive spread as evidence of Monte Carlo instability, or a small MCSE as evidence that the clinical prediction is trustworthy. The correction is to report these separately: describe predictive variability with the simulated distribution, and describe numerical stability with MCSEs or run-to-run comparisons. For example, a broad distribution of trial effects can coexist with a very stable estimated mean; conversely, a seemingly narrow distribution can have an unstable 99th percentile if only a few simulated trials determine its tail.