Categorical Missing Data & Nonresponse

Partially Classified, Not Missing

The arc has run on continuous outcomes throughout; this one turns to categorical data — contingency tables and surveys where a respondent answers one question and skips another. Such a unit is not missing, it is partially classified: we know which cells of the table it could belong to, just not which one. The Bayesian tool is the categorical twin of the Gaussian data augmentation that opened the arc — a Dirichlet–multinomial Gibbs sampler alternating an allocation step with a conjugate draw.

The I-step allocates each partial unit across the cells it is compatible with, in proportion to the current $p$ — a unit with $Y_1=i$ and $Y_2$ missing is split between $(i,1)$ and $(i,2)$, one missing both is spread over all four. The P-step then completes the table and draws $p\sim\text{Dirichlet}(\text{counts}+\alpha)$. Under MAR the two steps are the exact conditionals of the right posterior, so units missing both answers wash out entirely — the allocation and the marginalisation cancel, as they must.

Does it recover a known table?

Before any of that is trusted it is checked against a known answer. Four thousand units are drawn from a 2×2 with a true odds ratio of 9.33, then $Y_2$ is knocked out at random for 1,794 of them — nearly half the sample, which a complete-case analysis would simply delete. The augmentation keeps all 4,000 and returns a posterior median of 10.05 with a 95% interval of [8.30, 12.27], covering the truth, alongside a marginal of 0.56 against a true 0.55. The sampler does what it claims; the interesting question is what the retained units are actually worth, which the real data answer rather more sharply.

What the partial answers are worth

Congdon's survey has 315 respondents answering two yes/no questions, of whom only 224 answered both: 75 answered $Y_1$ alone, 2 answered $Y_2$ alone, and 14 skipped both. The natural claim is that folding in those 91 partial answers buys a better estimate. Put the complete-case analysis on the same Bayesian footing and it turns out to be half right, in a way worth being precise about. The two questions are strongly associated — an odds ratio near 8 — and that is settled entirely by the 224. Not one partial unit answered both questions, so none of them carries any information about the association, and the odds-ratio posterior is unmoved: median 7.72 against 7.77, intervals of width 11.9 and 11.6.

Where the partial answers do pay is the marginal of the question they answered. $P(Y_1=\text{yes})$ moves from 0.544 to 0.571 and narrows by 15% on 75 extra $Y_1$ answers — a genuine shift, since 65% of the $Y_2$-skippers said yes to $Y_1$ against 54% of the completers. $P(Y_2=\text{yes})$, with only 2 extra answers, does not move at all. Partial records pay off exactly where they carry data, and nowhere else. R's cat reaches the same verdict by a different route: em.cat on all 315 returns an odds ratio of 7.79 against 7.81 on the 224, and 20 proper multiple imputations pooled by Rubin's rules give 7.64.

When the skipping is not at random

All of that assumes the skips were at random. People may skip a question because of their answer, which is untestable, so it becomes a knob: a sensitivity parameter $\psi$ tilting the non-responders' category distribution away from the responders', with $\psi=1$ recovering MAR. Sweeping it from 0.4 to 3.0 moves $P(Y_2=\text{yes})$ across 0.30 to 0.44 around a MAR value of 0.35 — the categorical echo of the pattern-mixture $\delta$, and the same conclusion: the honest report is the curve, not the single number.

Where the arc meets latent-class analysis

The payoff is a connection. This is how latent class analysis handles a missing item: each subject's likelihood is a product over items, a missing item drops out of the product, and class membership is inferred by marginalising over the unobserved response — the same sum over compatible cells the I-step performs. The model's allocation for a $Y_1=\text{yes}$ unit with $Y_2$ missing comes out at 0.529 toward yes; among those who answered both, 0.533 did say yes. MAR is precisely the assumption that those two numbers agree, which makes the assumption visible rather than buried.

Three engines on the same table

Three independent engines agree. The from-scratch sampler is checked against PyMC, where the same ignorable model is written not as an augmentation but as a Dirichlet feeding one multinomial per missingness pattern — the completers over $p$, the $Y_2$-missing units over the row margins, the $Y_1$-missing over the column margins, and the both-missing group contributing nothing at all, which is the marginalisation made explicit. It returns an odds ratio of 7.67 against the augmentation's 7.77. In R, Schafer's cat — the categorical counterpart of the norm package that opened the arc — runs the whole prelim.cat → em.cat → da.cat pipeline to a posterior mean of 8.36, 95% CrI [4.05, 15.50], and imp.cat supplies proper multiple imputations whose pooled log odds ratio gives 7.64, CI [3.75, 15.55], with a fraction of missing information of 0.32.

Notebooks

Downloads

Categorical Missing-Data Module — Source Code

