5.3 Missing data

What to do about missing values depends on why they are missing. If they are missing at random given what we did observe, multiple imputation fills them in honestly, with their uncertainty. If they are missing because of the values themselves, no method can tell from the data alone, and only a sensitivity analysis can show how much it matters.

Code: data and style used in this chapter

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes

from statly import style

style.use()
rng = np.random.default_rng(1)

# Simulated patients: age (always recorded) and cholesterol (sometimes missing).
N = 500
age = rng.normal(50, 12, N)
cholesterol = 120 + 1.6 * age + rng.normal(0, 18, N)
TRUE_MEAN = cholesterol.mean()


def scatter_missing(ax: Axes, missing: np.ndarray, title: str) -> None:
    """Observed cholesterol (blue) and the values that went missing (open circles)."""
    ax.scatter(age[~missing], cholesterol[~missing], s=7, color=style.BLUE, alpha=0.6)
    ax.scatter(
        age[missing],
        cholesterol[missing],
        s=9,
        facecolors="none",
        edgecolors=style.CORAL,
        alpha=0.6,
    )
    gap = cholesterol[~missing].mean() - TRUE_MEAN
    style.note(
        ax,
        f"{missing.mean():.0%} missing\nmean of observed: {gap:+.1f}\nfrom the truth",
        x=0.03,
        y=0.97,
        color=style.INK,
    )
    ax.set(title=title, xlabel="age", ylabel="cholesterol")

5.3.1 Why values are missing

Three situations, with names that are easy to mix up:

Missing completely at random (MCAR): whether a value is missing has nothing to do with any value, observed or not. A test tube is dropped.
Missing at random (MAR): whether a value is missing depends only on things we did observe. Older patients are tested more often, and age is recorded for all.
Missing not at random (MNAR): whether a value is missing depends on the missing value itself. Patients with high cholesterol skip the follow-up because they feel unwell.

Code

def logistic(x: np.ndarray) -> np.ndarray:
    """The S-shaped curve from any number to a probability."""
    return 1 / (1 + np.exp(-x))


MISSING = {
    "MCAR": rng.uniform(size=N) < 0.3,
    "MAR (depends on age)": rng.uniform(size=N) < logistic(-0.9 - 0.12 * (age - 50)),
    "MNAR (depends on itself)": rng.uniform(size=N)
    < logistic(-0.9 + 0.06 * (cholesterol - cholesterol.mean())),
}

fig, axes = style.subplots(3, height=3.1, width=10.5)
for ax, (name, missing) in zip(axes, MISSING.items()):
    scatter_missing(ax, missing, name)
fig.tight_layout()
plt.show()

Figure 5.3.1: The same 500 simulated patients, with cholesterol missing for roughly a third (open red circles) under each mechanism. MCAR: the observed values are a random subset, and their mean is right. MAR: the young are more often missing; the observed mean is too high, but age explains why, and the relation with age is intact. MNAR: high values are more often missing; the observed mean is too low, and nothing observed reveals it.

No test can tell MAR from MNAR, because the difference lies in values we do not have. The choice rests on knowledge of how the data were collected.

5.3.2 Simple fixes and what they cost

Complete-case analysis drops every row with a missing value. Unbiased under MCAR (just wasteful); biased under MAR for many questions, and under MNAR.
Mean imputation fills each gap with the column mean. It keeps the mean (under MCAR) but shrinks the spread and weakens every relation with other variables.
Regression imputation fills each gap with its prediction from other variables. Better, but the filled-in values lie exactly on the line, so the data look more certain than they are.
Stochastic regression imputation adds random noise of the right size to each prediction, restoring the spread.

Code

missing = MISSING["MAR (depends on age)"]
fit = sm.OLS(cholesterol[~missing], sm.add_constant(age[~missing])).fit()
predicted = fit.params[0] + fit.params[1] * age
residual_sd = np.sqrt(fit.scale)
FILLS = {
    "mean imputation": np.full(N, cholesterol[~missing].mean()),
    "regression imputation": predicted,
    "regression + noise": predicted + rng.normal(0, residual_sd, N),
}

fig, axes = style.subplots(3, height=3.1, width=10.5)
for ax, (name, fill) in zip(axes, FILLS.items()):
    completed = np.where(missing, fill, cholesterol)
    ax.scatter(age[~missing], cholesterol[~missing], s=7, color=style.BLUE, alpha=0.5)
    ax.scatter(age[missing], fill[missing], s=9, color=style.CORAL, alpha=0.8)
    slope = np.polyfit(age, completed, 1)[0]
    style.note(
        ax,
        f"slope {slope:.2f} (true 1.6)\nsd {completed.std():.0f} (true {cholesterol.std():.0f})",
        x=0.03,
        y=0.97,
        color=style.INK,
    )
    ax.set(title=name, xlabel="age", ylabel="cholesterol")
fig.tight_layout()
plt.show()

Figure 5.3.2: The MAR data of Figure 2.1.3.1, filled in three ways (filled-in values in red). Mean imputation puts all of them on one horizontal line and flattens the relation with age. Regression imputation puts them on the regression line: the relation is right, the scatter is too small. Adding noise to the predictions gives filled-in values that look like real ones.

5.3.3 Multiple imputation

One filled-in data set, however good, treats guesses as if they were measurements, so standard errors come out too small. Multiple imputation fills the gaps m times (often 20 or more), each time with different random draws, including draws of the imputation model’s own coefficients. Each completed data set is analysed as usual, and the m results are pooled with Rubin’s rules: the estimate is the average, and its variance adds the average within-imputation variance $\bar W$ to the spread between imputations B:

