Missing Data — Foundations

Python · PyMC · R  ·  Download missing-data module

Treat the Gap as a Parameter

The Bayesian view of missing data is almost anticlimactic: a missing value is just another unknown, so give it a distribution and sample it alongside the parameters. That is data augmentation — the same move already used for the latent zz in Albert–Chib probit, the censored draws in Tobit, and the class label in latent class analysis. Those were special cases; missing data is the general one. Under a multivariate normal model the loop is two conjugate steps: impute each row's missing entries from their conditional normal given what was observed, then draw (μ,Σ)(\mu,\Sigma) from the completed data.

ymis∣yobs∼N ⁣(μmis+B(yobs−μobs),  Σmis−BΣobs,mis),B=Σmis,obsΣobs,obs−1y_{\text{mis}}\mid y_{\text{obs}}\sim N\!\big(\mu_{\text{mis}}+B(y_{\text{obs}}-\mu_{\text{obs}}),\;\Sigma_{\text{mis}}-B\Sigma_{\text{obs},\text{mis}}\big),\qquad B=\Sigma_{\text{mis},\text{obs}}\Sigma_{\text{obs},\text{obs}}^{-1}

Why it is allowed

What makes the loop legitimate is Rubin's taxonomy. Missingness is MCAR if it is independent of everything, MAR if it depends only on observed data, and MNAR if it depends on the values that are missing. Under MCAR and MAR the mechanism is ignorable — you may impute from the data model alone and never model why data went missing. That single word is what the whole first half of this arc rests on.

What deletion costs

The cost of ignoring it is shown rather than asserted. Under MAR with 45% missingness on a correlated pair, complete-case deletion estimates a mean of −0.30 against a truth of 0.00, discarding 1,838 of 4,000 rows to get there. Both imputation routes — the Bayesian augmentation and its frequentist twin, EM — land on −0.02. Deletion is unbiased only under MCAR, and it throws away every observed entry that happens to share a row with a gap.

estimate of a mean whose truth is 0.00MCARMARMNAR
complete-case deletionsafe−0.30−0.51
data augmentationsafe−0.02−0.37
EM (maximum likelihood)safe−0.02biased

Where ignorability runs out

The chapter is equally clear about where the method stops. Under MNAR the same sampler estimates −0.37 where complete-case gives −0.51: both still badly biased. When missingness depends on the unobserved value, the observed data are simply silent about the distribution of what is missing, and no ignorable method — Bayesian or frequentist — can recover it. That failure is the reason the arc continues into selection and pattern-mixture models, which require an explicit and untestable assumption about the mechanism.

Real data

On the real `airquality` data — 153 days, 37 ozone and 7 solar values missing — the joint model fills 44 cells, each gap borrowing strength from the variables measured on the same day. Schafer's norm reproduces the posterior mean in R, EM matches the maximum likelihood, and Amelia supplies five completed datasets by a different algorithm entirely.

Reading the Evidence at Its Actual Strength

Two corrections. The PyMC comparison labelled its mean vector "log-scale" when only two of its four entries are logged — 9.95 and 77.9 are plainly wind speed in mph and temperature in Fahrenheit. And the notebook claimed the imputed mean "shifts" on the strength of 42.1 → 41.7 ppb, which is a change of 0.4 that the 116 observed days swamp. The informative quantity is the imputed days alone, so that is what is reported: 40.5 ppb, 95% CrI [33.0, 50.9], against 42.1 observed. That is a shift in the direction MAR predicts, but by 1.6 ppb with an interval that comfortably contains the observed mean — suggestive, not decisive. Amelia puts the same quantity nearer 39. Thirty-seven imputed days cannot pin down much, and the sign alone should not be read as confirmation.

Notebooks

Python The three mechanisms simulated and contrasted, data augmentation step by step against EM, complete-case bias under MAR, the airquality imputation with the imputed cells reported separately from the overall mean, the MNAR failure, and a PyMC masked-array cross-check (Python) Open → R Independent R implementation — Schafer's norm for EM followed by data augmentation, and Amelia's bootstrap-EM multiple imputation as a second algorithm (R) Open →

Downloads

md_missdata.py Mechanism simulation for MCAR, MAR and MNAR; Schafer data-augmentation Gibbs for an incomplete multivariate normal, grouped by missingness pattern; EM for the same likelihood; and complete-case moments (NumPy / SciPy) airquality.csv New York air quality — 153 days with 37 missing ozone and 7 missing solar readings, the standard imputation benchmark

Missing-Data Module — Source Code

