Multiple Imputation with mice and miceforest

Single imputation’s flaw is not the guess. It is that the guess is filed as a measurement, so the model prices it at full confidence. Multiple imputation fabricates m plausible data sets instead of one, fits the analysis model m times, and pools the fits with Rubin’s rules: the average within-data-set variance, plus the disagreement between data sets, scaled by 1 + 1/m. The between-imputation disagreement is the price of not knowing, and it belongs in the standard error.

The mice run

The imputation engine below is mice, the reference implementation in R. Give it the cohort with its factor types set and it chooses a method per variable: predictive mean matching for the numeric ones, logistic regression for the binary one, a proportional odds model for the ordered categories. The Python tab builds the same design with miceforest, whose mean-matching step plays the same role: imputed values always come back as observed values, so there are no fractional smokers and no negative cholesterol.

R
Python

suppressPackageStartupMessages(library(mice))

missing_data <- read.csv("/opt/data/cohort_missing.csv", na.strings = "")
complete_data <- read.csv("/opt/data/cohort_complete.csv", na.strings = "")

missing_data$activity <- factor(
  missing_data$activity,
  levels = c("low", "medium", "high"),
  ordered = TRUE
)
missing_data$smoking <- factor(missing_data$smoking, levels = c(0, 1))
complete_data$activity <- factor(
  complete_data$activity,
  levels = c("low", "medium", "high"),
  ordered = TRUE
)
complete_data$smoking <- factor(complete_data$smoking, levels = c(0, 1))

complete_fit <- lm(sbp ~ age + bmi + chol + smoking, data = complete_data)
complete_chol <- coef(summary(complete_fit))["chol", ]
cat(
  sprintf(
    "Complete data: chol coefficient %.4f (SE %.4f)\n",
    complete_chol["Estimate"],
    complete_chol["Std. Error"]
  )
)

cohort_data <- missing_data[, setdiff(names(missing_data), "patient_id")]
imp <- mice(cohort_data, m = 20, seed = 20260820, printFlag = FALSE)
print(imp$method)

Complete data: chol coefficient 3.0370 (SE 0.8299)
     age      bmi activity  smoking     chol      sbp
      ""       ""   "polr" "logreg"    "pmm"    "pmm"

import numpy as np
import pandas as pd
import miceforest
import statsmodels.api as sm

missing_data = pd.read_csv("/opt/data/cohort_missing.csv")
complete_data = pd.read_csv("/opt/data/cohort_complete.csv")
activity_codes = {"low": 1, "medium": 2, "high": 3}

# Mean matching only returns observed values, so these codes remain categories.
missing_data["activity"] = missing_data["activity"].map(activity_codes)
complete_data["activity"] = complete_data["activity"].map(activity_codes)
frame = missing_data.drop(columns="patient_id")

complete_design = sm.add_constant(
    complete_data[["age", "bmi", "chol", "smoking"]],
)
complete_fit = sm.OLS(complete_data["sbp"], complete_design).fit()
print(
    "Complete data: chol coefficient "
    f"{complete_fit.params['chol']:.4f} "
    f"(SE {complete_fit.bse['chol']:.4f})"
)

# mean_match_candidates is the predictive mean matching analogue: each imputed
# value is drawn from the observed values most likely to be the true one.
kernel = miceforest.ImputationKernel(
    frame,
    num_datasets=20,
    random_state=20260820,
    mean_match_candidates=5,
    save_all_iterations_data=False,
)
# One thread keeps the lightgbm fit deterministic across gate re-runs.
kernel.mice(iterations=5, verbose=False, num_threads=1)
imputed_cells = int(frame.isna().sum().sum())
print(f"Cells imputed per dataset: {imputed_cells}")

Complete data: chol coefficient 3.0370 (SE 0.8299)
Cells imputed per dataset: 210

The benchmark stands at 3.0370 with a standard error of 0.8299, and the truth target of this page is that pair. The R tab also reports the method mice chose: nothing for the two complete variables, pmm for cholesterol and blood pressure, logreg for smoking, polr for activity. Choosing the model per variable type is mice’s job one; the next page takes it apart.

Pooling with Rubin’s rules

Fit the model on each of the 20 completed data sets and pool. Rubin’s rules sum the mean within-imputation variance and the between-imputation variance with a finite-m correction. The pooled standard error is then honest: it includes a term for the uncertainty that the imputed values introduced. mice reports a t interval on Barnard-Rubin degrees of freedom; the Python tab here pools by hand with the same variance decomposition and a normal approximation, a difference flagged in the code comments.

