Missing Covariates in a Regression

Impute the Predictor Using the Outcome

The commonest missing-data situation in practice is a fully observed outcome and predictors with holes in them. The Bayesian remedy is a joint model: write a distribution for the covariates, and augment the missing ones inside the regression. Gibbs then alternates three blocks — draw the regression coefficients from the completed data, draw the covariate moments from the completed covariates, and impute each missing predictor from its full conditional given the observed predictors and the outcome.

That last clause is the whole project. A missing $x_i$ is informed by two things: the covariate model says it should look like the other predictors, and the regression says it must be consistent with the $y$ that was actually observed. Drop the second and you dilute precisely the association you are trying to estimate. The module exposes a use_y switch so the failure can be produced on demand rather than described.

The attenuation, produced on demand

And it is dramatic. On simulated data with a true coefficient of 1.5, the joint model recovers 1.50 and complete-case deletion 1.51 — while imputing the predictor from the other predictors alone collapses it to 0.84, and mean imputation lands at 1.33. R's mice reproduces the same failure on the real data: the ozone–temperature coefficient falls from 7.65 to 5.39 when the outcome is removed from the imputation model. The diagnostic plot makes the mechanism visible — outcome-blind imputations are flat across $y$, erasing the very information the slope depends on, while the joint model's imputations rise with $y$ and match the observed cloud.

When deletion is actually safe

The project also makes a point that runs against the usual advice, and makes it correctly. For missing covariates, complete-case deletion is unbiased whenever missingness does not depend on $y$ given $x$ — a weaker condition than MCAR — which is why it recovers 1.51 above. The notebook shows the contrast directly: with missingness on a covariate, deletion is fine; with missingness on the outcome, it is biased and only the joint model recovers the slope. So here deletion is more defensible than usual — just less efficient, discarding 42 of 153 days on the real data.

Real data

On `airquality` — modelling temperature on log ozone, log solar radiation and wind — the joint model keeps all 153 days with credible intervals where complete-case would keep 111. Ozone carries a strong positive association (7.42, CrI [5.56, 9.25]), wind a negative one, and solar radiation adds little once ozone is in the model, its interval spanning zero. R's mice gives 7.65 on the same log scale, and a PyMC model with the missing predictors as latent variables reproduces the from-scratch fit at correlation 1.000 — masked-array imputation is the automatic version of imputing conditional on the outcome.

Notebooks

Downloads

Missing-Covariate Module — Source Code

"""
misscov.py -- MISSING COVARIATES in a regression, the joint-model way (from scratch).

Backs the notebooks in  "Missing Covariates in a Regression".

The commonest missing-data situation is a fully observed OUTCOME y and PREDICTORS x with holes. The
Bayesian remedy (Congdon; Ibrahim et al.) is to write a model for the covariates and augment the
missing ones inside the regression -- a JOINT model of outcome and covariates:

    y_i = beta0 + beta' x_i + eps_i,   eps_i ~ N(0, sigma^2)          (the regression of interest)
    x_i ~ N(mu, Sigma)                                                (a model for the covariates)

Gibbs alternates three blocks: draw (beta, sigma^2) from the completed regression, draw (mu, Sigma)
from the completed covariates (conjugate NIW), and IMPUTE each missing covariate from its full
conditional given the observed covariates AND the outcome.

That last "and the outcome" is the crux of missing-covariate imputation. A missing x_i is informed by
two things -- the covariate model (x correlates with the observed covariates) and the regression
(x must be consistent with the y actually seen). Combining a Gaussian covariate prior with the
Gaussian y-likelihood gives a Gaussian conditional, and the y-term is what keeps the estimated
beta unbiased. Impute the covariates WITHOUT the outcome (the common mistake -- filling predictors
from the other predictors only) and the x-y association is diluted, ATTENUATING beta toward zero.
The `use_y` switch below turns that term on and off so the attenuation can be shown directly.

A second lesson specific to missing covariates: complete-case regression is unbiased for beta
whenever missingness does not depend on y given x (a weaker condition than MCAR) -- so here deletion
is more defensible than for missing outcomes, though still less efficient than imputation.
"""

import numpy as np
from scipy.stats import invwishart


def simulate_misscov(N, beta, sigma, mu_x, Sig_x, rng):
    """y = beta0 + beta[1:]'x + N(0,sigma^2);  x ~ N(mu_x, Sig_x)."""
    X = rng.multivariate_normal(np.asarray(mu_x, float), np.asarray(Sig_x, float), N)
    y = beta[0] + X @ np.asarray(beta[1:], float) + sigma * rng.standard_normal(N)
    return y, X


def _patterns(Mask):
    g = {}
    for i, row in enumerate(Mask):
        g.setdefault(tuple(np.where(row)[0]), []).append(i)
    return {k: np.array(v) for k, v in g.items()}


