Missing Data — Foundations (R)¶
norm (Schafer's data augmentation) and Amelia (EMB multiple imputation)¶
The R counterpart to foundations_python.ipynb. Where the Python notebook builds the multivariate-normal data-augmentation sampler from scratch, this notebook uses the packages that made it standard: norm — Schafer's implementation of EM and data augmentation for the incomplete multivariate normal — and Amelia — King et al.'s EMB (expectation-maximisation with bootstrapping) engine for fast multiple imputation. Same data: the classic airquality set, imputed under the ignorable (MAR) assumption.
options(repr.plot.width=11, repr.plot.height=4.5, warn=-1)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages({library(norm); library(Amelia)})
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; GREY<-"#718096"
aq<-read.csv("airquality.csv")
cat("airquality:", nrow(aq), "rows; NA per column:\n"); print(colSums(is.na(aq)))
airquality: 153 rows; NA per column:
Ozone Solar.R Wind Temp
37 7 0 0
1. Schafer's norm — EM then data augmentation¶
norm is the reference implementation of the method the Python notebook builds by hand. prelim.norm sorts the missing pattern, em.norm finds the maximum-likelihood $(\mu,\Sigma)$, and da.norm runs the data-augmentation Gibbs from there. We log-scale the two positive, skewed concentrations first, exactly as in Python.
Z<-aq; Z$Ozone<-log(Z$Ozone); Z$Solar.R<-log(Z$Solar.R)
x<-as.matrix(Z); s<-prelim.norm(x)
thetahat<-em.norm(s, showits=FALSE) # ML estimate of (mu, Sigma)
mle<-getparam.norm(s, thetahat)
cat("EM maximum-likelihood mean (log Ozone, log Solar.R, Wind, Temp):\n"); print(round(mle$mu,2))
rngseed(1234)
theta<-da.norm(s, thetahat, steps=2000, showits=FALSE) # data augmentation from the MLE
da<-getparam.norm(s, theta)
cat("\ndata-augmentation posterior mean:\n"); print(round(da$mu,2))
cat("\nEM and DA agree, and match the from-scratch sampler in the Python notebook. imp.norm draws a completed\n")
cat("dataset; back-transforming exp() gives imputed ozone on the original ppb scale.\n")
ximp<-imp.norm(s, theta, x)
oz_obs<-aq$Ozone[!is.na(aq$Ozone)]; oz_imp<-exp(ximp[is.na(aq$Ozone),"Ozone"])
hist(oz_obs, breaks=20, col=adjustcolor(BLUE,.6), xlim=range(c(oz_obs,oz_imp)), freq=FALSE, xlab="Ozone (ppb)", main="airquality: observed vs imputed Ozone (norm)")
hist(oz_imp, breaks=15, col=adjustcolor(RED,.6), freq=FALSE, add=TRUE)
legend("topright", c("observed","imputed"), fill=c(BLUE,RED), bty="n")
EM maximum-likelihood mean (log Ozone, log Solar.R, Wind, Temp):
[1] 3.42 5.00 9.96 77.88
data-augmentation posterior mean:
[1] 3.38 4.92 9.95 78.17
EM and DA agree, and match the from-scratch sampler in the Python notebook. imp.norm draws a completed
dataset; back-transforming exp() gives imputed ozone on the original ppb scale.
2. Amelia — multiple imputation by EMB¶
Amelia produces $m$ complete datasets at once by bootstrapping the EM estimate (the EMB algorithm), the practical multiple-imputation workflow. Its overimputation diagnostic hides observed values, re-imputes them, and checks that the true values fall inside the imputation intervals — a direct test of the imputation model.
set.seed(1)
am<-amelia(aq, m=5, logs=c("Ozone","Solar.R"), p2s=0) # 5 imputations; log the skewed concentrations
cat("Amelia produced", am$m, "completed datasets.\n")
cat("imputed-Ozone means across the 5 imputations:", round(sapply(am$imputations, function(d) mean(d$Ozone[is.na(aq$Ozone)])),1), "\n")
cat("observed-only Ozone mean:", round(mean(aq$Ozone,na.rm=TRUE),1), "-- the imputed means differ because the\n")
cat("missing days were not a random slice (MAR), and vary across imputations, reflecting imputation uncertainty.\n")
overimpute(am, var="Ozone")
Amelia produced 5 completed datasets.
imputed-Ozone means across the 5 imputations: 39.2 38.2 41.2 38.8 38.9
observed-only Ozone mean: 42.1 -- the imputed means differ because the
missing days were not a random slice (MAR), and vary across imputations, reflecting imputation uncertainty.
3. Summary¶
The standard R packages reproduce the from-scratch results. norm is Schafer's implementation of exactly the method the Python notebook builds — em.norm for the maximum-likelihood $(\mu,\Sigma)$, da.norm for the data-augmentation posterior, imp.norm for a completed dataset — and it recovers the same mean vector and imputes the same ozone distribution. Amelia wraps the practical multiple-imputation workflow: its EMB algorithm bootstraps the EM fit to yield several completed datasets quickly, with an overimputation diagnostic that validates the model against held-out observed values.
norm and Amelia are the field-standard tools for multivariate-normal missing data and the frequentist/likelihood mirror of the Bayesian data augmentation in foundations_python.ipynb. The multiple imputations Amelia produces here are the raw material for Rubin's rules, which the next project (Multiple Imputation by Chained Equations) formalises with mice. All of this rests on the ignorable (MAR) assumption; non-ignorable (MNAR) missingness is the subject of the later selection and pattern-mixture projects.