Categorical Missing Data & Nonresponse¶

Partially classified tables, and where the arc meets latent-class analysis¶

The whole arc has run on continuous outcomes; the final project turns to categorical data — contingency tables and surveys where a unit may answer one question and skip another. The Bayesian tool is the categorical twin of Project 1's Gaussian data augmentation: a Dirichlet–multinomial Gibbs sampler. For a two-way table with cell probabilities $p$,

  • I-step (augment): allocate each partially-classified unit to the cells it is compatible with, in proportion to the current $p$ (a multinomial draw). A unit with $Y_1=i$ but $Y_2$ missing is split between $(i,1)$ and $(i,2)$; a unit missing both is spread over all four cells.
  • P-step (update): with the table completed, draw $p\sim\text{Dirichlet}(\text{counts}+\alpha)$.

Under MAR this is exact and ignorable, and using the partial units sharpens the estimate over a complete-case analysis that discards them. Nonignorable (MNAR) nonresponse — the chance of answering depends on the true category — is handled as in the pattern-mixture project: a sensitivity parameter $\psi$ that tilts the missing units' category distribution away from the responders', swept to test the conclusion. And there is a payoff connection: this is how latent class analysis copes with a missing item, closing the loop between the missing-data and latent-class arcs. We validate the sampler, run Congdon's 2×2 nonresponse table, do the MNAR sensitivity, make the LCA link concrete, and cross-check in PyMC.

In [1]:
import numpy as np, matplotlib.pyplot as plt
import catmiss as C
rng = np.random.default_rng(7)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("Categorical missing data: allocate partially-classified units across the cells they are compatible with,")
print("in proportion to the current cell probabilities, then redraw the probabilities. Dirichlet-multinomial DA.")

Categorical missing data: allocate partially-classified units across the cells they are compatible with,
in proportion to the current cell probabilities, then redraw the probabilities. Dirichlet-multinomial DA.

1. Does it work? — recovering a table from partial classifications¶

Simulate 4000 units from a known 2×2 distribution, then let $Y_2$ go missing at random for some. The augmentation should recover the true association (odds ratio) and marginals using all the units — the partially-classified ones included — where a complete-case analysis would discard them.

In [2]:
p_true=np.array([[0.35,0.15],[0.10,0.40]]); OR_true=(p_true[0,0]*p_true[1,1])/(p_true[0,1]*p_true[1,0])
nb,n2,n1,nbm=C.simulate_table(4000,p_true,resp2=[0.55,0.55],rng=rng)   # MAR: response independent of Y2
P=C.augment_gibbs(nb,n2,n1,nbm,rng,draws=2500,burn=1200)
orr=C.odds_ratio(P); lo,hi=np.percentile(orr,[2.5,97.5])
print(f"true odds ratio {OR_true:.2f},  true P(Y2=2) {p_true[:,1].sum():.2f}")
print(f"augmented: OR median {np.median(orr):.2f} [{lo:.2f},{hi:.2f}],  P(Y2=2) {C.pmarg(P,1,1).mean():.2f}")
print(f"units used: {int(nb.sum()+n2.sum())} total ({int(nb.sum())} complete + {int(n2.sum())} with Y2 missing, kept not dropped)")
fig,ax=plt.subplots(1,2,figsize=(11,3.6))
ax[0].hist(orr,bins=40,color=BLUE,alpha=.8,density=True); ax[0].axvline(OR_true,color=GREEN,lw=2,label="truth")
ax[0].set_xlabel("odds ratio"); ax[0].set_title("Recovered association"); ax[0].legend(frameon=False,fontsize=8)
ax[1].hist(C.pmarg(P,1,1),bins=40,color=PURP,alpha=.8,density=True); ax[1].axvline(p_true[:,1].sum(),color=GREEN,lw=2)
ax[1].set_xlabel(r"$P(Y_2=2)$"); ax[1].set_title("Recovered marginal"); plt.tight_layout(); plt.show()
print("The augmentation recovers the odds ratio and the marginal, using the partially-classified units instead of")
print("discarding them -- the categorical version of the multivariate-normal data augmentation from Project 1.")

true odds ratio 9.33,  true P(Y2=2) 0.55
augmented: OR median 10.05 [8.30,12.27],  P(Y2=2) 0.56
units used: 4000 total (2206 complete + 1794 with Y2 missing, kept not dropped)
No description has been provided for this image 
The augmentation recovers the odds ratio and the marginal, using the partially-classified units instead of
discarding them -- the categorical version of the multivariate-normal data augmentation from Project 1.

