Pattern-Mixture Models — MNAR by Sensitivity Analysis

The Other Factorisation

Selection models factor the joint density as $p(y)\,p(r\mid y)$. Pattern-mixture models take the other factorisation, $p(y\mid r)\,p(r)$ — stratify by the missingness pattern, give each stratum its own outcome distribution, and mix. The two describe the same joint distribution and cannot be told apart by the data, so the choice is about which assumption you would rather state out loud.

Pattern-mixture's virtue is that it makes the assumption impossible to hide. The outcome distribution among the missing is never observed, so it cannot be estimated — it must be assumed. That assumption becomes a single interpretable knob: a sensitivity parameter $\delta$ saying how much worse the dropouts would have been than otherwise-similar completers. $\delta=0$ is MAR; $\delta\neq0$ is MNAR. Since $\delta$ is not identified, the output is not an estimate but a sensitivity analysis: sweep it and watch the conclusion.

Where MAR is not a safe default

A simulated trial shows why this matters more than it might sound. Dropout is concentrated in the treatment arm and among the sicker treated patients — 30% against 3% in control — so the survivors flatter the drug. Complete-case analysis reports an effect of −1.79 against a true −1.19. Crucially, MAR imputation does not rescue it: at $\delta=0$ the estimate is −1.70, because imputing dropouts to resemble completers repeats exactly the optimistic assumption that caused the problem.

The tipping point

Sweeping $\delta$ gives the tipping point — how much worse the dropouts would have had to be for the conclusion to disappear. Here it is $\delta=4.4$: 2.5 standard deviations of the endpoint, 44% of its observed range, and 3.7 times the effect being estimated. That reads as fairly robust — only a large MNAR effect overturns it. The $\delta$ that would actually recover the truth is 1.8, and the data cannot reveal it, which is precisely what MNAR hides. That sentence, not a single number, is the output of an MNAR analysis, and it is the form regulators ask for in trials with dropout.

A real trial

The real antidepressant trial is smaller and the verdict correspondingly humbler: 28 patients, 18 on drug, with week-6 dropout of 5 in the drug arm against 1 on placebo. The week-6 effect is −1.28 ± 0.70 under MAR — a suggestion of benefit whose interval already includes zero before any MNAR shift is applied. R's mice delta-adjustment reproduces the curve independently at −1.21, CI [−2.64, +0.22]. Any assumption that the drug dropouts fared worse only erodes it further, so the honest reading is promising but unproven — which a complete-case p-value alone would have obscured.

Notebooks

Downloads

Pattern-Mixture Module — Source Code

"""
patternmix.py -- PATTERN-MIXTURE MODELS and delta-sensitivity for MNAR missingness (from scratch).

Backs the notebooks in  "Pattern-Mixture Models -- MNAR by Sensitivity Analysis".

Selection models (Project 4) factor the joint density of outcome y and response indicator r as
p(y) p(r|y) -- a model for the outcome and a model for being observed. PATTERN-MIXTURE models take
the OTHER factorisation,

    p(y, r) = p(y | r) p(r),

stratifying by the missingness PATTERN and giving each stratum its own outcome distribution, then
mixing over the patterns. Its virtue is honesty. The outcome distribution among the MISSING
(r = 0) is never observed, so it cannot be estimated -- it must be ASSUMED. Pattern-mixture makes
that assumption an explicit, interpretable knob: a SENSITIVITY PARAMETER delta linking the missing
stratum to the observed one,

    E[y | missing]  =  E[y | observed, same covariates]  +  delta.

delta = 0 is the MAR ("missing at random", i.e. missing like the observed) assumption; delta != 0 is
MNAR -- the dropouts differ from completers by delta even after adjusting for what we saw. Because
delta is not identified by the data, the analysis is a SENSITIVITY ANALYSIS: sweep delta and watch
the conclusion. The TIPPING POINT is the delta at which the conclusion (e.g. a treatment benefit)
just reverses -- "how much worse would the dropouts have had to be for the result to go away?" This
is exactly the delta-adjustment / reference-based sensitivity analysis regulators ask for in trials
with dropout, and the transparent complement to the selection model's hard-to-interpret rho.

The engine is simple: fit a Bayesian regression to the COMPLETERS, impute each dropout from that
model shifted by delta (applied to a chosen arm), build multiple completed datasets, and combine
the estimand by Rubin's rules -- repeating across a grid of delta.
"""

import numpy as np


def _blr_predict(Xo, yo, Xm, rng):
    """Bayesian linear regression on completers (Xo,yo); return a posterior-predictive draw for Xm."""
    n, p = Xo.shape; XtXi = np.linalg.inv(Xo.T @ Xo + 1e-8 * np.eye(p))
    bhat = XtXi @ (Xo.T @ yo); resid = yo - Xo @ bhat
    s2 = max(resid @ resid, 1e-8) / rng.chisquare(max(n - p, 1))
    beta = bhat + np.linalg.cholesky(s2 * XtXi) @ rng.standard_normal(p)
    return Xm @ beta + np.sqrt(s2) * rng.standard_normal(Xm.shape[0])