"""
missdata.py -- FOUNDATIONS OF BAYESIAN MISSING-DATA IMPUTATION (from scratch).

Backs the notebooks in  "Missing Data -- Foundations".

The Bayesian view of missing data is disarmingly simple: a missing value is just another
unknown, so give it a distribution and sample it alongside the parameters. Under a multivariate
normal model that is Schafer's DATA AUGMENTATION (Little & Rubin; Tanner & Wong) -- the exact
"treat the unknown as a parameter and draw it" move already used for the latent z in Albert-Chib
probit, the censored draws in Tobit, and the latent class in LCA. Missing data is the general case;
those were special cases.

Rubin's taxonomy sets the ground rules. Missingness is MCAR if it is independent of all data, MAR
if it depends only on the OBSERVED data, and MNAR if it depends on the values that are missing.
Under MCAR/MAR the missingness mechanism is IGNORABLE: we may impute from the data model alone and
need not model why data are missing. Complete-case (listwise) deletion is unbiased only under MCAR
and generally wastes information; Bayesian imputation uses every observed entry.

The sampler is a two-step Gibbs on an incomplete N x p multivariate normal Y ~ N(mu, Sigma):

  I-step (impute):   for each row, draw the missing entries from their conditional normal given the
                     observed entries,  y_mis | y_obs ~ N( mu_mis + B (y_obs - mu_obs),  S_mis.obs ),
                     with  B = Sigma_mis,obs Sigma_obs,obs^{-1}  and Schur-complement covariance.
  P-step (params):   given the completed data, draw (mu, Sigma) from the conjugate normal-inverse-
                     Wishart posterior (here the standard noninformative limit):
                     Sigma ~ InvWishart(N-1, S),  mu | Sigma ~ N(ybar, Sigma / N).

Rows are grouped by missingness PATTERN so each conditional regression B is formed once per pattern
(Schafer's efficiency trick). Storing spaced completed datasets yields the multiple imputations that
feed Rubin's rules in the next project.
"""

import numpy as np
from scipy.stats import invwishart


# --------------------------------------------------------------------------- #
#  induce missingness under a chosen mechanism                                  #
# --------------------------------------------------------------------------- #

def induce_missing(Y, rng, mechanism="MCAR", frac=0.3, target=0, driver=1):
    """Return a copy of Y with NaNs in column `target` under MCAR / MAR / MNAR.
      MCAR: missing completely at random (rate `frac`).
      MAR : probability of missing rises with the OBSERVED `driver` column.
      MNAR: probability of missing rises with the target column's OWN (would-be) value."""
    Y = np.array(Y, float); N = Y.shape[0]
    if mechanism == "MCAR":
        p = np.full(N, frac)
    elif mechanism == "MAR":
        z = (Y[:, driver] - Y[:, driver].mean()) / Y[:, driver].std()
        p = 1 / (1 + np.exp(-(z * 2)))                       # higher driver -> more missing
    else:  # MNAR
        z = (Y[:, target] - Y[:, target].mean()) / Y[:, target].std()
        p = 1 / (1 + np.exp(-(z * 2)))                       # higher own value -> more missing
    p = p * (frac / p.mean())                                # scale to the requested average rate
    Y[rng.random(N) < np.clip(p, 0, 1), target] = np.nan
    return Y


# --------------------------------------------------------------------------- #
#  simulation                                                                   #
# --------------------------------------------------------------------------- #

def simulate_mvn(N, mu, Sigma, rng):
    return rng.multivariate_normal(np.asarray(mu, float), np.asarray(Sigma, float), size=N)


# --------------------------------------------------------------------------- #
#  data-augmentation Gibbs sampler                                              #
# --------------------------------------------------------------------------- #

def _patterns(Mask):
    """group row indices by identical missingness pattern (tuple of missing columns)."""
    groups = {}
    for i, row in enumerate(Mask):
        groups.setdefault(tuple(np.where(row)[0]), []).append(i)
    return {k: np.array(v) for k, v in groups.items()}