2. Real data — a 2×2 survey with item nonresponse¶

Congdon's categorical missing-data example (BMCD Ch. 11) is a survey of 315 people answering two yes/no questions. Only 224 answered both; the rest skipped one or both: 75 answered $Y_1$ but not $Y_2$, 2 answered $Y_2$ but not $Y_1$, and 14 skipped both. We estimate the association between the two questions using every partial answer, and compare against the complete-case analysis on the same Bayesian footing — a Dirichlet posterior on the 224 — to see exactly what the partial answers are worth.

In [3]:
nboth=np.array([[89.,13.],[57.,65.]]); ny2mis=np.array([26.,49.]); ny1mis=np.array([2.,0.]); nbothmis=14
P=C.augment_gibbs(nboth,ny2mis,ny1mis,nbothmis,rng,draws=20000,burn=4000)
gcc=rng.gamma(nboth.ravel()+0.5,size=(20000,4)); Pcc=(gcc/gcc.sum(1,keepdims=True)).reshape(-1,2,2)
orr=C.odds_ratio(P); orcc=C.odds_ratio(Pcc)
W=lambda x: float(np.diff(np.percentile(x,[2.5,97.5]))[0])
q=lambda x: (np.median(x), np.percentile(x,2.5), np.percentile(x,97.5))
print(f"total {int(nboth.sum()+ny2mis.sum()+ny1mis.sum()+nbothmis)} people; {int(nboth.sum())} answered both, 91 answered only one")
print("odds ratio   complete-case %.2f [%.2f,%.2f]   augmented %.2f [%.2f,%.2f]" % (*q(orcc), *q(orr)))
for v,name,extra in ((0,"P(Y1=yes)","75 more answers"),(1,"P(Y2=yes)","2 more answers")):
    a=C.pmarg(Pcc,v,1); b=C.pmarg(P,v,1)
    print(f"{name}    complete-case {a.mean():.3f} (width {W(a):.3f})   augmented {b.mean():.3f} (width {W(b):.3f})   [{extra}]")
fig,ax=plt.subplots(1,2,figsize=(11,3.6))
ax[0].hist(orcc,bins=60,range=(0,30),color=RED,alpha=.45,density=True,label="complete-case (224)")
ax[0].hist(orr,bins=60,range=(0,30),color=BLUE,alpha=.45,density=True,label="augmented (315)")
ax[0].set_xlabel("odds ratio between the two questions"); ax[0].set_title("Association: the partial answers add nothing")
ax[0].legend(frameon=False,fontsize=8)
ax[1].hist(C.pmarg(Pcc,0,1),bins=50,color=RED,alpha=.45,density=True,label="complete-case")
ax[1].hist(C.pmarg(P,0,1),bins=50,color=BLUE,alpha=.45,density=True,label="augmented")
ax[1].set_xlabel(r"$P(Y_1=\mathrm{yes})$"); ax[1].set_title("Margin of the answered question: shifted and sharpened")
ax[1].legend(frameon=False,fontsize=8); plt.tight_layout(); plt.show()
a1=C.pmarg(Pcc,0,1); b1=C.pmarg(P,0,1)
print(f"The two questions are strongly associated (OR ~{np.median(orr):.0f}), and that is settled by the 224 who answered both.")
print(f"None of the 91 partial units answered BOTH questions, so none of them speaks to the association: the odds-ratio")
print(f"posterior is unmoved ({W(orcc):.1f} wide against {W(orr):.1f}). What they do inform is the margin of the question they")
print(f"did answer -- P(Y1=yes) shifts {a1.mean():.3f} -> {b1.mean():.3f} and narrows by {1-W(b1)/W(a1):.0%} on 75 extra Y1 answers,")
print(f"while P(Y2=yes), with only 2 extra, is unchanged. Partial records pay off exactly where they carry data.")

total 315 people; 224 answered both, 91 answered only one
odds ratio   complete-case 7.72 [4.03,15.89]   augmented 7.77 [4.04,15.66]
P(Y1=yes)    complete-case 0.544 (width 0.130)   augmented 0.571 (width 0.110)   [75 more answers]
P(Y2=yes)    complete-case 0.350 (width 0.124)   augmented 0.357 (width 0.120)   [2 more answers]
No description has been provided for this image 
The two questions are strongly associated (OR ~8), and that is settled by the 224 who answered both.
None of the 91 partial units answered BOTH questions, so none of them speaks to the association: the odds-ratio
posterior is unmoved (11.9 wide against 11.6). What they do inform is the margin of the question they
did answer -- P(Y1=yes) shifts 0.544 -> 0.571 and narrows by 15% on 75 extra Y1 answers,
while P(Y2=yes), with only 2 extra, is unchanged. Partial records pay off exactly where they carry data.

