Monte Carlo Error In The Pooled Point Estimate

Why the Pooled Estimate Has Monte Carlo Error

With missing data, different completed data sets can produce different estimates of the same quantity. Rubin’s pooled point estimate is their average. If only a finite number of imputations is generated, that average has some simulation variability: generating a different finite set of imputations could produce a slightly different pooled estimate, even when the observed data and imputation procedure stay fixed. This variability is the Monte Carlo error of the pooled point estimate.

Monte Carlo error is not the same as uncertainty about the underlying population quantity. The latter is represented by the pooled standard error, which combines within-imputation and between-imputation uncertainty. Monte Carlo error instead describes how precisely the finite set of imputations has approximated the average over the imputation distribution. Both matter, but they answer different questions.

Definition: Let \(Q_l\) be the estimate of the target quantity from imputation \(l\), for \(l=1,\ldots,m\). The pooled point estimate is \(\bar{Q}=m^{-1}\sum_{l=1}^m Q_l\). Its Monte Carlo error is the standard deviation, across repeated finite sets of imputations under the same imputation procedure, of this pooled average. For approximately independent imputations, it is estimated by \(\sqrt{B/m}\), where \(B=(m-1)^{-1}\sum_{l=1}^m(Q_l-\bar{Q})^2\) is the between-imputation variance.

The usual formula requires care about what is being held fixed. It concerns the simulation variability of \(\bar{Q}\), conditional on the observed data and the imputation model. It is not a measure of whether that model is correct, whether the missing-data assumptions are plausible, or whether the analysis targets the right estimand. It also relies on having draws that adequately represent the intended imputation distribution.

Deriving the Usual Approximation

Suppose, conditional on the observed data and the imputation procedure, the estimates \(Q_1,\ldots,Q_m\) are independent draws with a common finite variance \(\sigma_Q^2\). The variance of their average is then

\operatorname{Var}(\bar{Q})= \operatorname{Var}\left(\frac{1}{m}\sum_{l=1}^{m}Q_l\right) =\frac{1}{m^2}\sum_{l=1}^{m}\operatorname{Var}(Q_l) =\frac{\sigma_Q^2}{m}.

The middle equality uses independence: the covariances between different \(Q_l\)'s are zero. The sample variance \(B\) estimates \(\sigma_Q^2\), so substituting \(B\) gives the estimated variance \(B/m\) and estimated standard deviation \(\sqrt{B/m}\). This is an estimated simulation standard error, not the actual error in a particular run. The actual deviation of \(\bar{Q}\) from the infinite-imputation average is unknown.

The pooled standard error comes from a different variance calculation. Let \(U_l\) be the estimated variance of \(Q_l\) within completed data set \(l\), and let \(\bar{U}=m^{-1}\sum_{l=1}^{m}U_l\). Rubin’s total variance estimate is \(T=\bar{U}+(1+1/m)B\); the pooled standard error is \(\sqrt{T}\). The within-imputation component reflects uncertainty from analyzing a completed data set. The between-imputation component reflects variation across completed data sets, with a finite-\(m\) adjustment. By contrast, the approximate Monte Carlo error of the point estimate is only \(\sqrt{B/m}\). It does not include \(\bar{U}\).

Worked Example: Computing Monte Carlo Error and the Pooled Standard Error

Suppose \(m=25\) imputations produce estimates of a regression slope measured in outcome units per one-unit increase in a predictor. The estimates have \(\bar{Q}=2.40\) outcome units per predictor unit, \(B=0.040\) squared outcome units per squared predictor unit, and \(\bar{U}=0.360\) in the same squared units. Treat the estimates as approximately independent draws for this calculation.

The approximate Monte Carlo error is

\sqrt{\frac{B}{m}} =\sqrt{\frac{0.040}{25}} =\sqrt{0.0016} =0.0400

Thus, the estimated simulation standard deviation of the pooled slope due to using 25 imputations is \(0.0400\) outcome units per predictor unit. The pooled variance and pooled standard error are instead