def pm_impute(y, arm, Xbase, rng, delta=0.0, arms=(1,), m=50):
    """Pattern-mixture multiple imputation. y:(n,) endpoint with NaN for dropouts; arm:(n,) 0/1;
    Xbase:(n,k) fully observed baseline covariates (an intercept is added). Dropouts are imputed
    from the completer regression (fitted per arm); `delta` is ADDED to the imputed values of the
    dropouts in the arms listed in `arms` (the MNAR shift). Returns m completed y-vectors."""
    y = np.asarray(y, float); arm = np.asarray(arm); Xb = np.column_stack([np.ones(len(y)), Xbase])
    completed = []
    for _ in range(m):
        yc = y.copy()
        for a in (0, 1):
            inarm = arm == a; obs = inarm & ~np.isnan(y); mis = inarm & np.isnan(y)
            if mis.sum() == 0:
                continue
            pred = _blr_predict(Xb[obs], y[obs], Xb[mis], rng)
            if a in arms:
                pred = pred + delta
            yc[mis] = pred
        completed.append(yc)
    return completed


def effect(completed, arm):
    """treatment-arm minus control-arm mean of the completed endpoint, pooled over imputations by
    Rubin's rules. Returns (estimate, se)."""
    arm = np.asarray(arm); ests = np.array([yc[arm == 1].mean() - yc[arm == 0].mean() for yc in completed])
    # within-imputation variance of a difference of means
    within = np.mean([yc[arm == 1].var(ddof=1) / (arm == 1).sum() + yc[arm == 0].var(ddof=1) / (arm == 0).sum()
                      for yc in completed])
    m = len(completed); B = ests.var(ddof=1); T = within + (1 + 1 / m) * B
    return ests.mean(), np.sqrt(T)


def sensitivity(y, arm, Xbase, deltas, rng, arms=(1,), m=50):
    """sweep delta; return arrays of estimate and 95% CI bounds of the treatment effect."""
    est = np.empty(len(deltas)); lo = np.empty(len(deltas)); hi = np.empty(len(deltas))
    for i, dl in enumerate(deltas):
        e, se = effect(pm_impute(y, arm, Xbase, rng, delta=dl, arms=arms, m=m), arm)
        est[i] = e; lo[i] = e - 1.96 * se; hi[i] = e + 1.96 * se
    return est, lo, hi


def tipping_point(deltas, est, lo, hi):
    """smallest |delta| at which the 95% CI first includes 0 (the conclusion becomes non-significant)."""
    sig = (lo > 0) | (hi < 0)
    if sig.all():
        return None
    flip = np.where(~sig)[0]
    return deltas[flip[np.argmin(np.abs(deltas[flip]))]]


def complete_case_effect(y, arm):
    """difference in observed means (dropouts discarded)."""
    y = np.asarray(y, float); arm = np.asarray(arm); o = ~np.isnan(y)
    a1 = o & (arm == 1); a0 = o & (arm == 0)
    est = y[a1].mean() - y[a0].mean()
    se = np.sqrt(y[a1].var(ddof=1) / a1.sum() + y[a0].var(ddof=1) / a0.sum())
    return est, se


def simulate_trial(n, tau, rng, drop_mnar=1.2):
    """two-arm trial with MNAR dropout concentrated in the TREATMENT arm: the sicker (higher-y)
    treated patients drop out, so the observed treatment mean flatters the drug. Returns observed
    endpoint y (NaN=dropout), arm, baseline, and the full (pre-dropout) endpoint."""
    arm = (rng.random(n) < 0.5).astype(int); base = rng.normal(5, 1.5, n)
    yfull = 5 + tau * arm + 0.5 * (base - 5) + rng.normal(0, 1.5, n)          # lower = better; tau<0 = benefit
    pdrop = np.where(arm == 1, 1 / (1 + np.exp(-(yfull - 5) * drop_mnar)), 0.03)   # sicker treated -> drop
    return np.where(rng.random(n) < pdrop, np.nan, yfull), arm, base, yfull

References

• Little, R. J. A. (1993). Pattern-mixture models for multivariate incomplete data. JASA 88(421), 125–134. — the factorisation and its identifying restrictions
• Little, R. J. A. (1994). A class of pattern-mixture models for normal incomplete data. Biometrika 81(3), 471–483. — the normal case implemented here
• National Research Council (2010). The Prevention and Treatment of Missing Data in Clinical Trials. National Academies Press. — why sensitivity analysis, not a single estimate, is the required output
• Carpenter, J. R., Roger, J. H. & Kenward, M. G. (2013). Analysis of longitudinal trials with protocol deviation: a framework for relevant, accessible assumptions, and inference via multiple imputation. Journal of Biopharmaceutical Statistics 23(6), 1352–1371. — reference-based and δ-adjustment sensitivity analysis
• Molenberghs, G. & Kenward, M. G. (2007). Missing Data in Clinical Studies. Wiley. — selection versus pattern-mixture, and why the data cannot choose between them