3. Nonignorable nonresponse — a sensitivity analysis¶

The clean answer above assumes the skips were at random. But people may skip a question because of their answer — MNAR. We cannot test this, so we make it a knob: a sensitivity parameter $\psi$ that tilts the missing units toward one category ($\psi=1$ is MAR). Sweeping $\psi$ for the $Y_2$ question shows how much the conclusion about $P(Y_2=\text{yes})$ depends on the untestable nonresponse assumption — the categorical echo of the pattern-mixture $\delta$.

In [4]:
psis=np.linspace(0.4,3.0,20); est=[]; lo_=[]; hi_=[]
for ps in psis:
    Pm=C.augment_gibbs(nboth,ny2mis,ny1mis,nbothmis,rng,draws=1500,burn=800,psi2=ps)
    m=C.pmarg(Pm,1,1); est.append(m.mean()); lo_.append(np.percentile(m,2.5)); hi_.append(np.percentile(m,97.5))
est=np.array(est); lo_=np.array(lo_); hi_=np.array(hi_); mar=est[np.argmin(np.abs(psis-1))]
fig,ax=plt.subplots(figsize=(8,4)); ax.fill_between(psis,lo_,hi_,color=ORANGE,alpha=.2); ax.plot(psis,est,color=ORANGE,lw=2)
ax.axvline(1,color=GREEN,ls="--",lw=1.5,label="psi=1 (MAR)"); ax.scatter([1],[mar],color=GREEN,zorder=3)
ax.set_xlabel(r"$\psi$  =  nonresponse tilt toward $Y_2=$ yes"); ax.set_ylabel(r"$P(Y_2=\mathrm{yes})$")
ax.set_title("MNAR sensitivity: how the estimate moves with the untestable nonresponse assumption"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"P(Y2=yes) ranges from {est.min():.2f} to {est.max():.2f} as psi runs from {psis[0]:.1f} to {psis[-1]:.1f}; the MAR value is {mar:.2f}.")
print("If people who would answer 'yes' were the ones skipping (psi>1), the true rate is higher than MAR suggests.")
print("The data cannot decide psi -- so the honest report is this whole curve, not the single MAR number.")
No description has been provided for this image 
P(Y2=yes) ranges from 0.29 to 0.44 as psi runs from 0.4 to 3.0; the MAR value is 0.35.
If people who would answer 'yes' were the ones skipping (psi>1), the true rate is higher than MAR suggests.
The data cannot decide psi -- so the honest report is this whole curve, not the single MAR number.

4. Where the arc meets latent-class analysis¶

This is exactly how a latent class model (the Latent Class Analysis arc) handles a missing item. In LCA each subject's likelihood is a product over items; a missing item simply drops out of the product, so the class membership is inferred by marginalising over the unobserved response — the same "sum over the cells the unit is compatible with" that the augmentation I-step performs. A subject who skips an item is not discarded; they are classified on the items they did answer. The panel below puts the model's allocation for a $Y_2$-missing unit beside what the completers with the same $Y_1$ actually answered: MAR is precisely the assumption that the two coincide.