T=\bar{U}+\left(1+\frac{1}{m}\right)B =0.360+\left(1+\frac{1}{25}\right)(0.040) =0.360+0.0416=0.4016, $$ $$ \sqrt{T}=\sqrt{0.4016}=0.6337.

The pooled standard error, \(0.6337\) outcome units per predictor unit, describes inferential uncertainty under the pooling framework. The Monte Carlo error, \(0.0400\) in the same units, describes simulation variability in the pooled estimate due to the finite number of imputations. The latter is about \(0.0400/0.6337=0.0631\), or 6.31%, of the pooled standard error. That comparison helps describe its scale; it does not make the two quantities interchangeable.

Assessing the Size in Context

A Monte Carlo error is most useful when reported in the same units as its estimate and compared with a relevant scale. Possible reference points include the pooled standard error, the precision needed for the scientific interpretation, or a decision threshold relevant to the analysis. A small ratio to the pooled standard error suggests that simulation noise is relatively small compared with inferential uncertainty. It does not establish that the imputation model is appropriate or that the analysis is unbiased.

The absolute size also matters. An error of \(0.6\) might be negligible for an estimate measured in hundreds, but consequential for a threshold or effect measured in fractions of a unit. Rounding can hide this issue: a result reported to two decimal places may appear stable even if the Monte Carlo error is large enough to change the last reported digit. Conversely, a small difference between two runs does not prove that the imputation procedure has negligible Monte Carlo error; random differences can happen to be small.

Worked Example: Interpreting Monte Carlo Error on an Effect Scale

In a hypothetical analysis, \(m=40\) imputations give a pooled estimated mean difference of \(1.20\) days, between-imputation variance \(B=16\) days squared, and average within-imputation variance \(\bar{U}=9\) days squared. Assume the imputation estimates are approximately independent.

The Monte Carlo error of the pooled mean difference is

\sqrt{\frac{B}{m}}=\sqrt{\frac{16}{40}} =\sqrt{0.400}=0.6325\text{ days}.

The pooled standard error is

\sqrt{T} =\sqrt{\bar{U}+\left(1+\frac{1}{m}\right)B} =\sqrt{9+\left(1+\frac{1}{40}\right)16} =\sqrt{25.4}=5.040\text{ days}.

Here, the Monte Carlo error is \(0.6325/5.040=0.1255\), or 12.55% of the pooled standard error. It is also a substantial fraction of the point estimate: \(0.6325/1.20=0.5271\). The two comparisons communicate different things. The first compares simulation noise with inferential uncertainty; the second compares it with the estimated effect itself. Neither ratio is a probability that the estimate is wrong, and neither alone answers whether the effect is scientifically important.

Dependence and Diagnostics

The simple formula \(\sqrt{B/m}\) assumes that the imputation-specific estimates behave like independent draws. In chained-equation procedures, iterations can be dependent: values updated in one sweep may be related to values in later sweeps. If the \(Q_l\)'s have positive serial correlation, averaging them reduces variability less than the independent-draw formula suggests.

For a stationary sequence with variance \(\sigma_Q^2\) and lag-\(k\) correlation \(\rho_k\), the variance of its average is

\operatorname{Var}(\bar{Q}) =\frac{\sigma_Q^2}{m} \left[1+2\sum_{k=1}^{m-1}\left(1-\frac{k}{m}\right)\rho_k\right].

This follows by expanding the variance of the sum: there are \(m\) variance terms and, for each lag \(k\), \(m-k\) covariance pairs, each counted twice. When the correlations are positive, the bracketed factor exceeds one, so the independent-draw approximation understates Monte Carlo error. In practice, assessing dependence may require examining trace plots or autocorrelation for the quantities being pooled and using suitable chain-based methods, such as an effective sample size or batch means. A single \(B\) value describes spread across imputations; it does not by itself reveal the order or dependence structure of the \(Q_l\)'s.

Worked Example: Allowing for Serial Dependence

