Multiple Imputation by Chained Equations

One Regression per Variable

Project 1 imputed by fitting one joint model to everything at once, which is elegant while every column is continuous and useless the moment one is binary. Chained equations sidestep the joint model entirely: impute each incomplete variable from its own regression on all the others — a linear model for a continuous column, a probit for a binary one — cycling until the fills stabilise. After a dozen sweeps the imputations behave like draws from a joint distribution that was never written down.

Why multiple, and by how much

The other half is why imputation must be multiple. A single filled-in value is a guess treated as if it were known, and the downstream analysis becomes overconfident by exactly the amount of that pretence. Across repeated simulations, single imputation covers a nominal 95% interval only 77% of the time, while Rubin's rules restore 95%. The missing piece is the between-imputation variance $B$, which single imputation simply omits.

Rubin's combining rules add it back — total variance is within-imputation plus between-imputation, and the ratio gives the fraction of missing information, which on the nhanes regression runs near 0.45 per term.

273 Degrees of Freedom from 25 Observations

The degrees of freedom deserve care, and on small samples — precisely where this project's dataset lives — the choice of formula is not cosmetic. Rubin's original 1987 formula is derived assuming the complete-data degrees of freedom are effectively infinite, which is why Barnard & Rubin (1999) exists. On nhanes — 25 rows, 4 parameters, so 21 complete-data df — it returns between 213 and 273. A pooled analysis cannot have more degrees of freedom than the complete data would supply, so the number is impossible on its face.

Implementing the Barnard–Rubin correction, which combines the old value with an observed-data df, gives 9.6 to 10.7 instead. The point estimates and standard errors are unaffected — but a $t$ interval built on 273 df rather than 10 is far too narrow, and on 25 rows that is the difference between an honest interval and a misleading one. Both columns are printed side by side, because seeing the impossible number is the lesson.

Two engines, and why they differ

The two engines also part company on the coefficients themselves: age 34.4 against 49.2, bmi 5.7 against 7.1, intercept −18.0 against −80.0, and the gap reproduces across seeds. The cause is visible in mice's own methods row — it defaults to predictive mean matching, which imputes by copying an observed value from a similar donor and therefore cannot go outside the observed range, while the from-scratch imputer draws from a Bayesian normal regression. With 25 rows and 27 of 100 cells missing, that choice moves the answer more than the pooling rule does. What genuinely does agree is what the project claims: both use all 25 rows rather than the 13 complete cases, both report a fraction of missing information near 0.4–0.5, and both give the same signs and rough magnitudes. nhanes is a teaching dataset chosen for being small enough to read, not for pinning down coefficients.

Notebooks

Downloads

MICE Module — Source Code

"""
mice_fcs.py -- MULTIPLE IMPUTATION by FULLY CONDITIONAL SPECIFICATION / chained equations (scratch).

Backs the notebooks in  "Multiple Imputation by Chained Equations".

Project 1 imputed by fitting ONE joint model (a multivariate normal) to all variables at once. That
is elegant when everything is continuous, but real datasets mix continuous, binary and categorical
columns, and no single joint distribution fits them all. MULTIPLE IMPUTATION BY CHAINED EQUATIONS
(MICE, van Buuren; also "fully conditional specification", FCS) sidesteps the joint model: it
imputes each incomplete variable from ITS OWN regression on all the others, cycling through the
variables until the fills stabilise. A continuous column gets a linear regression, a binary column a
probit/logistic one -- each variable is imputed by the model that suits it.

Two ideas do the work.

  MULTIPLE IMPUTATION.  A single filled-in value is a guess treated as if it were known, which makes
  the downstream analysis over-confident. Instead we create m > 1 completed datasets, each a
  different plausible draw, analyse every one, and COMBINE the results by Rubin's rules:

      qbar = mean of the m estimates
      T    = Ubar + (1 + 1/m) B        (within-imputation variance + between-imputation variance)

  The between-imputation term B is the extra uncertainty from not knowing the missing values, which
  single imputation simply omits. The fraction of missing information is (1+1/m)B / T.

  CHAINED EQUATIONS.  Within each completed dataset we sweep the variables in turn; for the current
  one we draw a proper Bayesian regression on the others (conjugate linear regression for continuous
  targets, Albert-Chib probit augmentation for binary ones -- the same augmentation used throughout
  this portfolio) and replace its missing entries with a posterior-predictive draw. A dozen sweeps
  and the imputations are a sample from the implied joint distribution, never having specified one.
"""