R
Python

fit <- with(imp, lm(sbp ~ age + bmi + chol + smoking))
pooled <- pool(fit)
print(summary(pooled, conf.int = TRUE))

         term   estimate  std.error statistic       df      p.value      2.5 %
1 (Intercept) 61.2131164 5.02494766 12.181842 197.4217 7.998547e-26 51.3036533
2         age  0.9265704 0.05987685 15.474601 216.6339 7.071040e-37  0.8085546
3         bmi  0.6917429 0.17841997  3.877049 263.6103 1.335509e-04  0.3404333
4        chol  2.8645593 0.92607684  3.093220 171.3441 2.312390e-03  1.0365709
5    smoking1  7.9024109 1.45861066  5.417766 135.3542 2.677274e-07  5.0177961
     97.5 %   conf.low conf.high
1 71.122579 51.3036533 71.122579
2  1.044586  0.8085546  1.044586
3  1.043053  0.3404333  1.043053
4  4.692548  1.0365709  4.692548
5 10.787026  5.0177961 10.787026

def pool_rubins(estimates, standard_errors):
    """Pool estimates with Rubin's rules.

    Args:
        estimates: One estimate from each imputed data set.
        standard_errors: One standard error from each imputed data set.

    Returns:
        A tuple with the pooled estimate, standard error, lower bound, and
        upper bound of a normal-approximation 95 percent confidence interval.
    """
    estimate_array = np.asarray(estimates, dtype=float)
    error_array = np.asarray(standard_errors, dtype=float)
    imputation_count = len(estimate_array)
    pooled_estimate = float(np.mean(estimate_array))
    within_variance = float(np.mean(error_array**2))
    between_variance = float(np.var(estimate_array, ddof=1))
    total_variance = within_variance + (
        1 + 1 / imputation_count
    ) * between_variance
    pooled_error = float(np.sqrt(total_variance))
    lower_bound = pooled_estimate - 1.96 * pooled_error
    upper_bound = pooled_estimate + 1.96 * pooled_error
    return pooled_estimate, pooled_error, lower_bound, upper_bound

def format_float(value):
    """Format a floating-point value for the tutorial table.

    Args:
        value: The value to format.

    Returns:
        The value with four decimal places.
    """
    return f"{value:.4f}"

coefficient_names = ["const", "age", "bmi", "chol", "smoking"]
estimates_by_term = {name: [] for name in coefficient_names}
errors_by_term = {name: [] for name in coefficient_names}

for dataset_index in range(20):
    completed_data = kernel.complete_data(dataset=dataset_index)
    design = sm.add_constant(
        completed_data[["age", "bmi", "chol", "smoking"]],
    )
    model_fit = sm.OLS(completed_data["sbp"], design).fit()
    for term in coefficient_names:
        estimates_by_term[term].append(model_fit.params[term])
        errors_by_term[term].append(model_fit.bse[term])

# mice reports a t interval with Barnard-Rubin degrees of freedom, so the R tab
# intervals differ slightly in width.
pooled_rows = {}
for term in coefficient_names:
    pooled_estimate, pooled_error, lower_bound, upper_bound = pool_rubins(
        estimates_by_term[term],
        errors_by_term[term],
    )
    pooled_rows[term] = {
        "estimate": pooled_estimate,
        "std.error": pooled_error,
        "2.5 %": lower_bound,
        "97.5 %": upper_bound,
    }

pooled_table = pd.DataFrame.from_dict(pooled_rows, orient="index")
print(pooled_table.to_string(float_format=format_float))

         estimate  std.error   2.5 %  97.5 %
const     61.9191     4.4011 53.2929 70.5453
age        0.8803     0.0541  0.7743  0.9862
bmi        0.5616     0.1637  0.2408  0.8825
chol       3.7856     0.8254  2.1678  5.4034
smoking    8.5084     1.2181  6.1210 10.8958

mice pools the cholesterol term at 2.8646 with a standard error of 0.9261 and a Barnard-Rubin 171 degrees of freedom, and the confidence interval from 1.0366 to 4.6925 contains the benchmark. miceforest pools at 3.7856 with a standard error of 0.8254, and its interval from 2.1678 to 5.4034 contains the benchmark too. The point estimates bracket the truth from opposite sides, which is exactly the spread you pay for when the imputation model is a choice and not a fact.

The disagreement is the information