def da_gibbs(Y, rng, draws=2000, burn=1000, m_keep=20):
    """Schafer data-augmentation Gibbs for an incomplete multivariate normal.
    Y : (N,p) with np.nan at missing entries. Returns posterior draws of mu (draws,p) and Sigma
    (draws,p,p), `m_keep` completed datasets (multiple imputations) and the posterior-mean imputation."""
    Y = np.array(Y, float); N, p = Y.shape
    Mask = np.isnan(Y); pat = _patterns(Mask)
    Yc = Y.copy()
    col_mean = np.nanmean(Y, axis=0)
    for j in range(p):
        Yc[Mask[:, j], j] = col_mean[j]                      # start from mean imputation
    MU = np.empty((draws, p)); SIG = np.empty((draws, p, p))
    keep_idx = set(np.linspace(0, draws - 1, m_keep).astype(int))
    completed = []; imp_sum = np.zeros((N, p)); n_imp = 0
    for it in range(draws + burn):
        # ---- P-step: (mu, Sigma) | completed data  (noninformative NIW limit) ----
        ybar = Yc.mean(0); D = Yc - ybar; S = D.T @ D
        Sigma = invwishart.rvs(df=N - 1, scale=S, random_state=rng)
        Sigma = np.atleast_2d(Sigma)
        mu = ybar + np.linalg.cholesky(Sigma / N) @ rng.standard_normal(p)
        # ---- I-step: impute missing entries per pattern ----
        for mis, rows in pat.items():
            if not mis:
                continue
            mis = list(mis); obs = [j for j in range(p) if j not in mis]
            if obs:
                So_o = Sigma[np.ix_(obs, obs)]
                B = Sigma[np.ix_(mis, obs)] @ np.linalg.inv(So_o)
                cmean = mu[mis] + (Yc[np.ix_(rows, obs)] - mu[obs]) @ B.T      # (n_rows, n_mis)
                ccov = Sigma[np.ix_(mis, mis)] - B @ Sigma[np.ix_(obs, mis)]
            else:                                            # whole row missing -> marginal
                cmean = np.tile(mu[mis], (len(rows), 1)); ccov = Sigma[np.ix_(mis, mis)]
            Lc = np.linalg.cholesky(ccov + 1e-10 * np.eye(len(mis)))
            draw = cmean + rng.standard_normal((len(rows), len(mis))) @ Lc.T
            Yc[np.ix_(rows, mis)] = draw
        if it >= burn:
            i = it - burn; MU[i] = mu; SIG[i] = Sigma
            imp_sum += Yc; n_imp += 1
            if i in keep_idx:
                completed.append(Yc.copy())
    return dict(mu=MU, Sigma=SIG, completed=completed, imputed_mean=imp_sum / n_imp, mask=Mask)


# --------------------------------------------------------------------------- #
#  EM (frequentist counterpart) and complete-case moments                       #
# --------------------------------------------------------------------------- #

def em_mvn(Y, tol=1e-6, max_iter=500):
    """EM for the incomplete-multivariate-normal MLE of (mu, Sigma) -- the frequentist twin of the
    Bayesian data augmentation (expectation replaces the imputation draw, maximisation the parameter draw)."""
    Y = np.array(Y, float); N, p = Y.shape; Mask = np.isnan(Y); pat = _patterns(Mask)
    cc = ~Mask.any(axis=1)                                   # init from complete-case moments:
    if cc.sum() > p:                                         # a DIAGONAL Sigma is a degenerate fixed
        mu = Y[cc].mean(0); Sigma = np.cov(Y[cc], rowvar=False)   # point (B=0), so start full-rank
    else:
        mu = np.nanmean(Y, axis=0); Sigma = np.diag(np.nanvar(Y, axis=0))
    for _ in range(max_iter):
        T1 = np.zeros(p); T2 = np.zeros((p, p))
        for mis, rows in pat.items():
            obs = [j for j in range(p) if j not in mis]; mis = list(mis)
            Yr = Y[rows]
            if mis and obs:
                So_o = Sigma[np.ix_(obs, obs)]
                B = Sigma[np.ix_(mis, obs)] @ np.linalg.inv(So_o)
                yhat = Yr.copy()
                yhat[:, mis] = mu[mis] + (Yr[:, obs] - mu[obs]) @ B.T
                ccov = Sigma[np.ix_(mis, mis)] - B @ Sigma[np.ix_(obs, mis)]
            elif mis:
                yhat = np.tile(mu, (len(rows), 1)); ccov = Sigma[np.ix_(mis, mis)]
            else:
                yhat = Yr; ccov = None
            T1 += yhat.sum(0); T2 += yhat.T @ yhat
            if mis and ccov is not None:
                add = np.zeros((p, p)); add[np.ix_(mis, mis)] = ccov * len(rows); T2 += add
        mu_new = T1 / N; Sigma_new = T2 / N - np.outer(mu_new, mu_new)
        if np.max(np.abs(mu_new - mu)) < tol:
            mu, Sigma = mu_new, Sigma_new; break
        mu, Sigma = mu_new, Sigma_new
    return mu, Sigma


def complete_case(Y):
    """listwise-deletion moments -- unbiased only under MCAR."""
    Y = np.array(Y, float); keep = ~np.isnan(Y).any(axis=1); Yc = Y[keep]
    return Yc.mean(0), np.cov(Yc, rowvar=False), keep.sum()

References