Written and maintained by CASRAI Editorial Board
Last updated
Rubin’s rules are the formulas for turning m separate analyses of m multiply-imputed datasets back into one point estimate and one standard error — and the reason they exist is that averaging the point estimates is the easy part. Once a dataset with missing values has been multiply imputed (typically via MICE — multiple imputation by chained equations, or a joint-model equivalent), you have m completed datasets, not one. Fitting your analysis model separately on each gives m point estimates and m standard errors. Simply averaging the m point estimates is valid, but averaging the m standard errors is not: doing that throws away the disagreement between imputations, which is itself real information about how uncertain the missing values were. Rubin (1987) worked out the correct combination rule, splitting the pooled variance into a within-imputation component (ordinary sampling variance, present even with no missing data) and a between-imputation component (variance introduced specifically by not knowing the missing values). This guide works through both formulas, the small-sample degrees-of-freedom adjustment, and a fully reproducible worked example pooling five imputed datasets by hand. It assumes you already have m completed datasets in hand from a multiple-imputation procedure — whatever caused the missingness in the first place, whether that’s non-response, study attrition, or a general selection mechanism, is a separate question this guide doesn’t re-litigate; see this cluster’s research methods hub for the broader set of quantitative-analysis guides it sits alongside.
Why You Can’t Just Average the Standard Errors
Say your analysis model produces a point estimate Q (a regression coefficient, a mean difference, a log-odds ratio — Rubin’s rules work on any scalar estimate with an approximately normal sampling distribution) and a standard error for Q. Fit that same model on each of the m imputed datasets and you get Q1…Qm and their variances U1…Um (variance, not SE — Rubin’s rules combine on the variance scale, so square each SE before pooling and take the square root at the end).
The pooled point estimate is the simple mean:
Q̄ = (1/m) × Σ Qi
That part is uncontroversial — it’s an average of m unbiased estimates of the same quantity. The standard error is where a naive average goes wrong. If every imputed dataset agreed almost exactly on Qi, the missing data cost you very little extra uncertainty. If the m estimates are scattered widely, that scatter is telling you the missing values could plausibly have taken quite different values — and an honest standard error has to reflect that spread, not just the average of m within-dataset standard errors that each pretend the imputed values were the real ones.
The Two Variance Components, Combined Into a Third
Rubin’s rules define three quantities. The first two are computed directly from your m analyses; the third is the pooled total you actually report.
- Within-imputation variance (Ū) — the average of the m ordinary sampling variances, i.e. what your standard error would be if there had been no missing data at all: Ū = (1/m) × Σ Ui.
- Between-imputation variance (B) — the sample variance of the m point estimates around their own mean, i.e. how much disagreement the imputation process itself introduced: B = (1/(m−1)) × Σ (Qi − Q̄)2.
- Total variance (T) — Ū plus B, plus one further correction term B/m that accounts for using a finite number of imputations rather than an infinite number: T = Ū + B + B/m = Ū + (1 + 1/m) × B.
The pooled standard error is √T. Two things about T are easy to miss on a first read. First, B is scaled up by the factor (1 + 1/m), not added in raw — with only m=5 imputations that factor is 1.2, a meaningful inflation; with m=50 it’s 1.02, barely different from just adding B. This is exactly why the number of imputations matters (see below): a small m doesn’t just give you a noisier B, it also amplifies whatever B you get. Second, if B is zero — every imputed dataset happened to produce the identical estimate — T collapses to exactly Ū, which is the correct behaviour: if the missing values genuinely didn’t affect the estimate, Rubin’s rules shouldn’t inflate the standard error at all.
A useful diagnostic falls directly out of these two components: the fraction of missing information (λ), roughly λ ≈ (1 + 1/m) × B / T. It ranges from 0 (missingness had no effect on this estimate’s precision) to close to 1 (the estimate is dominated by imputation uncertainty). λ is specific to each estimate in your model, not a single number for the whole dataset — a covariate with heavy missingness can still produce a low-λ estimate for a parameter it barely influences, and a nearly-complete covariate can still produce a high-λ estimate if it’s central to the model. Note that Rubin’s rules combine standard errors on the variance scale regardless of what kind produced them — conventional OLS SEs, heteroscedasticity-robust SEs, or cluster-robust SEs can all be pooled the same way, as long as the same SE type is used consistently across all m imputed datasets.
Degrees of Freedom: Rubin (1987) and the Barnard-Rubin (1999) Small-Sample Fix
Because T is built from a finite, often small, number of imputations, the pooled estimate’s reference distribution for confidence intervals and p-values is a t-distribution, not a normal, and it needs its own degrees of freedom — separate from whatever degrees of freedom your analysis model would normally use. Rubin’s original (1987) formula is:
vold = (m − 1) / λ2
This works well when the complete-data sample is large relative to m, but it has an awkward property in small samples: it can produce a degrees-of-freedom value larger than the complete-data degrees of freedom would ever allow, which overstates precision. Barnard and Rubin (1999, Biometrika) proposed an adjusted degrees of freedom that caps vold using the complete-data degrees of freedom vcom (the degrees of freedom your analysis model would have if there were no missing data at all):
vobs = ((vcom + 1) / (vcom + 3)) × vcom × (1 − λ), then v̂ = (1/vold + 1/vobs)−1
v̂ is what current software actually reports as the pooled test’s degrees of freedom (Stata’s mi estimate, R’s mice::pool(), and SAS’s PROC MIANALYZE all implement the Barnard-Rubin adjustment, not the plain 1987 formula, specifically because of this small-sample problem). You don’t need to compute v̂ by hand — every MI package does it — but knowing it exists explains why a pooled p-value from a small study with several imputed covariates can look conservative (wider CI, higher p) compared to naively plugging √T into a z-test.
A Reproducible Worked Example: Pooling Five Imputed Datasets by Hand
To make the mechanics concrete, here is one complete pooling calculation from real, independently computed numbers rather than an illustrative round-number example. We simulated one dataset of n=30 subjects with a fully observed predictor X and an outcome Y following a true linear relationship (Y = 2 + 0.8X + error), then set roughly a third of Y values to missing at random. The missing values were multiply imputed m=5 times using stochastic regression imputation (fit on the observed cases, with fresh random noise drawn for each of the five completed datasets — this is what makes the five completed datasets actually differ from each other). The analysis model — a simple linear regression of Y on X — was then fit separately on each of the five completed datasets, using deterministic seeded random-number generation throughout so every figure below is exactly reproducible from the same seed (script and method in the Reproducibility Note at the end).
The five separate analyses of the slope on X produced:
| Imputation | Point estimate (Qi) | Standard error | Variance (Ui) |
|---|---|---|---|
| 1 | 0.7543 | 0.0375 | 0.00141 |
| 2 | 0.6999 | 0.0451 | 0.00203 |
| 3 | 0.7602 | 0.0500 | 0.00250 |
| 4 | 0.7873 | 0.0478 | 0.00228 |
| 5 | 0.7563 | 0.0380 | 0.00144 |
Pooling these five by hand:
- Pooled point estimate: Q̄ = (0.7543 + 0.6999 + 0.7602 + 0.7873 + 0.7563) / 5 = 0.7516
- Within-imputation variance: Ū = mean(U1…U5) = 0.0019
- Between-imputation variance: B = sample variance of the five Qi values around 0.7516 = 0.0010
- Total variance: T = Ū + (1 + 1/5) × B = 0.0019 + 1.2 × 0.0010 = 0.0032
- Pooled standard error: √T = 0.0561
- Fraction of missing information: λ ≈ 0.386
- Rubin (1987) degrees of freedom: vold = (5−1) / 0.3862 ≈ 26.85
A 95% confidence interval for the pooled slope is then 0.7516 ± t0.975,df × 0.0561, using whichever of vold or the Barnard-Rubin v̂ your software reports (with a complete-data df of 28 here, the two are close in this particular run since λ is moderate, not extreme). Notice that Ū alone (0.0436 on the SE scale, √0.0019) understates the true uncertainty — the between-imputation spread contributes real additional width, exactly the correction a naive average-the-SEs approach would have missed.
Pool on the Scale Where the Sampling Distribution Is Normal
Rubin’s rules assume Q’s sampling distribution is approximately normal on the scale you pool it on. A raw odds ratio, risk ratio, or hazard ratio is not approximately normal — it’s bounded at zero and right-skewed — but its log is. The standard practice is: fit the analysis model on each imputed dataset, extract the log-odds-ratio (or log-RR, log-HR) and its standard error on the log scale, pool those m log-scale values with Q̄/Ū/B/T exactly as above, and only exponentiate the pooled log-estimate and its confidence-interval bounds back to the ratio scale at the very end. Pooling five already-exponentiated odds ratios directly, then reporting their plain average as if it were the pooled odds ratio, is a common and quietly wrong shortcut — it both biases the point estimate (the mean of exponentials isn’t the exponential of the mean) and produces a standard error with no defensible interpretation. Any quantity with a naturally skewed or bounded distribution — variances, correlation coefficients (usually pooled via Fisher’s z-transform), R² values — needs the same treatment: transform to a scale where normality is a reasonable approximation, pool there, then back-transform only the final pooled estimate and interval.
How Many Imputations Should You Use?
Rubin’s original efficiency argument suggested that even a small number of imputations — m=3 to 5 — recovers most of the achievable precision for a point estimate under moderate missingness, and for years that was the conventional default in textbooks and software. Two things have shifted contemporary practice toward using more: first, while the point estimate itself is not very sensitive to m, the pooled standard error, degrees of freedom, and especially p-values close to a significance threshold are noticeably less stable at m=5 than at m=20-50 — rerunning the same imputation with a different random seed at low m can visibly move a borderline p-value. Second, computation is no longer the constraint it was when Rubin was writing in the 1980s, which removes the main practical reason to keep m small. The more directly useful benchmark is λ itself: an estimate with a higher fraction of missing information needs more imputations to pool stably than one with a low λ, since B is estimated from only m data points and a variance estimated from 5 points is far noisier than one estimated from 50. If you don’t know λ in advance, running a first pass at a moderate m (e.g. 20), checking the reported λ for your key estimates, and increasing m for any estimate with substantial missing information is a reasonable iterative approach, rather than committing to a fixed m before seeing how much information was actually missing. Whatever m you settle on, record the choice and the imputation model in the trial or study’s statistical analysis plan before unblinding or final analysis, not after — choosing m post hoc based on which value produces a more favourable p-value is a form of analytic flexibility reviewers increasingly ask about directly.
Rubin’s Rules in Stata, R, and Python
- Stata: after
mi imputehas created the m imputed datasets and the data aremi set, run your analysis model throughmi estimate: regress y x(substituting any estimation command) — Stata fits the model on every imputed dataset automatically and reports the pooled coefficient, pooled SE, and the Barnard-Rubin degrees of freedom without any manual combination step. - R: with the
micepackage, impute viamice(data, m = 20), fit the analysis model on the imputed object withwith(imp, lm(y ~ x)), then pool withpool(fit)—summary(pool(fit))returns the pooled estimate, SE, t-statistic, Barnard-Rubin df, and p-value directly. - Python:
statsmodels’MICEData/MICEclasses handle imputation and pooling together (MICE(model_formula, model_class, mice_data).fit()); for imputations produced elsewhere, the pooling itself is a short manual calculation of Q̄, Ū, B, and T following the formulas above, since there is no single dedicated pooling function as commonly used as R’spool().
In every case, check that the software is actually applying Rubin’s rules across your m imputations and not, e.g., stacking the m completed datasets into one long dataset and running a single analysis on the combined rows — that produces a dataset m times too large, artificially shrinks every standard error, and is a different (incorrect) error from the naive-SE-averaging mistake this guide opened with.
Frequently Asked Questions
What is the difference between within-imputation and between-imputation variance?
Within-imputation variance (Ū) is the ordinary sampling variance you’d have even with no missing data — it’s the average of the m separate analyses’ own standard errors, squared. Between-imputation variance (B) is variance specifically caused by not knowing the missing values — how much the point estimate itself shifts depending on which imputed dataset you happened to analyze. Total variance combines both, with an extra correction term for using a finite m.
Can I just average the standard errors from each imputed dataset instead of using Rubin’s rules?
No. Averaging the SEs (or their squares) directly captures only the within-imputation variance and ignores the between-imputation variance entirely — it systematically understates the true uncertainty whenever the imputed datasets disagree with each other, which is precisely the situation multiple imputation was designed to represent honestly.
Why is the between-imputation variance multiplied by (1 + 1/m) instead of just added?
The extra 1/m term corrects for using a finite number of imputations to estimate B rather than an infinite number. With few imputations, B itself is a noisier estimate of the true between-imputation variance, so the correction inflates T slightly to account for that added uncertainty; as m grows, the (1 + 1/m) factor shrinks toward 1 and the correction becomes negligible.
Do Rubin’s rules apply to any statistic, or only regression coefficients?
They apply to any scalar estimate with an approximately normal sampling distribution — means, proportions (on the logit scale), regression coefficients, and log-transformed ratios all qualify. They are not valid, without modification, for statistics whose sampling distribution is not approximately normal on the scale being pooled, which is exactly why ratio-type estimates get pooled on the log scale rather than the raw ratio scale.
Is pooling multiply-imputed estimates the same as pooling effect sizes in a meta-analysis?
No, even though both are called “pooling” and both combine multiple estimates of the same underlying quantity into one. Rubin’s rules combine m analyses of the same sample across m versions of one imputed dataset. Meta-analytic pooling combines effect estimates from different studies, each with its own sample, and weights them by between-study heterogeneity rather than by missing-data uncertainty — the Mantel-Haenszel method for pooling stratum-specific odds ratios is one example of that different family. The formulas are unrelated even though the goal — one combined estimate and one combined standard error — sounds similar.
Reproducibility Note
The worked example above is not an illustrative round-number example — every value was independently computed via a seeded Node.js simulation: a Mulberry32 pseudorandom generator (fixed seed) drove Box-Muller normal draws to build the n=30 base dataset (Y = 2 + 0.8X + N(0,4) error), induce ~33% missingness in Y, fit a stochastic regression imputation model on the observed cases, generate m=5 completed datasets (each with independently drawn imputation noise), fit ordinary least squares on each, and combine the five slope estimates and variances via the exact Ū/B/T formulas given above. The full script and seed are available on request; every number in the worked-example table and the pooling calculation beneath it is exactly reproducible from that same seed, not asserted or rounded from an external source.