Consider a hypothetical stationary sequence of \(m=30\) imputation-specific estimates of a treatment effect, with estimated marginal variance \(B=0.090\) squared outcome units. Suppose their lag correlations are adequately represented by an AR(1) pattern, \(\rho_k=0.5^k\). The independent-imputation calculation would give

\sqrt{\frac{B}{m}}=\sqrt{\frac{0.090}{30}} =\sqrt{0.003}=0.05477

outcome units. Under the stated correlation pattern, the variance inflation factor for a finite sequence is

1+2\sum_{k=1}^{29}\left(1-\frac{k}{30}\right)(0.5)^k \approx 2.8667.

Therefore, the dependence-adjusted approximate Monte Carlo error is

\sqrt{\frac{0.090}{30}(2.8667)} =\sqrt{0.008600} =0.09274

outcome units. Under this illustrative dependence pattern, the independent formula would report a noticeably smaller error. This calculation is only as credible as the stationarity and correlation assumptions; real chains need diagnostics rather than an assumed AR(1) pattern.

Common Mistakes and Limits

Calling \(\sqrt{B/m}\) the pooled standard error. It omits within-imputation variance. Use \(\sqrt{T}\) for the pooled standard error and describe \(\sqrt{B/m}\) as approximate Monte Carlo error of the point estimate.
Treating Monte Carlo error as an observed discrepancy. It estimates the standard deviation of the pooled estimate over repeated finite imputation runs; it does not reveal the exact difference between the current estimate and an infinite-imputation value.
Assuming a small Monte Carlo error validates the analysis. It only addresses finite simulation precision under the imputation procedure. It cannot diagnose misspecified imputation models, unsuitable assumptions about missingness, or a mismatch between the estimand and analysis.
Ignoring dependence. A long sequence of highly correlated estimates does not provide as much independent information as the same number of independent estimates. Check behavior across the sequence and use methods that account for dependence when needed.

The main alternative to the plug-in calculation is a direct assessment of simulation variability across repeated runs or blocks of a chain. Such checks can reveal instability that a single \(B\) value obscures, but they also require adequate runs and careful treatment of dependence. None changes the interpretation: Monte Carlo error concerns computational approximation, while the pooled standard error concerns statistical uncertainty. Model checking and sensitivity analysis address different limitations and remain necessary.

Key takeaway: For approximately independent imputations, estimate the Monte Carlo error of the pooled point estimate by \(\sqrt{B/m}\). Keep it distinct from Rubin’s pooled standard error \(\sqrt{\bar{U}+(1+1/m)B}\), report it in context, and investigate dependence before relying on the independent-imputation approximation.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For \(m=16\) and \(B=0.25\) squared units, what is the approximate Monte Carlo error, in the units of \(Q\)?
Solution: \(\sqrt{B/m}=\sqrt{0.25/16}=\sqrt{0.015625}=0.125\) units.
If \(\bar{U}=4\), \(B=1\), and \(m=9\), calculate the pooled standard error and Monte Carlo error. Are they the same?
Solution: \(T=4+(1+1/9)(1)=4+1.1111=5.1111\), so \(\sqrt{T}=2.261\). The Monte Carlo error is \(\sqrt{1/9}=0.3333\). They are not the same: the pooled standard error also includes within-imputation uncertainty.
What does a Monte Carlo error of \(0.08\) estimate? Does it mean there is an 8% probability that the pooled estimate is wrong?
Solution: It estimates the standard deviation, in the estimate’s units, of the pooled average across repeated finite imputation runs under the stated procedure. It is not a probability and does not mean the estimate is wrong with any particular probability.
Why can positive autocorrelation make \(\sqrt{B/m}\) too small?
Solution: The independent formula sets all covariances between \(Q_l\)'s to zero. With positive autocorrelation, the variance of the average includes positive covariance terms, so its standard error is larger than the independent-draw calculation.
If Monte Carlo error is tiny relative to the pooled standard error, does that establish that the imputation model is valid?
Solution: No. It suggests that finite simulation noise in the pooled estimate is small relative to the estimated inferential uncertainty. It does not assess model specification, missingness assumptions, or the appropriateness of the estimand.