Missing Covariates in a Regression¶
Impute the predictors — but condition on the outcome¶
The most common missing-data situation in applied work: the outcome $y$ is fully observed but some predictors $x$ have holes. It has two twists that make it worth its own project.
Impute the covariates conditional on the outcome. The Bayesian remedy (Congdon; Ibrahim) is a joint model — a regression for $y$ and a model for the covariates — with the missing predictors augmented inside it: $$y_i=\beta_0+\beta'x_i+\varepsilon_i,\ \varepsilon_i\sim N(0,\sigma^2),\qquad x_i\sim N(\mu,\Sigma).$$ A missing $x_i$ is then imputed from its conditional given the observed covariates and the outcome: combining the Gaussian covariate prior with the Gaussian $y$-likelihood gives a Gaussian draw. The "and the outcome" is essential — impute the predictors from the other predictors only (the frequent mistake — "just fill in the X's") and the $x$–$y$ association is diluted, attenuating $\beta$ toward zero.
Complete-case is unusually robust here. For a missing outcome, listwise deletion is generally biased. For a missing covariate it is unbiased whenever missingness does not depend on $y$ given $x$ — a weaker condition than MCAR — so deletion is more defensible, though still less efficient. We build the joint-model sampler from scratch, demonstrate the attenuation trap and the complete-case result, run the airquality data, and cross-check the covariate imputation in PyMC.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import misscov as MC
rng = np.random.default_rng(3)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("Missing predictors: model the covariates and augment the missing ones INSIDE the regression, imputing them")
print("conditional on the outcome too -- otherwise the coefficient is attenuated toward zero.")
Missing predictors: model the covariates and augment the missing ones INSIDE the regression, imputing them conditional on the outcome too -- otherwise the coefficient is attenuated toward zero.
1. The attenuation trap — impute with the outcome¶
The central demonstration. We simulate $y=1+1.5x_1-0.8x_2+\varepsilon$ with $x_1$ missing (MAR, on the observed $x_2$), and estimate the slope on $x_1$ four ways: the joint model imputing $x_1$ using the outcome; the same sampler imputing $x_1$ from the covariate model only (ignoring $y$); naive mean imputation; and complete-case deletion. Only the methods that respect the $x$–$y$ relationship recover the true slope.
beta=np.array([1.0,1.5,-0.8])
y,X=MC.simulate_misscov(3000,beta,1.0,np.array([0.,0.]),np.array([[1,.5],[.5,1]]),rng)
Xm=X.copy(); pmar=1/(1+np.exp(-X[:,1])); Xm[rng.random(3000)<pmar*0.45/pmar.mean(),0]=np.nan # x1 MAR on observed x2
ry =MC.misscov_gibbs(y,Xm,rng,draws=1500,burn=800,use_y=True)
rny=MC.misscov_gibbs(y,Xm,rng,draws=1500,burn=800,use_y=False)
mb,_=MC.mean_impute(y,Xm); cb,cse,ncc=MC.complete_case(y,Xm)
res={"joint model\n(impute WITH y)":ry["beta"].mean(0)[1], "impute WITHOUT y\n(covariates only)":rny["beta"].mean(0)[1],
"mean\nimputation":mb[1], "complete-case\n(deletion)":cb[1]}
cols=[GREEN,RED,ORANGE,BLUE]
fig,ax=plt.subplots(figsize=(8,3.4))
ax.bar(range(4),list(res.values()),color=cols,alpha=.85)
ax.axhline(1.5,color="k",ls="--",lw=1.5,label="true slope 1.5"); ax.set_xticks(range(4)); ax.set_xticklabels(list(res),fontsize=8)
ax.set_ylabel(r"estimated slope on $x_1$"); ax.set_title("Impute covariates WITH the outcome, or attenuate the slope")
for i,v in enumerate(res.values()): ax.text(i,v+0.03,f"{v:.2f}",ha="center",fontsize=9)
ax.legend(frameon=False,fontsize=8); plt.tight_layout(); plt.show()
print(f"joint model {ry['beta'].mean(0)[1]:.2f} and complete-case {cb[1]:.2f} recover 1.5; imputing the predictor")
print(f"from the other predictors only collapses it to {rny['beta'].mean(0)[1]:.2f}, and mean imputation to {mb[1]:.2f}. Filling")
print("in covariates without using the outcome dilutes exactly the relationship you are trying to estimate.")
joint model 1.50 and complete-case 1.51 recover 1.5; imputing the predictor from the other predictors only collapses it to 0.84, and mean imputation to 1.33. Filling in covariates without using the outcome dilutes exactly the relationship you are trying to estimate.
2. Why "without $y$" fails, and the joint model works¶
The mechanism is easy to see. Imputing $x_1$ from $x_2$ alone places every filled value on the $x_2$-regression line, so the imputed points carry no residual $x_1$ variation tied to $y$ — the extra information $x_1$ has about $y$ is erased, and the slope shrinks. The joint model instead pulls each imputed $x_1$ toward values consistent with the $y$ actually observed, preserving the association. The scatter below shows imputed $x_1$ values against $y$: without the outcome they ignore $y$; with it they track it.
mask=np.isnan(Xm[:,0])
fig,ax=plt.subplots(1,2,figsize=(12,4),sharey=True)
ax[0].scatter(y[~mask],X[~mask,0],s=8,color=GREY,alpha=.4,label="observed")
ax[0].scatter(y[mask],rny["Ximp"][mask,0],s=10,color=RED,alpha=.6,label="imputed (no y)")
ax[0].set_title("Imputed WITHOUT the outcome — flat in y"); ax[0].set_xlabel("outcome y"); ax[0].set_ylabel(r"$x_1$"); ax[0].legend(frameon=False,fontsize=8)
ax[1].scatter(y[~mask],X[~mask,0],s=8,color=GREY,alpha=.4,label="observed")
ax[1].scatter(y[mask],ry["Ximp"][mask,0],s=10,color=GREEN,alpha=.6,label="imputed (with y)")
ax[1].set_title("Imputed WITH the outcome — tracks y"); ax[1].set_xlabel("outcome y"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("Left: the outcome-blind imputations are flat across y, erasing x1's information about y -> attenuation.")
print("Right: the joint model's imputations rise with y, matching the observed cloud -> the slope is preserved.")
Left: the outcome-blind imputations are flat across y, erasing x1's information about y -> attenuation. Right: the joint model's imputations rise with y, matching the observed cloud -> the slope is preserved.
3. When is complete-case deletion safe? — a covariate-specific result¶
A result that surprises people: for a missing covariate, complete-case regression is unbiased whenever the chance of being missing does not depend on the outcome given the covariates — even if it depends on the covariates themselves. That is why deletion did fine above (missingness was driven by $x_2$). It fails only when missingness depends on the outcome. The two panels contrast the cases; the joint model is correct in both, deletion only in the first.
fig,ax=plt.subplots(1,2,figsize=(12,3.6))
# (a) missingness depends on a covariate -> complete-case OK
yA,XA=MC.simulate_misscov(4000,beta,1.0,np.array([0.,0.]),np.array([[1,.5],[.5,1]]),rng)
XmA=XA.copy(); p=1/(1+np.exp(-XA[:,1])); XmA[rng.random(4000)<p*0.45/p.mean(),0]=np.nan
cbA,_,_=MC.complete_case(yA,XmA); rjA=MC.misscov_gibbs(yA,XmA,rng,draws=1000,burn=600)
# (b) missingness depends on the outcome y (with noisy y) -> complete-case biased
yB,XB=MC.simulate_misscov(4000,beta,3.0,np.array([0.,0.]),np.array([[1,.5],[.5,1]]),rng)
XmB=XB.copy(); q=1/(1+np.exp(-yB/3)); XmB[rng.random(4000)<q*0.5/q.mean(),0]=np.nan
cbB,_,_=MC.complete_case(yB,XmB); rjB=MC.misscov_gibbs(yB,XmB,rng,draws=1000,burn=600)
for k,(ttl,cc,jm) in enumerate([("missingness depends on a COVARIATE",cbA[1],rjA["beta"].mean(0)[1]),
("missingness depends on the OUTCOME",cbB[1],rjB["beta"].mean(0)[1])]):
ax[k].bar(["complete\ncase","joint\nmodel"],[cc,jm],color=[BLUE,GREEN],alpha=.85)
ax[k].axhline(1.5,color="k",ls="--",lw=1.2); ax[k].set_ylim(0,1.8); ax[k].set_title(ttl,fontsize=9)
for i,v in enumerate([cc,jm]): ax[k].text(i,v+0.05,f"{v:.2f}",ha="center")
ax[0].set_ylabel(r"slope on $x_1$ (truth 1.5)"); plt.tight_layout(); plt.show()
print("Left: missingness on a covariate -> complete-case is unbiased (and so is the joint model).")
print("Right: missingness on the outcome -> complete-case is biased, but the joint model still recovers the slope.")
print("So for MISSING COVARIATES deletion is more defensible than usual -- but imputation is both safe and efficient.")
Left: missingness on a covariate -> complete-case is unbiased (and so is the joint model). Right: missingness on the outcome -> complete-case is biased, but the joint model still recovers the slope. So for MISSING COVARIATES deletion is more defensible than usual -- but imputation is both safe and efficient.
4. Real data — airquality, with missing predictors¶
We reuse the airquality data from Project 1, but reframed as a regression: model Temp (fully observed) on Ozone, Solar.R and Wind — where Ozone (37 missing) and Solar.R (7 missing) are the incomplete covariates (associational; warm days tend to be sunny, high-ozone, low-wind). Project 1 imputed all four variables jointly; here three of them are predictors of the fourth, and we impute the missing predictors conditional on Temp and each other.
aq=pd.read_csv("airquality.csv")
Xr=np.column_stack([np.log(aq["Ozone"]), np.log(aq["Solar.R"]), aq["Wind"]]) # missing covariates: logOzone, logSolar
yr=aq["Temp"].to_numpy(float)
r=MC.misscov_gibbs(yr,Xr,rng,draws=3000,burn=1500); cb,cse,ncc=MC.complete_case(yr,Xr)
names=["intercept","log Ozone","log Solar.R","Wind"]; bm=r["beta"].mean(0); bl,bh=np.percentile(r["beta"],[2.5,97.5],axis=0)
print(f"Temp ~ log Ozone + log Solar.R + Wind (joint model uses all {len(yr)} rows; complete-case only {ncc})\n")
print(f"{'term':13s}{'joint est':>10s}{'95% CrI':>18s} {'complete-case':>14s}")
for i,nm in enumerate(names):
print(f"{nm:13s}{bm[i]:10.2f} [{bl[i]:6.2f},{bh[i]:6.2f}] {cb[i]:14.2f}")
fig,ax=plt.subplots(figsize=(7,3)); idx=[1,2,3]
ax.errorbar(bm[idx],np.arange(3)+.1,xerr=[bm[idx]-bl[idx],bh[idx]-bm[idx]],fmt="o",color=GREEN,capsize=3,label=f"joint model (n={len(yr)})")
ax.scatter(cb[idx],np.arange(3)-.1,color=BLUE,marker="s",label=f"complete-case (n={ncc})")
ax.axvline(0,color="k",lw=.8); ax.set_yticks(range(3)); ax.set_yticklabels(names[1:]); ax.set_xlabel("coefficient on Temp")
ax.set_title("airquality: warm days are high-ozone, sunny, low-wind"); ax.legend(frameon=False,fontsize=8); plt.tight_layout(); plt.show()
print("Ozone carries the strong positive association with temperature and wind a negative one; solar radiation adds")
print("little once ozone is in the model (its interval spans zero). The joint model keeps all 153 days with credible")
print(f"intervals, where complete-case would discard the {len(yr)-ncc} days that have a missing predictor.")
Temp ~ log Ozone + log Solar.R + Wind (joint model uses all 153 rows; complete-case only 111) term joint est 95% CrI complete-case intercept 55.38 [ 46.53, 63.65] 55.71 log Ozone 7.42 [ 5.56, 9.25] 7.51 log Solar.R 0.09 [ -1.57, 1.83] -0.08 Wind -0.33 [ -0.69, 0.04] -0.32
Ozone carries the strong positive association with temperature and wind a negative one; solar radiation adds little once ozone is in the model (its interval spans zero). The joint model keeps all 153 days with credible intervals, where complete-case would discard the 42 days that have a missing predictor.
5. Cross-check in PyMC — automatic covariate imputation¶
PyMC imputes missing covariates the same joint way: model the covariate vector as a masked multivariate normal (so the missing entries become latent variables) and use it in the regression, so the outcome informs the imputations automatically. This is the pattern in PyMC's missing values in covariates example. We compare its coefficients with the from-scratch joint model.
import pymc as pm
Xr_m=np.ma.masked_invalid(Xr)
with pm.Model() as mod:
mx=pm.Normal("mx",[3,5,10],5,shape=3)
chol,_,_=pm.LKJCholeskyCov("C",n=3,eta=2,sd_dist=pm.HalfNormal.dist(5),compute_corr=True)
Xlat=pm.MvNormal("X",mu=mx,chol=chol,observed=Xr_m) # imputes missing Ozone/Solar covariates
a0=pm.Normal("a0",70,20); b=pm.Normal("b",0,10,shape=3); s=pm.HalfNormal("s",10)
pm.Normal("Temp", a0+Xlat@b, s, observed=yr) # outcome informs the imputations
idata=pm.sample(800,tune=1200,chains=4,target_accept=0.9,random_seed=5,progressbar=False)
b_pm=idata.posterior["b"].mean(("chain","draw")).values
print("coefficient from-scratch PyMC")
for i,nm in enumerate(names[1:]): print(f"{nm:13s}{bm[i+1]:10.2f} {b_pm[i]:6.2f}")
print(f"\ncorrelation of the three coefficients: {np.corrcoef(bm[1:],b_pm)[0,1]:.3f}")
print("PyMC's masked-covariate model reproduces the from-scratch joint fit: making the missing predictors latent")
print("variables inside the regression is the automatic version of imputing them conditional on the outcome.")
g++ not available, if using conda: `conda install gxx`
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pymc\model\core.py:1337: ImputationWarning: Data in X contains missing values and will be automatically imputed from the sampling distribution. warnings.warn(impute_message, ImputationWarning)
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [mx, C, X_unobserved, a0, b, s]
Sampling 4 chains for 1_200 tune and 800 draw iterations (4_800 + 3_200 draws total) took 10 seconds.
coefficient from-scratch PyMC log Ozone 7.42 7.31 log Solar.R 0.09 0.03 Wind -0.33 -0.36 correlation of the three coefficients: 1.000 PyMC's masked-covariate model reproduces the from-scratch joint fit: making the missing predictors latent variables inside the regression is the automatic version of imputing them conditional on the outcome.
6. Summary¶
Missing predictors are the everyday case, and the joint model handles them by augmenting the missing covariates inside the regression — a model for $y$ given $x$ paired with a model for $x$. The decisive detail is that each missing covariate is imputed conditional on the outcome as well as the other covariates: do that and $\beta$ is recovered; impute the predictors from the other predictors alone (or by mean substitution) and the coefficient is attenuated toward zero, because the imputations erase exactly the $x$–$y$ information being estimated. A useful covariate-specific result balances this: complete-case deletion is unbiased for $\beta$ unless missingness depends on the outcome — more forgiving than for a missing outcome, though it still throws away data, as the 40 discarded airquality days showed. A masked-covariate PyMC model reproduced the from-scratch joint fit automatically.
The connections: this is Congdon's missing-covariate model (BMCD 11.3) and the PyMC missing values in covariates example; the "impute $x$ using $y$" rule is exactly why MICE (Project 2) puts the outcome in every imputation model, and the joint covariate model is the multivariate-normal augmentation of Project 1 now doing double duty. It also sits next to errors-in-variables / measurement-error models, where a covariate is likewise a partially known latent quantity. All of this remains within MAR; when missingness depends on the unobserved value itself, the next projects — selection (Project 4) and pattern-mixture (Project 5) models — take over.