Categorical Missing Data & Nonresponse (R)¶
cat — Schafer's data augmentation for incomplete categorical data¶
The R counterpart to catmiss_python.ipynb. Where the Python notebook builds the Dirichlet–multinomial augmentation from scratch, this notebook uses cat (Schafer) — the categorical analog of the norm package from Project 1's foundations. It implements exactly the same method: prelim.cat sorts the missingness patterns, em.cat finds the maximum-likelihood cell probabilities, and da.cat runs the data-augmentation Gibbs. Same data: Congdon's 2×2 survey table with item nonresponse.
options(repr.plot.width=8, repr.plot.height=4)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages(library(cat))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; PURP<-"#6b46c1"
# build the incomplete 2x2 data (Y1,Y2 in {1,2}, NA = not answered)
x<-rbind(matrix(rep(c(1,1),89),ncol=2,byrow=TRUE), matrix(rep(c(1,2),13),ncol=2,byrow=TRUE),
matrix(rep(c(2,1),57),ncol=2,byrow=TRUE), matrix(rep(c(2,2),65),ncol=2,byrow=TRUE),
matrix(rep(c(1,NA),26),ncol=2,byrow=TRUE), matrix(rep(c(2,NA),49),ncol=2,byrow=TRUE),
matrix(rep(c(NA,1),2),ncol=2,byrow=TRUE), matrix(rep(c(NA,NA),14),ncol=2,byrow=TRUE))
cat("total respondents:", nrow(x), " answered both:", sum(complete.cases(x)),
" partial:", sum(!complete.cases(x)), "\n")
total respondents: 315 answered both: 224 partial: 91
1. cat: EM then data augmentation¶
prelim.cat prepares the incomplete table, em.cat returns the maximum-likelihood cell probabilities (using the partially-classified units), and da.cat draws from the posterior by data augmentation. We read off the association (odds ratio) between the two questions and compare with the complete-case value.
s<-prelim.cat(x); th<-em.cat(s, showits=FALSE)
OR_em<-(th[1,1]*th[2,2])/(th[1,2]*th[2,1])
cc<-table(x[complete.cases(x),1], x[complete.cases(x),2]); OR_cc<-(cc[1,1]*cc[2,2])/(cc[1,2]*cc[2,1])
cat("EM cell probabilities:\n"); print(round(th,3))
cat(sprintf("\nodds ratio: complete-case %.2f (n=%d) EM (all %d) %.2f\n", OR_cc, sum(cc), nrow(x), OR_em))
# Bayesian DA posterior of the odds ratio
rngseed(1234); ths<-th; ORs<-numeric(2000)
for(k in 1:2000){ ths<-da.cat(s, ths, steps=4); ORs[k]<-(ths[1,1]*ths[2,2])/(ths[1,2]*ths[2,1]) }
cat(sprintf("DA (Bayesian) odds ratio: %.2f 95%% CrI [%.2f, %.2f]\n", mean(ORs), quantile(ORs,.025), quantile(ORs,.975)))
hist(ORs, breaks=40, col=adjustcolor(BLUE,.7), xlab="odds ratio between the two questions", main="cat: data-augmentation posterior of the association")
abline(v=OR_cc, col=RED, lwd=2, lty=2); abline(v=mean(ORs), col=BLUE, lwd=2)
legend("topright", c("complete-case","DA posterior mean"), col=c(RED,BLUE), lwd=2, lty=c(2,1), bty="n")
cat("cat's EM and DA reproduce the from-scratch augmentation: the two questions are strongly associated (OR ~8).\n")
cat(sprintf("Note how little the 91 partial records move it -- EM on all 315 gives %.2f against %.2f on the 224 who\n", OR_em, OR_cc))
cat("answered both. None of the partial units answered BOTH questions, so none of them carries information about\n")
cat("the association; what they sharpen is the marginal of the question they did answer.\n")
EM cell probabilities:
[,1] [,2] [1,] 0.375 0.054 [2,] 0.269 0.302
odds ratio: complete-case 7.81 (n=224) EM (all 315) 7.79
DA (Bayesian) odds ratio: 8.36 95% CrI [4.05, 15.50]
cat's EM and DA reproduce the from-scratch augmentation: the two questions are strongly associated (OR ~8).
Note how little the 91 partial records move it -- EM on all 315 gives 7.79 against 7.81 on the 224 who
answered both. None of the partial units answered BOTH questions, so none of them carries information about
the association; what they sharpen is the marginal of the question they did answer.
2. Multiple imputation with imp.cat¶
cat also produces proper multiple imputations: imp.cat fills the missing categories with a draw from their predictive distribution, and mi.inference combines analyses across imputations by Rubin's rules. We impute the missing answers repeatedly, confirm the imputed marginal matches the model, and pool the log odds ratio across imputations by Rubin's rules.
rngseed(99); ths<-em.cat(s, showits=FALSE); py2<-numeric(20); est<-numeric(20); se<-numeric(20)
for(k in 1:20){
ths<-da.cat(s, ths, steps=20); xi<-imp.cat(s, ths); py2[k]<-mean(xi[,2]==2)
n<-table(factor(xi[,1],1:2), factor(xi[,2],1:2)) # completed 2x2 for this imputation
est[k]<-log((n[1,1]*n[2,2])/(n[1,2]*n[2,1])); se[k]<-sqrt(sum(1/n))
}
cat(sprintf("imputed P(Y2 = yes) across imputations: mean %.3f (EM model %.3f)\n", mean(py2), sum(th[,2])))
mi<-mi.inference(as.list(est), as.list(se))
cat(sprintf("Rubin-pooled log odds ratio %.3f (se %.3f) -> OR %.2f, 95%% CI [%.2f, %.2f], FMI %.2f\n",
mi$est, mi$std.err, exp(mi$est), exp(mi$lower), exp(mi$upper), mi$fminf))
cat(sprintf("complete-case OR %.2f; the imputations move the association by %.0f%%, because none of the 91 partial\n",
OR_cc, 100*abs(exp(mi$est)-OR_cc)/OR_cc))
cat("units answered both questions -- multiple imputation recovers the margins, not the association.\n")
cat("cat is the categorical counterpart of norm (Project 1) and Amelia -- Schafer's data augmentation for\n")
cat("incomplete categorical data, with imp.cat/mi.inference supplying proper Rubin-rules multiple imputation.\n")
imputed P(Y2 = yes) across imputations: mean 0.355 (EM model 0.356)
Rubin-pooled log odds ratio 2.033 (se 0.360) -> OR 7.64, 95% CI [3.75, 15.55], FMI 0.32
complete-case OR 7.81; the imputations move the association by 2%, because none of the 91 partial
units answered both questions -- multiple imputation recovers the margins, not the association.
cat is the categorical counterpart of norm (Project 1) and Amelia -- Schafer's data augmentation for
incomplete categorical data, with imp.cat/mi.inference supplying proper Rubin-rules multiple imputation.
3. Summary¶
cat is the field-standard tool for incomplete categorical data and the direct categorical analog of the norm package used in Project 1. Its prelim.cat → em.cat → da.cat pipeline is the same maximum-likelihood-then-data-augmentation workflow, here on a contingency table: it estimates the cell probabilities in Congdon's 2×2 survey from all 315 respondents rather than the 224 who answered both — sharpening the marginals, though the association itself is carried entirely by the completers — reproducing the from-scratch Dirichlet–multinomial sampler in catmiss_python.ipynb, and imp.cat/mi.inference deliver Rubin-rules multiple imputation.
This closes the missing-data arc on the frequentist/package side: norm and Amelia for continuous data (Project 1), mice for chained equations (Projects 2–3, 5), sampleSelection for selection models (Project 4), and cat for categorical data here. All rest on the ignorable (MAR) assumption; the nonignorable case is a sensitivity analysis, as the Python notebook's $\psi$ sweep shows. And the augmentation that underlies cat is the same operation a latent class model uses for a missing item — the point where this arc joins the latent-class arc.