"""
catmiss.py -- MISSING CATEGORICAL DATA and nonresponse (from scratch).

Backs the notebooks in  "Categorical Missing Data & Nonresponse".

The whole arc has run on continuous outcomes; this final project moves to CATEGORICAL data --
contingency tables and surveys where some units are only PARTIALLY CLASSIFIED (one item answered,
another left blank). 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 (a 2x2 here), the conjugate loop is:

  I-step (augment): allocate each partially-classified unit to the full cells it is compatible with,
                    in proportion to the current p (a multinomial draw). A unit with Y1 = i but Y2
                    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 now complete, draw p ~ Dirichlet(counts + prior).

Under MAR this is exact and ignorable -- the partial units are allocated exactly as the observed
conditional distribution dictates, and using them SHARPENS the estimate over a complete-case
analysis that would throw them away.

NONIGNORABLE (MNAR) nonresponse -- the subject of Congdon's categorical missing-data chapter -- is
when the chance of answering DEPENDS on the true category (people with something to hide skip the
question). We encode it with a sensitivity parameter psi that TILTS the allocation of the missing
units away from the observed conditional: psi = 1 is MAR, psi != 1 makes the non-responders'
category distribution differ from the responders'. As in the pattern-mixture project, psi is not
identified by the data -- it is a knob for a sensitivity analysis.

CONNECTION TO LATENT CLASS ANALYSIS. An LCA (the Latent Class Analysis arc) handles a missing item on
exactly this principle: the item simply drops out of that subject's likelihood, i.e. the class
membership is inferred by marginalising over the unobserved response -- the same "sum over the cells
the unit is compatible with" that the I-step performs. Categorical missing data is where the
missing-data arc and the latent-class arc meet.
"""

import numpy as np


def augment_gibbs(nboth, ny2mis, ny1mis, nbothmis, rng, draws=4000, burn=2000,
                  psi1=1.0, psi2=1.0, alpha=0.5):
    """Dirichlet-multinomial augmentation for a 2x2 table with item nonresponse.
      nboth   : 2x2 complete-case counts (rows Y1, cols Y2)
      ny2mis  : length-2, units with Y1=i observed but Y2 missing
      ny1mis  : length-2, units with Y2=j observed but Y1 missing
      nbothmis: scalar, units missing both items
      psi1,psi2: MNAR sensitivity tilts (1 = MAR) for Y1, Y2 nonresponse
    Returns posterior draws of the 2x2 cell-probability matrix p."""
    nboth = np.asarray(nboth, float); ny2mis = np.asarray(ny2mis, float); ny1mis = np.asarray(ny1mis, float)
    t1 = np.array([1.0, psi1]); t2 = np.array([1.0, psi2])
    p = np.ones((2, 2)) / 4; P = np.empty((draws, 2, 2))
    for it in range(draws + burn):
        counts = nboth.copy()
        for i in (0, 1):                                    # Y1=i known, Y2 missing -> split over Y2
            w = p[i, :] * t2; counts[i, :] += rng.multinomial(int(ny2mis[i]), w / w.sum())
        for j in (0, 1):                                    # Y2=j known, Y1 missing -> split over Y1
            w = p[:, j] * t1; counts[:, j] += rng.multinomial(int(ny1mis[j]), w / w.sum())
        w = (p * np.outer(t1, t2)).ravel()                  # both missing -> spread over all 4
        counts += rng.multinomial(int(nbothmis), w / w.sum()).reshape(2, 2)
        g = rng.gamma(counts + alpha); p = g / g.sum()      # p ~ Dirichlet(counts + alpha)
        if it >= burn:
            P[it - burn] = p
    return P


def odds_ratio(P):
    """odds ratio of the 2x2 association from probability draws P:(draws,2,2)."""
    return (P[:, 0, 0] * P[:, 1, 1]) / (P[:, 0, 1] * P[:, 1, 0])


def pmarg(P, var, level):
    """marginal probability P(var = level); var in {0(Y1),1(Y2)}, level in {0,1}."""
    return P[:, level, :].sum(1) if var == 0 else P[:, :, level].sum(1)


def complete_case(nboth):
    """complete-case cell proportions (partial units discarded)."""
    nboth = np.asarray(nboth, float); return nboth / nboth.sum()


def simulate_table(N, p, resp2, rng):
    """simulate N units from a 2x2 p, then make Y2 missing with probability that depends on the true
    Y2 value (resp2 = [P(respond|Y2=1), P(respond|Y2=2)]; equal => MAR, unequal => MNAR).
    Returns the four partially-classified count groups."""
    p = np.asarray(p, float); flat = p.ravel()
    draw = rng.multinomial(N, flat).reshape(2, 2)
    nboth = np.zeros((2, 2)); ny2mis = np.zeros(2)
    for i in (0, 1):
        for j in (0, 1):
            nij = draw[i, j]; resp = rng.random(nij) < resp2[j]
            nboth[i, j] = resp.sum(); ny2mis[i] += (~resp).sum()
    return nboth, ny2mis, np.zeros(2), 0

References

• Schafer, J. L. (1997). Analysis of Incomplete Multivariate Data. Chapman & Hall. — chapters 7–8 on categorical data; the method behind the cat package
• Congdon, P. (2010). Applied Bayesian Hierarchical Methods. Chapman & Hall/CRC. — the 2×2 nonresponse table and the nonignorable treatment
• Little, R. J. A. & Rubin, D. B. (2019). Statistical Analysis with Missing Data, 3rd ed. Wiley. — chapter 13 on partially classified contingency tables
• Chen, T. & Fienberg, S. E. (1974). Two-dimensional contingency tables with both completely and partially cross-classified data. Biometrics 30(4), 629–642. — the original treatment of exactly this problem
• Rubin, D. B. (1987). Multiple Imputation for Nonresponse in Surveys. Wiley. — the pooling rules used by mi.inference