import numpy as np
from scipy.special import ndtr as Phi, ndtri as Phinv


def _rtrunc_sign(mean, pos, rng):
    lo = np.where(pos, Phi(-mean), 0.0); hi = np.where(pos, 1.0, Phi(-mean))
    u = lo + rng.random(mean.shape) * (hi - lo)
    return mean + Phinv(np.clip(u, 1e-12, 1 - 1e-12))


# --------------------------------------------------------------------------- #
#  per-variable Bayesian imputation models                                      #
# --------------------------------------------------------------------------- #

def _blr_impute(y, X, obs, mis, rng):
    """Bayesian linear regression (conjugate reference prior): fit on observed rows, draw
    (beta, sigma^2) from the posterior, return a posterior-PREDICTIVE draw for the missing rows."""
    Xo = X[obs]; yo = y[obs]; n, p = Xo.shape
    XtXi = np.linalg.inv(Xo.T @ Xo + 1e-8 * np.eye(p))
    bhat = XtXi @ (Xo.T @ yo); resid = yo - Xo @ bhat
    sig2 = (resid @ resid) / rng.chisquare(max(n - p, 1))            # sigma^2 ~ RSS / chi^2_{n-p}
    beta = bhat + np.linalg.cholesky(sig2 * XtXi) @ rng.standard_normal(p)
    return X[mis] @ beta + np.sqrt(sig2) * rng.standard_normal(len(mis))


def _probit_impute(y, X, obs, mis, rng, aug=20):
    """Bayesian probit regression via Albert-Chib augmentation: fit on observed rows, then draw the
    missing entries as Bernoulli(Phi(x'beta)) -- the binary-variable imputer."""
    Xo = X[obs]; yo = y[obs]; n, p = Xo.shape
    XtXi = np.linalg.inv(Xo.T @ Xo + 1e-6 * np.eye(p)); L0 = np.linalg.cholesky(XtXi)
    beta = np.zeros(p)
    for _ in range(aug):
        z = _rtrunc_sign(Xo @ beta, yo > 0.5, rng)
        beta = XtXi @ (Xo.T @ z) + L0 @ rng.standard_normal(p)
    return (rng.random(len(mis)) < Phi(X[mis] @ beta)).astype(float)


# --------------------------------------------------------------------------- #
#  the chained-equations sampler                                                #
# --------------------------------------------------------------------------- #

def fcs_impute(Y, types, rng, m=20, iters=12):
    """Fully conditional specification. Y:(N,p) with NaN; types: list of 'cont' / 'bin' / 'complete'
    per column ('complete' = never imputed, used only as a predictor). Returns a list of m completed
    (N,p) datasets -- the multiple imputations."""
    Y = np.array(Y, float); N, p = Y.shape; Mask = np.isnan(Y)
    completed = []
    for _ in range(m):
        Yc = Y.copy()
        for j in range(p):                                          # initialise from random observed values
            if Mask[:, j].any():
                Yc[Mask[:, j], j] = rng.choice(Yc[~Mask[:, j], j], Mask[:, j].sum())
        for _ in range(iters):                                      # cycle through variables
            for j in range(p):
                if types[j] == "complete" or not Mask[:, j].any():
                    continue
                mis = np.where(Mask[:, j])[0]; obs = np.where(~Mask[:, j])[0]
                others = [k for k in range(p) if k != j]
                X = np.column_stack([np.ones(N)] + [Yc[:, k] for k in others])
                imp = _blr_impute(Yc[:, j], X, obs, mis, rng) if types[j] == "cont" \
                    else _probit_impute(Yc[:, j], X, obs, mis, rng)
                Yc[mis, j] = imp
        completed.append(Yc)
    return completed


