Multiple Imputation

September 2026 · 12 min read

Suppose we want a population mean, but some sampled measurements are missing. A single completed table lets us calculate a mean. An ordinary complete-data standard error then treats those filled entries as observed, leaving out uncertainty about what could have occupied the gaps.

Multiple imputation carries several plausible completed datasets through the analysis. If different plausible fills lead to different answers, that disagreement becomes part of the reported uncertainty. This works under suitable imputation and analysis assumptions; randomizing cells alone is not enough.

Impute, analyze, then pool

For example, fit the same regression to every completed dataset and pool a particular coefficient. Do not first average the filled cells and fit only one regression: that generally loses the variation we wanted to carry forward. Nor should you average p-values or treat the completed datasets as independent new samples.

See the two sources of uncertainty before the formula

Imagine three completed-data analyses report estimates 4, 5, and 6 for the same population quantity. Each reports a standard error of 0.5, so each estimated variance is 0.25. These are constructed analysis outputs for learning the pooling calculation.

The pooled point estimate is their average, 5. The usual sampling uncertainty remains: averaging the three variances gives 0.25. We also need to account for the estimates disagreeing. Their sample variance is 1, since their deviations from 5 are −1, 0, and 1.

We call these the within-imputation variance and between-imputation variance. “Within” means uncertainty reported by each analysis; it is not the variance of the raw feature values themselves.

Rubin’s rule adds those components

Let $m$ be the number of imputations, $\bar U$ the average estimated variance within the analyses, and $B$ the sample variance between their estimates. The total estimated variance, $T$, is

The extra $B/m$ accounts for using finitely many imputations. In our example, $T=0.25+(1+1/3)\times1\approx1.5833$. The pooled standard error is its square root, approximately 1.2583, rather than the complete-data standard error of 0.5.

If the estimates agree, the between component is zero, but ordinary sampling uncertainty remains. Increasing the number of imputations reduces Monte Carlo error; it does not make missing information disappear.

Plausible randomness has a statistical model behind it

A stochastic imputer needs to reflect what is uncertain about the unobserved measurements. For a normal model, that includes both the variation of individual measurements and uncertainty in the model’s fitted mean and variance. Holding estimated parameters fixed and merely adding arbitrary noise can miss part of that uncertainty.

The same point applies to resampling: randomly copying observed values into gaps creates different tables, but is not by itself a justification for valid multiple-imputation inference. Proper procedures specify how the draws account for the information available and how they relate to the intended analysis.

import numpy as np

# Teaching model: independent normally distributed measurements;
# missingness does not depend on the measurement (MCAR).
raw = np.array([2., 4., np.nan, 6., np.nan, 8.])
missing = np.isnan(raw)
observed = raw[~missing]
n_obs = len(observed)
if n_obs < 2 or observed.var(ddof=1) <= 0:
    raise ValueError("this example needs at least two varying observations")
rng = np.random.default_rng(42)
m = 100
estimates, variances = [], []

for _ in range(m):
    # Reference prior: flat in the mean, inverse in the variance.
    # Draw population variance, then population mean, then missing values.
    variance = ((n_obs - 1) * observed.var(ddof=1)
                / rng.chisquare(n_obs - 1))
    mean = rng.normal(observed.mean(), np.sqrt(variance / n_obs))
    completed = raw.copy()
    completed[missing] = rng.normal(
        mean, np.sqrt(variance), size=missing.sum()
    )
    assert np.array_equal(completed[~missing], observed)
    estimates.append(completed.mean())
    variances.append(completed.var(ddof=1) / len(completed))

q = np.asarray(estimates)
u = np.asarray(variances)
pooled = q.mean()
within = u.mean()
between = q.var(ddof=1)
total = within + (1 + 1 / m) * between
print(f"pooled mean: {pooled:.3f}")
print(f"pooled standard error: {np.sqrt(total):.3f}")

pooled mean: 4.971
pooled standard error: 1.889

What makes the resulting inference defensible?

Specify the quantity to estimate and its analysis model before choosing the imputation model. Preserve important relationships, including interactions, nonlinear terms, group structure, and variable bounds when they matter. Examine whether observed and imputed values are plausible in comparable subgroups.

Standard applications often assume MAR conditional on the included observed variables. That assumption and adequate model specification do substantive work. Multiple imputation can also be part of an MNAR sensitivity analysis, but only after explicitly introducing assumptions about how the hidden values differ. Pooling alone does not correct an unspecified selection process.

For chained equations, inspect chain behavior and the stability of the quantities you intend to report. Choose the number of imputations based on the missing information and desired Monte Carlo precision; a fixed “5 to 50 always suffices” rule is too coarse. More datasets help stabilize simulation, not cure model misspecification.

Keep prediction and inference workflows distinct

For predictive modeling, fit the imputation procedure inside each training fold and reproduce the intended missing-input behavior on validation rows. Several imputed training sets can lead to an ensemble, but assess its predictive accuracy, calibration, and cost directly. AUC, accuracy, and neural-network weights are not automatically suitable inputs to the scalar pooling formula above.

For inference on an existing study, the workflow need not involve a train/test split at all. An observed study outcome can appropriately help impute covariates when the imputation and analysis models are compatible. In a deployed predictor, an unknown future outcome is unavailable as an imputer input. The difference is what information the task permits.

C++: pool estimates that have already been analyzed

The following function handles the pooling arithmetic only. Its inputs are per-completion estimates and variances, not raw filled values or standard errors. Keeping this step separate makes it clear that a correct formula cannot rescue unsuitable imputation draws.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

struct Pooled {
    double estimate, within, between, totalVariance;
};

Pooled pool(const std::vector<double>& estimates,
            const std::vector<double>& variances) {
    if (estimates.size() < 2 || estimates.size() != variances.size())
        throw std::invalid_argument("matching arrays and at least two imputations required");
    const double m = static_cast<double>(estimates.size());
    double qbar = 0, within = 0;
    for (std::size_t j = 0; j < estimates.size(); ++j) {
        if (!std::isfinite(estimates[j]) || !std::isfinite(variances[j])
            || variances[j] < 0)
            throw std::invalid_argument("finite estimates and nonnegative variances required");
        qbar += estimates[j] / m;
        within += variances[j] / m;
    }
    double between = 0;
    for (double q : estimates) between += (q - qbar) * (q - qbar);
    between /= m - 1;
    return {qbar, within, between, within + (1 + 1/m) * between};
}

int main() {
    const auto result = pool({4, 5, 6}, {0.25, 0.25, 0.25});
    std::cout << std::fixed << std::setprecision(4)
              << "estimate: " << result.estimate << "\n"
              << "within: " << result.within << "\n"
              << "between: " << result.between << "\n"
              << "total variance: " << result.totalVariance << "\n"
              << "standard error: " << std::sqrt(result.totalVariance) << "\n";
}

estimate: 5.0000
within: 0.2500
between: 1.0000
total variance: 1.5833
standard error: 1.2583

References and further reading

• Van Buuren: why and when multiple imputation works — the variance components, finite-imputation correction, and proper imputation.
• mice: pooling analysis estimates — the impute–analyze–pool order and adjusted degrees of freedom.
• mice: Bayesian normal regression imputation — drawing model parameters and residual variation for stochastic normal imputations.
• Van Buuren: parameter pooling — appropriate scales for pooling different statistical quantities.
• Van Buuren: fully conditional specification — model compatibility and monitoring chained imputation.