Monte Carlo Error and Simulation Size

Why Simulation Size Matters

In the parametric g-formula, a fitted set of longitudinal models can be used to simulate counterfactual histories under a specified treatment regime. Averaging an outcome summary over those simulated histories approximates the corresponding model-based expectation. Because the simulated histories are random, the average changes if the random-number seed changes, even when the data and fitted models stay fixed.

That variation is Monte Carlo error. It is not the same as sampling uncertainty: sampling uncertainty arises because the observed study data, and hence the fitted models and target estimate, would change in another sample. This tutorial quantifies the first source of variation and explains how to make it small enough for the calculation at hand.

Let \(D\) denote the observed data and let the fitted models be held fixed conditional on \(D\). For a chosen regime \(g\), write \(Z_{im}(g)\) for the outcome summary from the \(m\)th simulated trajectory associated with baseline individual \(i\). The summary may be a scalar quantity of interest; the simulation may generate a full longitudinal history before that summary is calculated. Suppose there are \(n\) baseline individuals and \(M\) independent simulated trajectories per individual. The simulated estimate is

Define \(\mu_i(g)=E[Z_{im}(g)\mid D]\) and \(\sigma_i^2(g)=\operatorname{Var}(Z_{im}(g)\mid D)\). The corresponding infinite-simulation, model-based estimate for these baseline individuals is

This is the quantity the simulation is approximating. It is still based on the observed data and fitted models; “infinite simulation” removes only the random error from approximating its expectation.

Conditional Monte Carlo Error

For a simple variance formula, assume that, conditional on \(D\), simulated trajectories are independent across \(m\) and \(i\); that each individual’s trajectories have the same conditional mean \(\mu_i(g)\) and variance \(\sigma_i^2(g)<\infty\); and that \(n\) and the baseline individuals are fixed during simulation. Independence here concerns whole simulated trajectories, not the time points within a trajectory. Outcomes at different times may be dependent within a trajectory.

Proof. By linearity of conditional expectation,

For the variance, conditional independence makes the covariance between distinct simulated trajectories zero. Therefore,

Taking the square root gives the MCSE. The factor \(1/M\) in the variance means that multiplying the simulation size by four divides the MCSE by two, provided the simulation design and conditional variances stay the same. This is a precision trade-off, not a guarantee that an estimate is close to the causal quantity of interest: model misspecification and other sources of error are not removed by increasing \(M\). \(\square\)

Estimating the MCSE and Choosing \(M\)

The conditional variances \(\sigma_i^2(g)\) are generally unknown, but they can be estimated from the simulated values. Let

For \(M\geq 2\), the conditional sample variance \(S_i^2(g)\) is unbiased for \(\sigma_i^2(g)\) when the simulated values are independent and identically distributed conditional on \(D\), with finite variance. Substitution gives the usual estimated MCSE:

To target an MCSE no larger than \(\delta\), use a pilot simulation to estimate the variance component and then choose \(M\) accordingly. If the same number of trajectories is used for every baseline individual, the plug-in calculation is

The pilot estimate is uncertain, so this calculation is a planning guide rather than a certainty about the final MCSE. A practical workflow is to run a pilot, calculate the estimated MCSE, increase \(M\) if it exceeds the target, and verify the MCSE in the final run. If several regimes or outcome summaries are reported, choose a target and check the simulation error for each reported quantity, not just one convenient estimate.

Repeated Runs and Data Uncertainty

Repeated simulation runs are a direct diagnostic for Monte Carlo variation. Keep the data, fitted models, baseline individuals, regime, and simulation size fixed; then repeat the simulation with independent random numbers. The spread of the resulting estimates reflects simulation variability under those fixed conditions. For \(R\geq2\) independent runs with estimates \(\hat{\psi}_{g,M}^{(1)},\ldots,\hat{\psi}_{g,M}^{(R)}\), calculate their sample standard deviation:

This estimates the MCSE for that simulation design, subject to the variability from having only \(R\) runs. A handful of runs can reveal an obvious instability, but their sample standard deviation may itself be imprecise. Repeating runs does not cause the underlying data or fitted models to vary, so it cannot estimate their sampling uncertainty.

Proof. Add and subtract \(E[\hat{\psi}_{g,M}\mid D]\) inside the deviation from the overall mean. Equivalently, write the centered estimate as the sum of two terms:

On squaring and taking expectations, the cross-product is zero. Conditional on \(D\), the first bracket has conditional expectation zero, while the second is fixed. The expected square of the first bracket is \(E[\operatorname{Var}(\hat{\psi}_{g,M}\mid D)]\); the expected square of the second is \(\operatorname{Var}(E[\hat{\psi}_{g,M}\mid D])\). Their sum is the stated identity. \(\square\)

When the conditional simulation estimate is unbiased for \(\hat{\psi}_{g,\infty}\), its conditional mean is \(\hat{\psi}_{g,\infty}\). The first term in the decomposition is then variation of that model-based estimate across data sets; the second is the average Monte Carlo variance. The two terms answer different questions and should not be reported as if they were interchangeable.

Diagnostics, Alternatives, and Limitations

A useful diagnostic is to compare the reported MCSE with the scale of the estimate and with the precision needed for the analysis. Also check that the MCSE decreases as \(M\) increases: under independent simulation and finite variance, doubling \(M\) should reduce it by a factor of \(1/\sqrt{2}\), approximately. A result that changes substantially across independent runs, or an estimated MCSE that remains large relative to the desired precision, indicates that more simulation may be needed.

The simple formula assumes independent simulated trajectories. If simulations are correlated—for example, because random numbers are reused across runs or trajectories—the covariance terms do not vanish, and the formula can understate or overstate the actual MCSE. Estimate variability using independent complete runs, or use a variance calculation that accounts for the dependence created by the design. For a contrast between regimes simulated with shared random numbers, calculate uncertainty for the paired contrast itself; separate marginal MCSEs do not include its covariance.

Monte Carlo integration is useful when the required expectation is difficult to compute directly. If exact or deterministic integration is feasible, it avoids random simulation error, although it may require other computational approximations. Monte Carlo error is only one limitation: increasing \(M\) does not repair misspecified conditional models, unsupported extrapolation, or a mismatch between the target population and the baseline individuals being averaged.

A common mistake is to treat a small MCSE as evidence that the causal estimate is precise or unbiased. The reasoning is wrong because the MCSE measures only variation from the simulated draws conditional on the fitted analysis. The correction is to assess simulation error separately from data and model uncertainty, and to report what the MCSE actually quantifies.