def misscov_gibbs(y, X, rng, draws=2000, burn=1000, use_y=True):
    """Joint-model Gibbs for a regression with missing covariates. X:(N,p) with NaN, y:(N,) observed.
    If use_y=False the missing covariates are imputed from the covariate model ONLY (ignoring the
    outcome) -- which attenuates beta. Returns posterior draws of beta:(draws,p+1), sigma^2, the
    covariate moments, and the posterior-mean imputed X."""
    y = np.asarray(y, float); X = np.array(X, float); N, p = X.shape
    Mask = np.isnan(X); pat = _patterns(Mask)
    Xc = X.copy()
    for j in range(p):
        Xc[Mask[:, j], j] = np.nanmean(X[:, j])
    beta = np.zeros(p + 1); sig2 = np.var(y)
    B = np.empty((draws, p + 1)); S2 = np.empty(draws); MU = np.empty((draws, p)); imp_sum = np.zeros((N, p)); nimp = 0
    for it in range(draws + burn):
        # ---- (beta, sigma^2) | completed data ----
        D = np.column_stack([np.ones(N), Xc]); DtDi = np.linalg.inv(D.T @ D + 1e-8 * np.eye(p + 1))
        bhat = DtDi @ (D.T @ y); resid = y - D @ bhat
        sig2 = (resid @ resid) / rng.chisquare(max(N - (p + 1), 1))
        beta = bhat + np.linalg.cholesky(sig2 * DtDi) @ rng.standard_normal(p + 1)
        # ---- (mu, Sigma) | completed covariates  (noninformative NIW) ----
        xbar = Xc.mean(0); Sx = (Xc - xbar).T @ (Xc - xbar)
        Sigma = np.atleast_2d(invwishart.rvs(df=N - 1, scale=Sx, random_state=rng))
        mu = xbar + np.linalg.cholesky(Sigma / N) @ rng.standard_normal(p)
        # ---- impute missing covariates | observed covariates (+ outcome) ----
        b0 = beta[0]; bc = beta[1:]
        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:
                Sooi = np.linalg.inv(Sigma[np.ix_(obs, obs)]); Bm = Sigma[np.ix_(mis, obs)] @ Sooi
                a = mu[mis] + (Xc[np.ix_(rows, obs)] - mu[obs]) @ Bm.T          # covariate-model mean
                V = Sigma[np.ix_(mis, mis)] - Bm @ Sigma[np.ix_(obs, mis)]
            else:
                a = np.tile(mu[mis], (len(rows), 1)); V = Sigma[np.ix_(mis, mis)]
            Vi = np.linalg.inv(V); bm = bc[mis]
            if use_y:
                r = y[rows] - b0 - (Xc[np.ix_(rows, obs)] @ bc[obs] if obs else 0.0)   # y residual
                P = Vi + np.outer(bm, bm) / sig2                              # combine prior + y-likelihood
                Pi = np.linalg.inv(P)
                rhs = a @ Vi.T + np.outer(r, bm) / sig2
                postmean = rhs @ Pi.T; Lp = np.linalg.cholesky(Pi)
            else:                                                            # ignore the outcome
                postmean = a; Lp = np.linalg.cholesky(V)
            Xc[np.ix_(rows, mis)] = postmean + rng.standard_normal((len(rows), len(mis))) @ Lp.T
        if it >= burn:
            i = it - burn; B[i] = beta; S2[i] = sig2; MU[i] = mu; imp_sum += Xc; nimp += 1
    return dict(beta=B, sigma2=S2, mu=MU, Ximp=imp_sum / nimp, mask=Mask)


def complete_case(y, X):
    """listwise-deletion OLS -- unbiased for beta if missingness is independent of y given x."""
    y = np.asarray(y, float); X = np.asarray(X, float)
    keep = ~np.isnan(X).any(axis=1); D = np.column_stack([np.ones(keep.sum()), X[keep]])
    DtDi = np.linalg.inv(D.T @ D); b = DtDi @ (D.T @ y[keep])
    resid = y[keep] - D @ b; s2 = (resid @ resid) / (keep.sum() - D.shape[1])
    return b, np.sqrt(np.diag(s2 * DtDi)), keep.sum()


def mean_impute(y, X):
    """fill each missing covariate with its column mean, then OLS -- biased (attenuates beta)."""
    y = np.asarray(y, float); Xc = np.array(X, float)
    for j in range(Xc.shape[1]):
        Xc[np.isnan(Xc[:, j]), j] = np.nanmean(Xc[:, j])
    D = np.column_stack([np.ones(len(y)), Xc]); DtDi = np.linalg.inv(D.T @ D); b = DtDi @ (D.T @ y)
    resid = y - D @ b; s2 = (resid @ resid) / (len(y) - D.shape[1])
    return b, np.sqrt(np.diag(s2 * DtDi))

References

• Ibrahim, J. G., Chen, MH., Lipsitz, S. R. & Herring, A. H. (2005). Missing-data methods for generalized linear models: a comparative review. JASA 100(469), 332–346. — the joint-model approach to missing covariates
• Little, R. J. A. (1992). Regression with missing X's: a review. JASA 87(420), 1227–1237. — including the condition under which complete-case regression stays unbiased
• Moons, K. G. M., Donders, R. A. R. T., Stijnen, T. & Harrell, F. E. (2006). Using the outcome for imputation of missing predictor values was preferred. Journal of Clinical Epidemiology 59(10), 1092–1101. — the attenuation this project demonstrates, argued in the applied literature
• Congdon, P. (2010). Applied Bayesian Hierarchical Methods. CRC Press. — augmenting missing covariates within the regression
• van Buuren, S. (2018). Flexible Imputation of Missing Data (2nd ed.). CRC Press. — why the outcome belongs in the imputation model, and mice's defaults