$\hat\theta = \frac{1}{m}\sum_{j=1}^m \hat\theta_j, \qquad T = \bar W + \Big(1 + \frac{1}{m}\Big) B$

The between-imputation part is exactly the uncertainty caused by the missing values.

Code

M = 20


def one_imputation(seed: int) -> tuple[float, float]:
    """Mean and its variance from one completed data set, with coefficients and noise drawn at random."""
    local = np.random.default_rng(seed)
    coef = local.multivariate_normal(
        fit.params, fit.cov_params()
    )  # uncertainty in the imputation model
    fill = coef[0] + coef[1] * age + local.normal(0, residual_sd, N)
    completed = np.where(missing, fill, cholesterol)
    return completed.mean(), completed.var(ddof=1) / N


draws = np.array([one_imputation(s) for s in range(M)])
pooled = draws[:, 0].mean()
total_var = draws[:, 1].mean() + (1 + 1 / M) * draws[:, 0].var(ddof=1)
single = np.where(missing, predicted, cholesterol)
observed = cholesterol[~missing]
RESULTS = {
    "complete cases": (observed.mean(), observed.std(ddof=1) / np.sqrt(len(observed))),
    "single imputation": (single.mean(), single.std(ddof=1) / np.sqrt(N)),
    "multiple imputation": (pooled, np.sqrt(total_var)),
}

fig, axes = style.subplots(2, height=3.1, width=10)
axes[0].scatter(range(1, M + 1), draws[:, 0], color=style.GREY, s=18)
axes[0].axhline(pooled, color=style.BLUE, lw=2)
axes[0].axhspan(
    pooled - 1.96 * np.sqrt(total_var),
    pooled + 1.96 * np.sqrt(total_var),
    color=style.BLUE,
    alpha=0.15,
)
axes[0].axhline(TRUE_MEAN, color=style.CORAL, ls="--", lw=1.2)
axes[0].set(xlabel="imputation", ylabel="mean cholesterol", title=f"{M} imputations")
for i, (name, (est, se)) in enumerate(RESULTS.items()):
    axes[1].errorbar(
        est,
        i,
        xerr=1.96 * se,
        fmt="o",
        color=(style.GREY, style.GOLD, style.BLUE)[i],
        capsize=4,
        lw=2,
    )
axes[1].axvline(TRUE_MEAN, color=style.CORAL, ls="--", lw=1.2)
axes[1].set_yticks(range(3), list(RESULTS))
axes[1].invert_yaxis()
axes[1].set(xlabel="mean cholesterol, 95% interval", title="Three answers")
fig.tight_layout()
plt.show()

Figure 5.3.3: Estimating mean cholesterol from the MAR data. Left: the 20 imputed means (dots), each from a different completed data set, with their pooled estimate and 95% interval (blue). Right: the estimates with 95% intervals from complete cases (biased upwards), single regression imputation (right centre, but too narrow), and multiple imputation (right centre, honest width). Dashed: the true mean.

With several incomplete variables, chained equations (MICE) impute each variable in turn from all the others, cycling until the imputations settle. The imputation model should include every variable of the analysis, the outcome too, and any variable that predicts missingness.

5.3.4 When data are missing not at random: sensitivity analysis

Under MNAR, imputing from observed variables still leaves a bias, and the data cannot say how large. A sensitivity analysis asks: if the missing values were systematically δ higher than the imputation model predicts, how would the answer change? If the conclusion survives every plausible δ, it is robust; if not, the honest answer is that the data cannot settle the question.

Code

missing_nr = MISSING["MNAR (depends on itself)"]
fit_nr = sm.OLS(cholesterol[~missing_nr], sm.add_constant(age[~missing_nr])).fit()
fill_nr = fit_nr.params[0] + fit_nr.params[1] * age
deltas = np.linspace(0, 40, 21)
means = np.array([np.where(missing_nr, fill_nr + d, cholesterol).mean() for d in deltas])
se = cholesterol[~missing_nr].std() / np.sqrt(N)

fig, ax = plt.subplots(figsize=(8, 3.2))
ax.fill_between(deltas, means - 1.96 * se, means + 1.96 * se, color=style.BLUE, alpha=0.2)
ax.plot(deltas, means, color=style.BLUE, lw=2)
ax.axhline(TRUE_MEAN, color=style.CORAL, ls="--", lw=1.2)
ax.set(xlabel="assumed shift δ of the missing values", ylabel="estimated mean cholesterol")
plt.show()

Figure 5.3.4: The MNAR data of Figure 2.1.3.1: imputed from age, then every imputed value shifted by δ. The estimated mean (with its 95% interval) rises with δ. Here the truth (dashed) corresponds to δ of about 16; the data alone give no way to know that, so the analysis reports the range of answers for the range of δ thought plausible.

5.3.5 Missing values in prediction

For prediction (Part 6), the aim is accuracy, not unbiased coefficients, and the fact that a value is missing can itself be informative. Common practice: impute inside the cross-validation folds (Chapter 6.2), and add a missingness indicator, a column that is 1 where the value was missing. Many tree methods handle missing values directly. Either way, the imputation must be learned from training data only.

5.3.6 In one picture

MCAR: missingness unrelated to anything; MAR: related to observed values only; MNAR: related to the missing values themselves. MAR and MNAR cannot be told apart from the data.
Complete cases are biased unless MCAR; mean imputation distorts spreads and relations; single imputation is overconfident.
Multiple imputation + Rubin’s rules carry the uncertainty of the missing values into the standard errors.
Under possible MNAR, run a sensitivity analysis over assumed shifts.
In prediction, impute within folds and consider a missingness indicator.

Next: cause and effect, and how to estimate it when an experiment is not possible.