Missing Data — Foundations¶
Mechanisms, ignorability, and imputation as data augmentation¶
Real datasets have holes. What you do about them decides whether your analysis is honest. This first project of the missing-data arc lays the groundwork: why data go missing (Rubin's taxonomy), when you can ignore the reason, and how the Bayesian answer — treat each missing value as an unknown and sample it — falls straight out of machinery this portfolio already uses.
Rubin's mechanisms. Let $R$ mark which entries are observed. Missingness is MCAR (completely at random) if $R$ is independent of all data; MAR (at random) if $R$ depends only on the observed values; MNAR (not at random) if it depends on the values that are missing. Under MCAR and MAR the mechanism is ignorable — we can impute from the data model alone without modelling why data are missing. Complete-case (listwise) deletion, the default in most software, is unbiased only under MCAR and always throws away information.
The Bayesian idea. A missing value $y_{\text{mis}}$ is just another unknown. Give the data a model $y\sim N(\mu,\Sigma)$ and the posterior over $(\mu,\Sigma,y_{\text{mis}})$ is sampled by a two-step Gibbs — data augmentation (Tanner–Wong; Schafer): $$\textbf{I-step: } y_{\text{mis}}\mid y_{\text{obs}},\mu,\Sigma \sim N\!\big(\mu_{\text{mis}}+B(y_{\text{obs}}-\mu_{\text{obs}}),\,\Sigma_{\text{mis}\mid\text{obs}}\big),\qquad \textbf{P-step: } \mu,\Sigma\mid y_{\text{complete}}\sim \text{NIW}.$$ This is the same "impute the unknown, then update the parameters" loop behind the latent $z$ of Albert–Chib probit, the censored draws of Tobit, and the latent class of LCA. Missing data is the general case. We validate the sampler, show how complete-case deletion goes wrong under MAR, impute the classic airquality data, flag the MNAR limit that later projects address, and cross-check in PyMC.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import missdata as M
rng = np.random.default_rng(1)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("A missing value is an unknown: give it a model and sample it. Under MCAR/MAR the mechanism is ignorable;")
print("complete-case deletion is safe only under MCAR. The engine is data augmentation -- impute, update, repeat.")
A missing value is an unknown: give it a model and sample it. Under MCAR/MAR the mechanism is ignorable; complete-case deletion is safe only under MCAR. The engine is data augmentation -- impute, update, repeat.
1. The three mechanisms¶
A quick simulation makes the taxonomy concrete. Two correlated variables; we knock out entries of the first three ways — MCAR (at random), MAR (missing when the observed second variable is high), MNAR (missing when the first variable's own value is high) — and look at what the surviving values of the first variable look like. Under MCAR the observed subsample is representative; under MAR and MNAR it is skewed, and that skew is what biases naive analyses.
mu=np.array([0.,0.]); S=np.array([[1,.75],[.75,1]]); Y=M.simulate_mvn(6000,mu,S,rng)
mechs={"MCAR":M.induce_missing(Y,rng,"MCAR",0.5,0),
"MAR (on observed X2)":M.induce_missing(Y,rng,"MAR",0.5,target=0,driver=1),
"MNAR (on own X1)":M.induce_missing(Y,rng,"MNAR",0.5,target=0)}
fig,ax=plt.subplots(1,3,figsize=(13,3.6),sharex=True,sharey=True)
for k,(name,Ym) in enumerate(mechs.items()):
obs=~np.isnan(Ym[:,0])
ax[k].scatter(Y[~obs,1],Y[~obs,0],s=6,color="#d0d0d0",label="missing X1")
ax[k].scatter(Ym[obs,1],Ym[obs,0],s=6,color=BLUE,label="observed X1")
ax[k].axhline(Y[:,0].mean(),color=GREEN,lw=1.5,ls="--"); ax[k].axhline(np.nanmean(Ym[:,0]),color=RED,lw=1.5)
ax[k].set_title(f"{name}\nobs. mean X1 = {np.nanmean(Ym[:,0]):+.2f}",fontsize=9); ax[k].set_xlabel("X2")
ax[0].set_ylabel("X1"); ax[0].legend(frameon=False,fontsize=7,loc="lower right")
plt.suptitle("green dashed = true mean (0); red = mean of the surviving X1 values",y=1.06,fontsize=9); plt.tight_layout(); plt.show()
print("MCAR: the observed X1 still averages ~0 -- deletion is safe. MAR/MNAR: the survivors are a biased slice")
print("(red shifted off green), so throwing away the incomplete rows corrupts the estimate. That is the problem.")
MCAR: the observed X1 still averages ~0 -- deletion is safe. MAR/MNAR: the survivors are a biased slice (red shifted off green), so throwing away the incomplete rows corrupts the estimate. That is the problem.
2. Imputation as data augmentation, step by step¶
Before running the sampler, it is worth being explicit about what it does, because the results below are labelled with three competing estimators — complete-case, EM, and DA — and the difference between them is the whole point.
A missing value is simply an unknown, so we give it a distribution and sample it alongside the parameters. Under the multivariate normal model that is data augmentation (Tanner–Wong; Schafer) — a two-step Gibbs loop repeated until it converges:
- Impute step (I). For every missing entry, draw a value from its conditional distribution given that row's observed entries and the current $(\mu,\Sigma)$: $y_{\text{mis}}\mid y_{\text{obs}}\sim N\!\big(\mu_{\text{mis}}+B(y_{\text{obs}}-\mu_{\text{obs}}),\,\Sigma_{\text{mis}\mid\text{obs}}\big)$ with $B=\Sigma_{\text{mis},\text{obs}}\Sigma_{\text{obs},\text{obs}}^{-1}$. This "augments" the incomplete data into a complete dataset — and crucially it borrows strength from the observed variables through their correlation, rather than guessing blindly.
- Parameter step (P). With the data now complete, draw fresh $(\mu,\Sigma)$ from the conjugate normal-inverse-Wishart posterior.
Iterating produces a joint posterior over both the parameters and the missing values. So an estimate labelled "DA imputed" is one computed after filling the gaps this way — reported as a posterior mean with a credible interval that honestly reflects the uncertainty of never having measured those entries. The two rivals it is compared against are:
| label | what it does | valid when |
|---|---|---|
| complete-case | delete every row with any missing entry, analyse the rest | MCAR only |
| EM | frequentist maximum likelihood — expectation replaces the I-step draw, maximisation the P-step draw | MCAR/MAR (point estimate, no built-in uncertainty) |
| DA | the Bayesian data-augmentation posterior over parameters and missing values | MCAR/MAR (estimate and honest uncertainty) |
This is the same impute-then-update loop as the latent $z$ of Albert–Chib probit and the censored draw of Tobit — here the "latent" quantity is just the missing data itself. The stored completed datasets it produces are, as we will see, exactly the multiple imputations that feed Rubin's rules in Project 2.
3. Does the sampler work? — recovering $\mu,\Sigma$¶
Under MCAR, imputing should simply recover the truth with the right uncertainty. We knock out 35% of one column of a three-variate normal and check that the data-augmentation posterior recovers $\mu$ and $\Sigma$, matching the frequentist EM maximum likelihood.
mu=np.array([0.,2.,-1.]); S=np.array([[1,.7,.4],[.7,1.5,.3],[.4,.3,2.]])
Y=M.simulate_mvn(4000,mu,S,rng); Ym=M.induce_missing(Y,rng,"MCAR",0.35,target=0)
res=M.da_gibbs(Ym,rng,draws=2000,burn=1000); em_mu,em_S=M.em_mvn(Ym)
print("mu : true",mu," DA",res["mu"].mean(0).round(2)," EM(MLE)",em_mu.round(2))
print("diag Sigma: true",np.diag(S)," DA",np.diag(res["Sigma"].mean(0)).round(2))
fig,ax=plt.subplots(1,2,figsize=(11,3.8))
ax[0].plot(res["mu"][:,0],color=BLUE,lw=.6); ax[0].axhline(mu[0],color=GREEN,ls="--")
ax[0].set_title(r"trace of $\mu_1$ (imputed column)"); ax[0].set_xlabel("draw")
lo,hi=np.percentile(res["mu"][:,0],[2.5,97.5])
ax[1].hist(res["mu"][:,0],bins=40,color=BLUE,alpha=.8); ax[1].axvline(mu[0],color=GREEN,lw=2,label="truth")
ax[1].axvline(em_mu[0],color=ORANGE,lw=2,ls=":",label="EM MLE"); ax[1].axvspan(lo,hi,color=BLUE,alpha=.15)
ax[1].set_title(r"posterior of $\mu_1$"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("The augmentation posterior centres on the truth and agrees with the EM maximum likelihood -- imputation")
print("recovers what deletion could not, because it uses the observed entries of the incomplete rows too.")
mu : true [ 0. 2. -1.] DA [-0.04 1.99 -0.98] EM(MLE) [-0.04 1.99 -0.98] diag Sigma: true [1. 1.5 2. ] DA [1. 1.47 2.03]
The augmentation posterior centres on the truth and agrees with the EM maximum likelihood -- imputation recovers what deletion could not, because it uses the observed entries of the incomplete rows too.
4. Where complete-case deletion goes wrong — MAR¶
The important failure. We make one column MAR: it goes missing when a correlated, fully observed column is high. Because the two are correlated, dropping those rows also drops the high values of the missing column, so the complete-case mean is biased. Imputation — Bayesian DA or frequentist EM — uses the observed correlate to fill the gaps correctly and recovers the truth.
Ymar=M.induce_missing(Y,rng,"MAR",0.45,target=0,driver=1)
cc_mu,_,ncc=M.complete_case(Ymar); da=M.da_gibbs(Ymar,rng,draws=2000,burn=1000); em2,_=M.em_mvn(Ymar)
lo,hi=np.percentile(da["mu"][:,0],[2.5,97.5])
print(f"true mu1 = {mu[0]:.2f}")
print(f"complete-case = {cc_mu[0]:+.2f} (uses only {ncc} of {len(Y)} rows) <-- BIASED")
print(f"DA imputed = {da['mu'].mean(0)[0]:+.2f} 95% CrI [{lo:+.2f}, {hi:+.2f}]")
print(f"EM (MLE) = {em2[0]:+.2f}")
fig,ax=plt.subplots(figsize=(7.5,2.6))
ests={"complete-case\n(listwise)":(cc_mu[0],RED),"EM\n(MLE)":(em2[0],ORANGE),"DA\n(Bayesian)":(da["mu"].mean(0)[0],BLUE)}
for i,(k,(v,c)) in enumerate(ests.items()):
ax.scatter(v,i,color=c,s=90,zorder=3)
if k.startswith("DA"): ax.plot([lo,hi],[i,i],color=c,lw=2)
ax.axvline(mu[0],color=GREEN,lw=2,label="truth (0)"); ax.set_yticks(range(3)); ax.set_yticklabels(list(ests))
ax.set_xlabel(r"estimate of $\mu_1$"); ax.legend(frameon=False,fontsize=8); ax.set_title("Under MAR, deletion is biased; imputation is not")
plt.tight_layout(); plt.show()
print("Complete-case deletion is pulled well off the truth; both imputation methods land on it. This is why 'just")
print("drop the incomplete rows' is dangerous -- MCAR is the only mechanism under which it is safe.")
true mu1 = 0.00 complete-case = -0.30 (uses only 2162 of 4000 rows) <-- BIASED DA imputed = -0.02 95% CrI [-0.06, +0.03] EM (MLE) = -0.02
Complete-case deletion is pulled well off the truth; both imputation methods land on it. This is why 'just drop the incomplete rows' is dangerous -- MCAR is the only mechanism under which it is safe.
5. Real data — imputing airquality¶
The airquality data (New York, 1973) is the textbook missing-data set: of 153 daily records, Ozone is missing on 37 days and Solar.R on 7. It is the same data behind the spline GAM in Bayesian Penalised Splines & Additive Models — there we smoothed it, here we complete it. We model the four continuous variables jointly (log-scaling the two positive, skewed concentrations toward normality) and impute the gaps, carrying the imputation uncertainty forward.
aq=pd.read_csv("airquality.csv"); cols=list(aq.columns)
fig,ax=plt.subplots(figsize=(9,2.6)); ax.imshow(aq.isna().T.values,aspect="auto",cmap="Greys",interpolation="none")
ax.set_yticks(range(4)); ax.set_yticklabels(cols); ax.set_xlabel("day"); ax.set_title("airquality missingness map (black = missing)")
for j,c in enumerate(cols): ax.text(155,j,f"{aq[c].isna().sum()} NA",va="center",fontsize=8)
plt.tight_layout(); plt.show()
Z=aq.copy(); Z["Ozone"]=np.log(Z["Ozone"]); Z["Solar.R"]=np.log(Z["Solar.R"]) # log the positive skewed concentrations
A=Z.to_numpy(); aqres=M.da_gibbs(A,rng,draws=3000,burn=1500,m_keep=60)
print(f"joint multivariate-normal imputation on [log Ozone, log Solar.R, Wind, Temp]; {np.isnan(A).sum()} cells filled.")
print("Correlations let each gap borrow strength from the observed variables on the same day (e.g. hot, low-wind")
print("days impute high ozone).")
joint multivariate-normal imputation on [log Ozone, log Solar.R, Wind, Temp]; 44 cells filled. Correlations let each gap borrow strength from the observed variables on the same day (e.g. hot, low-wind days impute high ozone).
# imputed Ozone values (back-transformed) with their imputation uncertainty vs the observed distribution
oz_mask=aq["Ozone"].isna().values; comp=np.array(aqres["completed"]) # (m,N,4)
oz_imp=np.exp(comp[:,oz_mask,0]) # imputed ozone draws
imp_mean=oz_imp.mean(0); lo=np.percentile(oz_imp,2.5,axis=0); hi=np.percentile(oz_imp,97.5,axis=0)
fig,ax=plt.subplots(1,2,figsize=(12,4))
ax[0].hist(aq["Ozone"].dropna(),bins=20,color=BLUE,alpha=.7,label="observed",density=True)
ax[0].hist(imp_mean,bins=20,color=RED,alpha=.6,label="imputed (posterior mean)",density=True)
ax[0].set_xlabel("Ozone (ppb)"); ax[0].set_ylabel("density"); ax[0].set_title("Imputed vs observed Ozone"); ax[0].legend(frameon=False,fontsize=8)
tmp=aq["Temp"].values
ax[1].scatter(tmp[~oz_mask],aq["Ozone"].values[~oz_mask],color=BLUE,s=18,label="observed")
order=np.argsort(tmp[oz_mask]); tm=tmp[oz_mask][order]
ax[1].errorbar(tm,imp_mean[order],yerr=[imp_mean[order]-lo[order],hi[order]-imp_mean[order]],fmt="o",color=RED,ms=4,elinewidth=1,capsize=2,label="imputed +/- 95%")
ax[1].set_xlabel("Temp (F)"); ax[1].set_ylabel("Ozone (ppb)"); ax[1].set_title("Imputed ozone follows the Temp relationship"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("The imputed ozone values track the temperature relationship the model learned, and each carries an honest")
print("uncertainty band -- wide, because ozone is the most-missing, hardest-to-pin variable. Imputation is not a")
print("single number substituted in; it is a distribution.")
The imputed ozone values track the temperature relationship the model learned, and each carries an honest uncertainty band -- wide, because ozone is the most-missing, hardest-to-pin variable. Imputation is not a single number substituted in; it is a distribution.
# a concrete estimate: mean Ozone -- complete-case throws away 42 days, DA keeps them (with honest uncertainty)
oz_full=aq["Ozone"].dropna().mean()
oz_da=np.exp(comp[:,:,0]).mean(1) # mean ozone per completed dataset
print(f"mean Ozone -- observed-only average: {oz_full:.1f} ppb (ignores the 37 missing days)")
print(f" -- DA imputed: {oz_da.mean():.1f} ppb 95% CrI [{np.percentile(oz_da,2.5):.1f}, {np.percentile(oz_da,97.5):.1f}]")
cell = np.exp(comp[:, np.isnan(A[:, 0]), 0]).mean(1) # the IMPUTED days only
print(" -- imputed days only: %.1f ppb 95%% CrI [%.1f, %.1f]"
% (cell.mean(), np.percentile(cell, 2.5), np.percentile(cell, 97.5)))
lo_c, hi_c = np.percentile(cell, [2.5, 97.5])
print("\nRead that carefully rather than as confirmation. The overall mean barely moves (%.1f against"
% oz_da.mean())
print("%.1f) because only %d of %d days are missing, so the observed days dominate either way. And the"
% (oz_full, int(np.isnan(A[:, 0]).sum()), len(A)))
print("imputed days themselves average %.1f ppb -- below the observed %.1f, in the direction MAR predicts,"
% (cell.mean(), oz_full))
print("but by only %.1f ppb and with a 95%% interval of [%.1f, %.1f] that comfortably contains the observed"
% (oz_full - cell.mean(), lo_c, hi_c))
print("mean. On these data that is suggestive, not decisive: 37 imputed days cannot pin down much.")
print("\nAmelia in the R notebook puts the same quantity nearer 39 ppb across its five imputations, a")
print("little lower than this sampler. Two correct implementations differing by a couple of ppb on a")
print("quantity this uncertain is the expected outcome, and it is worth more than either number alone --")
print("what neither engine supports is treating an imputed mean as a precise estimate.")
mean Ozone -- observed-only average: 42.1 ppb (ignores the 37 missing days)
-- DA imputed: 41.7 ppb 95% CrI [39.9, 44.3]
-- imputed days only: 40.5 ppb 95% CrI [33.0, 50.9]
Read that carefully rather than as confirmation. The overall mean barely moves (41.7 against
42.1) because only 37 of 153 days are missing, so the observed days dominate either way. And the
imputed days themselves average 40.5 ppb -- below the observed 42.1, in the direction MAR predicts,
but by only 1.6 ppb and with a 95% interval of [33.0, 50.9] that comfortably contains the observed
mean. On these data that is suggestive, not decisive: 37 imputed days cannot pin down much.
Amelia in the R notebook puts the same quantity nearer 39 ppb across its five imputations, a
little lower than this sampler. Two correct implementations differing by a couple of ppb on a
quantity this uncertain is the expected outcome, and it is worth more than either number alone --
what neither engine supports is treating an imputed mean as a precise estimate.
6. The limit of ignorability — MNAR¶
Data augmentation is correct only when the mechanism is ignorable (MCAR/MAR). If missingness depends on the unobserved value itself — MNAR — imputing from the observed data cannot recover the truth, because the observed data contain no information about how the missing values differ. Here we make a column MNAR (high values hidden) and watch even the Bayesian imputation stay biased.
Ymnar=M.induce_missing(Y,rng,"MNAR",0.45,target=0)
cc3,_,_=M.complete_case(Ymnar); da3=M.da_gibbs(Ymnar,rng,draws=1500,burn=800)
print(f"true mu1 = {mu[0]:.2f}")
print(f"complete-case = {cc3[0]:+.2f} DA imputed = {da3['mu'].mean(0)[0]:+.2f} <-- both still biased")
print("Under MNAR the observed data are silent about the missing values' distribution, so no ignorable method --")
print("Bayesian or frequentist -- can fix it. Recovering here needs a MODEL of the missingness itself: selection")
print("models (Project 4) and pattern-mixture models (Project 5), with an explicit, untestable assumption.")
true mu1 = 0.00 complete-case = -0.51 DA imputed = -0.37 <-- both still biased Under MNAR the observed data are silent about the missing values' distribution, so no ignorable method -- Bayesian or frequentist -- can fix it. Recovering here needs a MODEL of the missingness itself: selection models (Project 4) and pattern-mixture models (Project 5), with an explicit, untestable assumption.
7. Cross-check in PyMC¶
PyMC imputes automatically: pass the data as a masked array and any missing entry becomes a latent variable (y_unobserved) sampled with the rest. We fit the four-variable multivariate normal with an LKJ prior on the correlation and compare its imputed cells with the from-scratch sampler.
import pymc as pm
Am=np.ma.masked_invalid(A)
with pm.Model() as mod:
muv=pm.Normal("mu",0,10,shape=4)
chol,_,_=pm.LKJCholeskyCov("C",n=4,eta=2,sd_dist=pm.HalfNormal.dist(5),compute_corr=True)
pm.MvNormal("y",mu=muv,chol=chol,observed=Am)
idata=pm.sample(700,tune=1000,chains=4,target_accept=0.9,random_seed=3,progressbar=False)
mu_pm=idata.posterior["mu"].mean(("chain","draw")).values
print("posterior mean of mu, in the units each variable was modelled on")
print(" [log Ozone, log Solar.R, Wind (mph), Temp (F)]")
print("from-scratch", aqres["mu"].mean(0).round(2), " PyMC", mu_pm.round(2))
# compare imputed cells
miss_rc=np.argwhere(np.isnan(A)); da_cell=aqres["imputed_mean"][np.isnan(A)]
pm_cell=idata.posterior["y_unobserved"].mean(("chain","draw")).values
fig,ax=plt.subplots(figsize=(5.2,5)); ax.scatter(da_cell,pm_cell,color=BLUE,s=20)
lim=[min(da_cell.min(),pm_cell.min()),max(da_cell.max(),pm_cell.max())]; ax.plot(lim,lim,"k--",lw=1)
ax.set_xlabel("from-scratch imputed cell"); ax.set_ylabel("PyMC imputed cell"); ax.set_title(f"Imputed values agree (r={np.corrcoef(da_cell,pm_cell)[0,1]:.3f})")
plt.tight_layout(); plt.show()
print("The two samplers impute the same cell values and recover the same mean vector: PyMC's automatic masked-array")
print("imputation is the data-augmentation loop under the hood.")
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 y 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: [mu, C, y_unobserved]
Sampling 4 chains for 1_000 tune and 700 draw iterations (4_000 + 2_800 draws total) took 6 seconds.
posterior mean of mu, in the units each variable was modelled on
[log Ozone, log Solar.R, Wind (mph), Temp (F)]
from-scratch [ 3.42 5. 9.95 77.9 ] PyMC [ 3.39 4.98 10.02 77.46]
The two samplers impute the same cell values and recover the same mean vector: PyMC's automatic masked-array imputation is the data-augmentation loop under the hood.
8. Summary¶
Missing data is not a nuisance to be deleted but an inference problem. Rubin's taxonomy says when you can ignore the why: under MCAR/MAR the mechanism is ignorable and imputing from the data model is valid; under MNAR it is not. Complete-case deletion — the silent default — is unbiased only under MCAR and always discards information; under MAR it is biased, as the simulation and the airquality mean showed. The Bayesian remedy is data augmentation: model the data, treat each missing value as an unknown, and Gibbs-sample it from its conditional given the observed entries and the current parameters — recovering the truth under MAR with honest, correlation-aware uncertainty, and matching both the EM maximum likelihood and PyMC's automatic masked-array imputation.
The connections run throughout the portfolio. The impute-then-update loop is the Albert–Chib augmentation (latent $z$), the Tobit censored-value draw, and the LCA latent class — missing data is the general frame, those were special cases. The NIW parameter update is the same conjugate multivariate-normal posterior used in the asset-risk and multivariate-regression projects. And the MNAR wall we hit here is exactly what the next projects climb: Project 2 turns these single-model draws into practical multiple imputation by chained equations (Rubin's rules, mice); Projects 4–5 confront non-ignorable missingness with selection and pattern-mixture models — where, as with the Tobit selection you have already met, the fix requires an explicit and untestable assumption about why the data are missing.