# --------------------------------------------------------------------------- #
#  analysis on completed data + Rubin's rules                                   #
# --------------------------------------------------------------------------- #

def ols(y, X):
    """OLS fit -> (coefficients, coefficient covariance matrix)."""
    XtXi = np.linalg.inv(X.T @ X); b = XtXi @ (X.T @ y)
    resid = y - X @ b; s2 = (resid @ resid) / (len(y) - X.shape[1])
    return b, s2 * XtXi


def rubin_pool(betas, covs, dfcom=None):
    """Combine m complete-data fits by Rubin's rules. betas:(m,k), covs:(m,k,k).

    `dfcom` is the COMPLETE-DATA degrees of freedom (n - k). Supply it: Rubin's original (1987)
    df formula assumes dfcom is effectively infinite and can return more degrees of freedom than
    the complete data could ever provide -- on a 25-row dataset it happily reports 190. Barnard &
    Rubin (1999) correct that by combining it with the observed-data df; with dfcom given, `df`
    is the corrected value and `df_rubin` the uncorrected one for comparison."""
    betas = np.asarray(betas); covs = np.asarray(covs); m, k = betas.shape
    qbar = betas.mean(0)
    Ubar = covs.mean(0).diagonal()                                  # within-imputation variance
    B = betas.var(0, ddof=1)                                        # between-imputation variance
    T = Ubar + (1 + 1 / m) * B                                      # total variance
    fmi = (1 + 1 / m) * B / T                                       # fraction of missing information
    with np.errstate(divide="ignore", invalid="ignore"):
        df_rubin = (m - 1) * (1 + Ubar / ((1 + 1 / m) * B)) ** 2     # Rubin (1987), large-sample
        df = df_rubin
        if dfcom is not None:                                        # Barnard & Rubin (1999)
            nu_obs = (dfcom + 1) / (dfcom + 3) * dfcom * (1 - fmi)
            df = 1.0 / (1.0 / df_rubin + 1.0 / nu_obs)
    return dict(estimate=qbar, se=np.sqrt(T), within=Ubar, between=B, fmi=fmi,
                df=df, df_rubin=df_rubin)


def complete_case_ols(Y, ycol, xcols):
    """listwise-deletion regression: keep only fully observed rows across the used columns."""
    keep = ~np.isnan(Y[:, [ycol] + xcols]).any(axis=1)
    X = np.column_stack([np.ones(keep.sum())] + [Y[keep, c] for c in xcols])
    b, cov = ols(Y[keep, ycol], X)
    return b, np.sqrt(np.diag(cov)), keep.sum()

References

• Rubin, D. B. (1987). Multiple Imputation for Nonresponse in Surveys. Wiley. — the combining rules, and the original degrees-of-freedom formula this project had to correct
• Barnard, J. & Rubin, D. B. (1999). Small-sample degrees of freedom with multiple imputation. Biometrika 86(4), 948–955. — the correction, and exactly the situation it was written for
• van Buuren, S. & Groothuis-Oudshoorn, K. (2011). mice: Multivariate imputation by chained equations in R. Journal of Statistical Software 45(3). — the R engine, its pmm default, and the nhanes example
• van Buuren, S. (2018). Flexible Imputation of Missing Data (2nd ed.). CRC Press. — fully conditional specification, and why predictive mean matching behaves differently
• Little, R. J. A. (1988). Missing-data adjustments in large surveys. Journal of Business & Economic Statistics 6(3), 287–296. — predictive mean matching
• Albert, J. H. & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. JASA 88(422), 669–679. — the augmentation used by the binary imputer