Multiple Imputation by Chained Equations (MICE)¶
One regression per variable, m completed datasets, Rubin's rules¶
Project 1 imputed by fitting one joint model — a multivariate normal — to every variable at once. That is clean when all variables are continuous, but real data 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) is the practical workhorse that sidesteps the joint model. It rests on two ideas.
Multiple imputation. A single filled-in value is a guess treated as if it were known, which makes the downstream analysis over-confident. So we build $m>1$ completed datasets, analyse each, and combine by Rubin's rules: $$\bar q=\frac1m\sum_i \hat q_i,\qquad T=\bar U+\Big(1+\tfrac1m\Big)B,$$ the total variance being the average within-imputation variance $\bar U$ plus the between-imputation variance $B$ — the extra uncertainty of not knowing the missing values, which single imputation simply drops. The fraction of missing information is $(1+1/m)B/T$.
Chained equations. Rather than specify a joint model, impute each incomplete variable from its own regression on all the others, cycling through the variables until the fills stabilise — a linear regression for a continuous column, a probit/logistic one for a binary column. Each variable gets the model that suits it. We build the sampler from scratch (conjugate linear regression + Albert–Chib probit imputers), prove Rubin's rules restore correct coverage, run the canonical nhanes data, and cross-check against a single joint PyMC model.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import mice_fcs as F
rng = np.random.default_rng(2)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("MICE imputes each variable by its own regression on the others (mixed types welcome), builds m completed")
print("datasets, and pools the analyses by Rubin's rules so the missing-data uncertainty is carried, not lost.")
MICE imputes each variable by its own regression on the others (mixed types welcome), builds m completed datasets, and pools the analyses by Rubin's rules so the missing-data uncertainty is carried, not lost.
1. Why multiple imputation — single imputation lies about uncertainty¶
Fill a gap with one plausible value and analyse the result, and you have told the analysis that value was observed. The estimate may be fine, but its standard error is too small — you get false confidence. Multiple imputation fixes this by drawing the missing values several times: the spread across completed datasets is exactly the uncertainty a single fill throws away. The picture below shows one coefficient estimated from several completed datasets — each imputation gives a slightly different answer, and that scatter is real information.
# one MAR dataset, show the estimate wobbling across imputations
N=200; x1=rng.standard_normal(N); x2=0.6*x1+0.8*rng.standard_normal(N); y=1+0.5*x1+rng.standard_normal(N)
Y=np.column_stack([y,x1,x2]); pm=1/(1+np.exp(-x2)); Y[rng.random(N)<pm*0.45/pm.mean(),1]=np.nan
comp=F.fcs_impute(Y,["complete","cont","cont"],rng,m=25,iters=10)
slopes=[F.ols(d[:,0],np.column_stack([np.ones(N),d[:,1],d[:,2]]))[0][1] for d in comp]
fig,ax=plt.subplots(figsize=(8,3))
ax.scatter(slopes,np.zeros_like(slopes),color=BLUE,alpha=.6,s=50)
ax.axvline(np.mean(slopes),color=RED,lw=2,label=f"pooled mean {np.mean(slopes):.3f}"); ax.axvline(0.5,color=GREEN,lw=2,ls="--",label="truth 0.5")
ax.set_yticks([]); ax.set_xlabel(r"estimated slope on $x_1$ (one point per imputation)"); ax.set_title("Each imputation gives a different estimate — that scatter is the missing-data uncertainty")
ax.legend(frameon=False,fontsize=8); plt.tight_layout(); plt.show()
print(f"The {len(slopes)} imputations scatter around the truth; their spread (between-imputation variance) is precisely")
print("what single imputation ignores. Rubin's rules add it back into the standard error.")
The 25 imputations scatter around the truth; their spread (between-imputation variance) is precisely what single imputation ignores. Rubin's rules add it back into the standard error.
2. Rubin's rules restore correct coverage¶
The claim to verify: multiple imputation with Rubin's rules gives confidence intervals that cover at their nominal rate, while single imputation's intervals are too narrow. We repeat the whole exercise over many simulated MAR datasets and record how often the 95% interval for the slope actually contains the truth.
def one_rep(seed):
r=np.random.default_rng(seed); n=200
a=r.standard_normal(n); b=0.6*a+0.8*r.standard_normal(n); yy=1+0.5*a+r.standard_normal(n)
D=np.column_stack([yy,a,b]); p=1/(1+np.exp(-b)); D[r.random(n)<p*0.45/p.mean(),1]=np.nan
cc=F.fcs_impute(D,["complete","cont","cont"],r,m=20,iters=10)
bb=[]; vv=[]
for d in cc:
X=np.column_stack([np.ones(n),d[:,1],d[:,2]]); be,co=F.ols(d[:,0],X); bb.append(be); vv.append(co)
pr=F.rubin_pool(np.array(bb),np.array(vv)); est,se=pr["estimate"][1],pr["se"][1]
d0=cc[0]; X=np.column_stack([np.ones(n),d0[:,1],d0[:,2]]); b1,c1=F.ols(d0[:,0],X)
return abs(est-0.5)<1.96*se, abs(b1[1]-0.5)<1.96*np.sqrt(c1[1,1]), pr["fmi"][1]
R=[one_rep(s) for s in range(60)]
mi_cov=np.mean([r[0] for r in R]); si_cov=np.mean([r[1] for r in R]); fmi=np.mean([r[2] for r in R])
fig,ax=plt.subplots(figsize=(6.5,3.2))
ax.bar(["multiple\nimputation","single\nimputation"],[mi_cov,si_cov],color=[GREEN,RED],alpha=.85)
ax.axhline(0.95,color="k",ls="--",lw=1,label="nominal 95%"); ax.set_ylim(0,1); ax.set_ylabel("actual 95%-CI coverage")
ax.set_title("Rubin's rules cover correctly; single imputation is over-confident"); ax.legend(frameon=False,fontsize=8)
for i,v in enumerate([mi_cov,si_cov]): ax.text(i,v+0.02,f"{v:.0%}",ha="center")
plt.tight_layout(); plt.show()
print(f"multiple imputation covers at {mi_cov:.0%} (nominal 95%); single imputation only {si_cov:.0%} -- it understates the")
print(f"standard error by omitting the between-imputation variance. Mean fraction of missing information: {fmi:.2f}.")
multiple imputation covers at 95% (nominal 95%); single imputation only 77% -- it understates the standard error by omitting the between-imputation variance. Mean fraction of missing information: 0.39.
3. Chained equations with mixed variable types¶
The advantage over Project 1's joint normal: each variable is imputed by the model that fits it. We simulate three variables — two continuous and one binary — induce missingness in all three, and impute by cycling: linear regression for the continuous columns, Albert–Chib probit for the binary one. The recovered means and the binary column's rate should match the truth.
n=1500
z1=rng.standard_normal(n); z2=0.5*z1+rng.standard_normal(n)
pb=1/(1+np.exp(-(0.8*z1))); z3=(rng.random(n)<pb).astype(float) # binary, depends on z1
D=np.column_stack([z1,z2,z3])
for j,fr in [(0,0.25),(1,0.3),(2,0.25)]:
D[rng.random(n)<fr,j]=np.nan
comp=F.fcs_impute(D,["cont","cont","bin"],rng,m=15,iters=15)
imp_mean=np.mean([c.mean(0) for c in comp],axis=0)
print("means -- true [0, 0, %.2f] ; MICE imputed %s" % (pb.mean(), np.round(imp_mean,2)))
print("binary column: true rate %.2f, imputed rate %.2f (probit imputer keeps it 0/1 and calibrated)" % (z3.mean(), imp_mean[2]))
# show imputed binary distribution matches
fig,ax=plt.subplots(1,2,figsize=(11,3.4))
allc=np.array(comp); mis2=np.isnan(D[:,2])
ax[0].bar(["0","1"],[np.mean(1-allc[:,mis2,2]),np.mean(allc[:,mis2,2])],color=PURP,alpha=.8)
ax[0].set_title("imputed values of the binary column"); ax[0].set_ylabel("proportion")
ax[1].hist(D[~np.isnan(D[:,0]),0],bins=30,density=True,color=BLUE,alpha=.6,label="observed z1")
ax[1].hist(allc[:,np.isnan(D[:,0]),0].ravel(),bins=30,density=True,color=RED,alpha=.5,label="imputed z1")
ax[1].legend(frameon=False,fontsize=8); ax[1].set_title("continuous column: imputed matches observed"); ax[1].set_xlabel("z1")
plt.tight_layout(); plt.show()
print("Mixed types are handled seamlessly -- something the single joint-normal model of Project 1 cannot do. Each")
print("variable is imputed by its own regression, which is the whole idea of fully conditional specification.")
means -- true [0, 0, 0.51] ; MICE imputed [ 0.03 -0.03 0.5 ] binary column: true rate 0.49, imputed rate 0.50 (probit imputer keeps it 0/1 and calibrated)
Mixed types are handled seamlessly -- something the single joint-normal model of Project 1 cannot do. Each variable is imputed by its own regression, which is the whole idea of fully conditional specification.
4. Real data — nhanes¶
The nhanes set (25 people; age complete, bmi, hypertension and chlolesterol missing) is the canonical MICE example — and a cautionary one: only 13 rows are complete, so listwise deletion would nearly halve the data. We impute the three incomplete variables (bmi and chl continuous, hyp binary) and fit the regression chl ~ age + bmi + hyp, pooling by Rubin's rules, then compare with the complete-case fit.
d=pd.read_csv("nhanes.csv"); Y=d[["age","bmi","hyp","chl"]].to_numpy(float)
comp=F.fcs_impute(Y,["complete","cont","bin","cont"],rng,m=50,iters=15)
betas=[]; covs=[]
for c in comp:
X=np.column_stack([np.ones(25),c[:,0],c[:,1],c[:,2]]); b,cov=F.ols(c[:,3],X); betas.append(b); covs.append(cov)
# dfcom = complete-data degrees of freedom, n - k. Supplying it switches rubin_pool to the
# Barnard-Rubin (1999) small-sample correction; without it you get Rubin's 1987 formula,
# which on 25 rows reports well over a hundred degrees of freedom.
pr=F.rubin_pool(np.array(betas),np.array(covs),dfcom=25-4); names=["intercept","age","bmi","hyp"]
cc_b,cc_se,ncc=F.complete_case_ols(Y,3,[0,1,2])
print(f"chl ~ age + bmi + hyp (MICE m=50 uses all 25 rows; complete-case uses only {ncc})\n")
print(f"{'term':10s}{'MICE est':>10s}{'MICE SE':>9s}{'fmi':>6s}{'df (BR99)':>11s}{'df (R87)':>10s} {'complete-case est':>18s}{'SE':>7s}")
for i,nm in enumerate(names):
print(f"{nm:10s}{pr['estimate'][i]:10.2f}{pr['se'][i]:9.2f}{pr['fmi'][i]:6.2f}{pr['df'][i]:11.1f}{pr['df_rubin'][i]:10.1f} {cc_b[i]:18.2f}{cc_se[i]:7.2f}")
fig,ax=plt.subplots(figsize=(7,3))
idx=[1,2,3]; ax.errorbar(pr["estimate"][idx],np.arange(3)+.1,xerr=1.96*pr["se"][idx],fmt="o",color=GREEN,capsize=3,label="MICE (n=25)")
ax.errorbar(cc_b[idx],np.arange(3)-.1,xerr=1.96*cc_se[idx],fmt="s",color=RED,capsize=3,label="complete-case (n=13)")
ax.axvline(0,color="k",lw=.8); ax.set_yticks(range(3)); ax.set_yticklabels([names[i] for i in idx]); ax.set_xlabel("coefficient")
ax.set_title("chl regression: MICE keeps all the data"); ax.legend(frameon=False,fontsize=8); plt.tight_layout(); plt.show()
print("\nLook at the two df columns: BR99 is the Barnard-Rubin (1999) small-sample correction,")
print("R87 is Rubin's 1987 formula. The complete data would supply n - k = %d degrees of freedom, so" % (25-4))
print("any pooled df above that is impossible -- and Rubin's 1987 formula returns %.0f to %.0f, because it"
% (pr["df_rubin"].min(), pr["df_rubin"].max()))
print("is derived assuming the complete-data df is effectively infinite. Barnard and Rubin (1999) fix")
print("this by combining it with the observed-data df, giving %.1f to %.1f here. On 25 rows that is the"
% (pr["df"].min(), pr["df"].max()))
print("difference between a t interval that is roughly right and one that is far too narrow, so it is")
print("worth getting right even though the point estimates and standard errors are unaffected.")
print("MICE recovers similar coefficients but uses every row and reports an honest fraction of missing information")
print("per term; complete-case throws away 12 of 25 people. With larger missingness or MAR bias the two would")
print("diverge -- here the gain is efficiency and calibrated uncertainty.")
chl ~ age + bmi + hyp (MICE m=50 uses all 25 rows; complete-case uses only 13) term MICE est MICE SE fmi df (BR99) df (R87) complete-case est SE intercept -79.97 71.90 0.48 9.6 213.0 -87.19 66.99 age 49.18 13.65 0.46 9.9 228.5 55.21 14.29 bmi 7.14 2.19 0.47 9.7 221.3 7.07 2.05 hyp -6.39 23.08 0.42 10.7 272.8 -6.22 23.18
Look at the two df columns: BR99 is the Barnard-Rubin (1999) small-sample correction, R87 is Rubin's 1987 formula. The complete data would supply n - k = 21 degrees of freedom, so any pooled df above that is impossible -- and Rubin's 1987 formula returns 213 to 273, because it is derived assuming the complete-data df is effectively infinite. Barnard and Rubin (1999) fix this by combining it with the observed-data df, giving 9.6 to 10.7 here. On 25 rows that is the difference between a t interval that is roughly right and one that is far too narrow, so it is worth getting right even though the point estimates and standard errors are unaffected. MICE recovers similar coefficients but uses every row and reports an honest fraction of missing information per term; complete-case throws away 12 of 25 people. With larger missingness or MAR bias the two would diverge -- here the gain is efficiency and calibrated uncertainty.
5. Cross-check — one joint model in PyMC¶
MICE is modular: separate regressions, stitched together by Rubin's rules. The Bayesian ideal it approximates is a single joint model that imputes and estimates at once — the MCMC integrates over the missing values automatically, so no pooling step is needed. The two should agree when the data are informative enough that priors do not dominate (unlike nhanes' 25 rows). We simulate a larger MAR dataset, recover the slope both ways, and compare — MICE's pooled estimate against the PyMC joint model's coefficient posterior.
import pymc as pm
n=400; xa=rng.standard_normal(n); xb=0.6*xa+0.8*rng.standard_normal(n); yv=1.0+0.8*xa+0.5*xb+rng.standard_normal(n)
Dsim=np.column_stack([yv,xa,xb]); pmar=1/(1+np.exp(-xb)); Dsim[rng.random(n)<pmar*0.4/pmar.mean(),1]=np.nan # xa MAR on xb
# MICE + Rubin's rules
cs=F.fcs_impute(Dsim,["complete","cont","cont"],rng,m=30,iters=12)
bb=[]; cc=[]
for d in cs:
X=np.column_stack([np.ones(n),d[:,1],d[:,2]]); b,c=F.ols(d[:,0],X); bb.append(b); cc.append(c)
prm=F.rubin_pool(np.array(bb),np.array(cc))
# PyMC single joint model: model xa (imputes it) + regression for y
xa_m=np.ma.masked_invalid(Dsim[:,1]); xb_o=Dsim[:,2]; y_o=Dsim[:,0]
with pm.Model() as mod:
g0=pm.Normal("g0",0,5); g1=pm.Normal("g1",0,5); sx=pm.HalfNormal("sx",5)
xa_=pm.Normal("xa", g0+g1*xb_o, sx, observed=xa_m) # covariate model -> imputes missing xa
a0=pm.Normal("a0",0,5); b1=pm.Normal("b1",0,5); b2=pm.Normal("b2",0,5); s=pm.HalfNormal("s",5)
pm.Normal("y", a0+b1*xa_+b2*xb_o, s, observed=y_o)
idata=pm.sample(1000,tune=1500,chains=4,target_accept=0.9,random_seed=4,progressbar=False)
bp=idata.posterior["b1"]
print("slope on the MAR covariate (truth 0.80):")
print(f" MICE (Rubin pooled): {prm['estimate'][1]:.3f} +/- {prm['se'][1]:.3f}")
print(f" PyMC (joint model) : {float(bp.mean()):.3f} +/- {float(bp.std()):.3f}")
print("With informative data the modular MICE pooling and the single joint Bayesian model agree closely: MICE is a")
print("practical, mixed-type-friendly approximation to the one joint model PyMC fits directly. (On nhanes' 25 rows")
print("the two would differ, because there the priors -- not the data -- decide, which is a small-sample caveat, not")
print("a flaw in either method.)")
g++ not available, if using conda: `conda install gxx`
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pymc\model\core.py:1337: ImputationWarning: Data in xa contains missing values and will be automatically imputed from the sampling distribution. warnings.warn(impute_message, ImputationWarning)
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [g0, g1, sx, xa_unobserved, a0, b1, b2, s]
Sampling 4 chains for 1_500 tune and 1_000 draw iterations (6_000 + 4_000 draws total) took 5 seconds.
slope on the MAR covariate (truth 0.80): MICE (Rubin pooled): 0.843 +/- 0.079 PyMC (joint model) : 0.839 +/- 0.079 With informative data the modular MICE pooling and the single joint Bayesian model agree closely: MICE is a practical, mixed-type-friendly approximation to the one joint model PyMC fits directly. (On nhanes' 25 rows the two would differ, because there the priors -- not the data -- decide, which is a small-sample caveat, not a flaw in either method.)
6. Summary¶
Multiple imputation by chained equations is the practical face of Bayesian missing-data handling. Where Project 1 fit one joint normal, MICE imputes each variable from its own regression on the others — linear for continuous, probit for binary — cycling until stable, so it copes with the mixed variable types real data always have. Building $m$ completed datasets and pooling the analyses by Rubin's rules carries the missing-data uncertainty forward: the simulation confirmed that multiple imputation covers at its nominal 95% while single imputation, ignoring the between-imputation variance, sinks to ~77%. On nhanes it used all 25 people where listwise deletion kept only 13, reporting a fraction of missing information per coefficient. A single joint PyMC model reproduced the pooled estimate — MICE is a modular approximation to that one Bayesian model.
The connections: each conditional imputer is a Bayesian regression from this portfolio's toolkit, and the binary one is the Albert–Chib probit augmentation seen in the IRT, Tobit and LCA projects; the completed datasets are the same objects Project 1's data augmentation produced, now combined rather than just stored. MICE assumes MAR — each variable's missingness explained by the others. When missingness depends on the unobserved value itself (MNAR), no chained equation can fix it; that is the subject of Project 4 (selection models) and Project 5 (pattern-mixture). Next, though, Project 3 zooms in on the most common case — missing covariates in a regression — and the choice between imputing them and modelling them jointly.