The reason a pooled standard error can be honest is that the m fits tell you how loud the missing data speak. Collect the cholesterol coefficient from each of the 20 fits and look at the spread.

R
Python

chol_estimates <- sapply(fit$analyses, coef)["chol", ]
cat(sprintf("Chol coefficient minimum: %.4f\n", min(chol_estimates)))
cat(sprintf("Chol coefficient maximum: %.4f\n", max(chol_estimates)))
cat(sprintf("Chol coefficient standard deviation: %.4f\n", sd(chol_estimates)))

Chol coefficient minimum: 1.9353
Chol coefficient maximum: 3.6176
Chol coefficient standard deviation: 0.4248

chol_estimates = estimates_by_term["chol"]
print(f"Chol coefficient minimum: {min(chol_estimates):.4f}")
print(f"Chol coefficient maximum: {max(chol_estimates):.4f}")
print(
    "Chol coefficient standard deviation: "
    f"{np.std(chol_estimates, ddof=1):.4f}"
)

Chol coefficient minimum: 3.4040
Chol coefficient maximum: 4.2102
Chol coefficient standard deviation: 0.2049

mice’s twenty fits spread the cholesterol term from 1.9353 to 3.6176, a standard deviation of 0.4248 that feeds the 0.9261 standard error. miceforest’s twenty fits span 3.4040 to 4.2102 with a standard deviation of 0.2049, roughly half, so its pooled standard error lands close to the single-imputation level. A boosted mean-matching engine is more certain about its guesses than a parametric chain is, and this page cannot tell you which certainty is right. The evaluation page puts the question to the truth, which is the only referee the cohort has.

Single versus multiple

One more measurement before leaving single imputation behind, because it is the argument for everything above.

# A single fabricated data set is priced by lm as if it were measured.
single_imp <- mice(
  cohort_data,
  m = 1,
  maxit = 5,
  seed = 20260821,
  printFlag = FALSE
)
single_fit <- lm(
  sbp ~ age + bmi + chol + smoking,
  data = complete(single_imp)
)
single_chol <- coef(summary(single_fit))["chol", ]
cat(
  sprintf(
    "Single pmm imputation: chol coefficient %.4f (SE %.4f)\n",
    single_chol["Estimate"],
    single_chol["Std. Error"]
  )
)

Single pmm imputation: chol coefficient 2.6593 (SE 0.8142)

A single predictive-mean-matching imputation returns 2.6593 with a standard error of 0.8142, practically the benchmark’s 0.8299. The model cannot see which rows were fabricated, so it prices the fabricated ones as measurements. That is the mean-imputation disease from the previous page in its most polite form: the best single guess still lies about confidence, and only the pooled run prices the guesswork.

How many imputations

The cost of m is compute, so the only question is how large m has to be before the pooled answer stops moving.

imp_5 <- mice(cohort_data, m = 5, seed = 20260820, printFlag = FALSE)
imp_100 <- mice(cohort_data, m = 100, seed = 20260822, printFlag = FALSE)

for (imputation_count in c(5, 20, 100)) {
  current_imp <- switch(
    as.character(imputation_count),
    "5" = imp_5,
    "20" = imp,
    "100" = imp_100
  )
  current_fit <- with(
    current_imp,
    lm(sbp ~ age + bmi + chol + smoking)
  )
  current_pooled <- pool(current_fit)
  current_summary <- summary(current_pooled)
  current_chol <- current_summary[current_summary$term == "chol", ]
  cat(
    sprintf(
      "m = %d: pooled chol coefficient %.4f (SE %.4f)\n",
      imputation_count,
      current_chol[["estimate"]],
      current_chol[["std.error"]]
    )
  )
}

m = 5: pooled chol coefficient 2.9162 (SE 0.8420)
m = 20: pooled chol coefficient 2.8646 (SE 0.9261)
m = 100: pooled chol coefficient 2.7851 (SE 0.9536)

At m = 5 the standard error is 0.8420 and the estimate is 2.9162; at m = 100 the standard error has climbed to 0.9536 with the estimate at 2.7851, and m = 20 sits between them. The point estimate is stable from m = 20 on, but the standard error itself carries Monte Carlo noise, and at m = 5 that noise is larger than the margin between the honest 0.93 and the fake 0.83. The common guidance is to set m at least to the percentage of incomplete information, and here the largest per-variable share is 15.5 percent missing in cholesterol, so m = 20 is comfortable and m = 100 only steadies the standard error.

Single imputation is outpriced; multiple imputation is the default. The pages left are about what model fabricates the values and how to check its work.