In [5]:
phat=P.mean(0)                                                    # estimated 2x2 cell probs
alloc=phat[1,:]/phat[1,:].sum()          # augmentation: a Y1=yes unit with Y2 missing splits by p(Y2 | Y1=yes)
emp=nboth[1,:]/nboth[1,:].sum()          # what the completers with Y1=yes actually answered
print(f"model allocation for a Y1=yes unit with Y2 missing:  P(Y2=no,yes) = {alloc.round(3)}")
print(f"the same conditional among those who answered both:              = {emp.round(3)}")
print("Under MAR the non-responders are assumed to answer like the responders, which is why the two agree; an LCA")
print("marginalising the missing item over its categories performs exactly this allocation.")
fig,ax=plt.subplots(figsize=(6.5,3)); i=np.arange(2)
ax.bar(i-0.18,alloc,0.36,color=PURP,label="augmentation / LCA marginalisation")
ax.bar(i+0.18,emp,0.36,color=GREY,label="observed among completers")
ax.set_xticks(i); ax.set_xticklabels(["Y2=no","Y2=yes"]); ax.set_ylabel("allocation weight")
ax.set_title("Missing item -> marginalise: the LCA principle"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("Marginalising a missing item and augmenting a partially-classified unit are the same operation. Categorical")
print("missing data is the meeting point of this arc and the latent-class arc: both sum over what was not observed.")

model allocation for a Y1=yes unit with Y2 missing:  P(Y2=no,yes) = [0.471 0.529]
the same conditional among those who answered both:              = [0.467 0.533]
Under MAR the non-responders are assumed to answer like the responders, which is why the two agree; an LCA
marginalising the missing item over its categories performs exactly this allocation.
No description has been provided for this image 
Marginalising a missing item and augmenting a partially-classified unit are the same operation. Categorical
missing data is the meeting point of this arc and the latent-class arc: both sum over what was not observed.

5. Cross-check in PyMC¶

The ignorable (MAR) model is a Dirichlet over the four cell probabilities feeding several multinomials — one per missingness pattern, each over the cells that pattern is compatible with (the complete cells use $p$, the $Y_2$-missing units the row margins $P(Y_1)$, the $Y_1$-missing units the column margins). We fit it in PyMC and compare the odds ratio with the from-scratch augmentation.

In [6]:
import pymc as pm, pytensor.tensor as pt
with pm.Model() as mod:
    p=pm.Dirichlet("p",a=np.full(4,0.5)); pm2=p.reshape((2,2))
    pm.Multinomial("both", n=int(nboth.sum()), p=p, observed=nboth.ravel().astype(int))
    pm.Multinomial("y2mis", n=int(ny2mis.sum()), p=pm2.sum(axis=1), observed=ny2mis.astype(int))   # row margins P(Y1)
    pm.Multinomial("y1mis", n=int(ny1mis.sum()), p=pm2.sum(axis=0), observed=ny1mis.astype(int))   # col margins P(Y2)
    idata=pm.sample(1500,tune=1500,chains=4,target_accept=0.9,random_seed=8,progressbar=False)
ps=idata.posterior["p"].values.reshape(-1,4); orr_pm=(ps[:,0]*ps[:,3])/(ps[:,1]*ps[:,2])
print(f"odds ratio (median):  from-scratch {np.median(orr):.2f}   PyMC {np.median(orr_pm):.2f}")
print(f"P(Y2=yes):   from-scratch {C.pmarg(P,1,1).mean():.3f}   PyMC {(ps[:,1]+ps[:,3]).mean():.3f}")
print("PyMC's Dirichlet-multinomial with pattern-specific margins matches the from-scratch augmentation: two ways")
print("to write the same ignorable categorical missing-data model.")

g++ not available, if using conda: `conda install gxx`

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [p]

Sampling 4 chains for 1_500 tune and 1_500 draw iterations (6_000 + 6_000 draws total) took 6 seconds.

odds ratio (median):  from-scratch 7.77   PyMC 7.67
P(Y2=yes):   from-scratch 0.357   PyMC 0.357
PyMC's Dirichlet-multinomial with pattern-specific margins matches the from-scratch augmentation: two ways
to write the same ignorable categorical missing-data model.

6. Summary — and the arc complete¶

Categorical missing data closes the loop. The Dirichlet–multinomial augmentation allocates each partially-classified unit across the cells it is compatible with and redraws the cell probabilities — the categorical twin of the multivariate-normal data augmentation that opened the arc, recovering the association from Congdon's 2×2 survey while using every partial answer the complete-case analysis would throw away. When nonresponse may be nonignorable, a sensitivity parameter $\psi$ traces how the conclusion depends on the untestable assumption, just as the pattern-mixture $\delta$ did for continuous outcomes. And the augmentation's I-step is precisely how latent class analysis treats a missing item — marginalise it out — so this project is where the missing-data and latent-class arcs join.

The six-project journey. Missing data is one idea — treat the unknown as a parameter and sample it — worked through every regime. Ignorable (MCAR/MAR): the multivariate-normal data augmentation (Project 1), the practical MICE multiple imputation with Rubin's rules (Project 2), and missing covariates imputed conditional on the outcome (Project 3). Nonignorable (MNAR): selection models that model why data are missing (Project 4) and pattern-mixture models that make the untestable assumption an explicit sensitivity knob (Project 5). And finally categorical data (Project 6), where the same augmentation meets the latent-class arc. Throughout, the from-scratch samplers are the very augmentation used across this portfolio — Albert–Chib, Tobit, LCA — confirming that missing data was never a special topic but the general frame in which all of them sit.