Why Predictive Model Selection Fails for Causal Effect Estimation

A Question That Sounds Easy But Isn’t

The usual framing sets the wrong goal. When estimating the effect of an exposure on an outcome — say, air pollution on hospital admissions, or a drug versus time to recovery — we have observational data, one exposure variable, and a pile of candidate covariates. Some are confounders. Some are noise. The biggest question is: which ones go in the model?

The usual reflex hands the choice to a selection rule. Run stepwise routines, compare BIC values, or cross-validate Lasso. Full Bayesian averaging proceeds over candidate specifications.

The problem with these approaches is that every one of these tools optimizes predictive fit. A model can predict the outcome beautifully while being systematically wrong about the effect we care about.

This issue was the subject of a 2012 Biometrics paper by Chi Wang, Giovanni Parmigiani, and Francesca Dominici, who called the underlying issue adjustment uncertainty — the uncertainty about which variables to adjust for — and proposed a fix called Bayesian Adjustment for Confounding (BAC) [1]. It also has a famous frequentist twin that most people don’t realize is the same insight: post-double-selection and double machine learning [3][4].

This post demonstrates three things: the failure mode in runnable Python, the fix, and the fix’s own sharp edge — the case where “adjust for anything that predicts exposure” is exactly the wrong advice and makes your estimate worse.

If you have been reading about propensity scores, subclassification, and matching, this is the question those methods assume you have already answered: how do you decide, from data, which covariates a valid adjustment requires?

A Good Predictor and a Confounder Are Not the Same

Why does a specification tuned for prediction produce a flawed effect estimate? Because the two tasks reward distinct traits. Prediction favors a covariate that accounts for outcome variance. Effect estimation favors a covariate tied to both exposure and outcome.

The problem arises when a confounder is strongly associated with the exposure but only weakly with the outcome. To a predictive criterion, that variable looks nearly useless because it explains little outcome variance, and much of what it explains is already captured by the correlated exposure. In realistic finite samples, BIC may prefer dropping it. But dropping it changes what the exposure coefficient means. The coefficient now absorbs the confounder’s path, and the estimated effect is contaminated.

Wang and coauthors illustrate the point with five candidate covariates. U1, U2, and U3 each link to both X and Y. U4 predicts Y alone. U5 is noise.

Figure 1: The setup from Wang et al. (2012), redrawn. U1, U2, and U3 are true confounders. Each touches both X and Y. U4 predicts only Y; U5 is pure noise. The trap is U2: strong on exposure, weak on outcome.

Consider the argument reduced to two t-statistics. On the simulated data below, U2 shows a t-statistic of 1.2 inside the outcome model, yielding p = 0.23 — indistinguishable from noise. Inside the exposure model, the same variable shows a t-statistic of 31.5. Same variable, same dataset. One model calls it irrelevant. The other calls it central.

Every routine we normally invoke examines only the outcome model. Omitted variable bias from excluding U2 equals:

It is the product of U2’s coefficient in the outcome equation times its coefficient in the regression of X on U2 (i.e., delta). The first factor stays small — which is exactly why BIC remains unmoved. The second factor grows large because U2 tracks exposure closely. Their product contaminates the target estimate. No outcome fit rule ever inspects the second factor.

BIC remains selection consistent and would retain U2 with sufficient observations. The difficulty is finite samples. The asymptotic disappearance of bias offers little comfort when the dataset contains roughly 1,000 rows.

The Failure Mode, in Python

Simulating the setup in Figure 1 and fitting two outcome models — one fully adjusted, one that drops U2 — makes the problem concrete.

import numpy as np, pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(11)
n, M = 1000, 5

U = rng.normal(size=(n, M))                              # candidates U1..U5
X = U[:,0] + U[:,1] + 0.1*U[:,2] + rng.normal(size=n)   # exposure
Y = 0.1*X + U[:,0] + 0.1*U[:,1] + U[:,2] + U[:,3] \
    + rng.normal(size=n)                                 # true effect = 0.1

cols = [f"U{j+1}" for j in range(M)]
df = pd.DataFrame(U, columns=cols); df["X"] = X; df["Y"] = Y

def fit_outcome(includes):
    feats = ["X"] + [c for c, inc in zip(cols, includes) if inc]
    return sm.OLS(df["Y"], sm.add_constant(df[feats])).fit()

full  = fit_outcome((1,1,1,1,0))   # adjusts for U1..U4
noU2  = fit_outcome((1,0,1,1,0))   # silently drops confounder U2

The results (true effect is 0.1):

adjusts for U1,U2,U3,U4   β̂ = 0.144   95% CI (0.083, 0.206)   BIC = 2907.6
drops U2                  β̂ = 0.171   95% CI (0.127, 0.214)   BIC = 2902.1

Three things happened:

1. The model missing a true confounder has the better BIC. By the usual logic, it is the model we would choose.
2. Its confidence interval excludes the true value. We would report a significant effect of ~0.17 with confidence and be completely wrong.
3. The fully adjusted model, penalized for carrying a “useless” predictor, is the one whose interval actually covers the truth.

The two models fit the data equally well, but they estimate different parameters. Only one of them is the causal quantity we wanted.

Bayesian Model Averaging Has the Same Blind Spot

Suppose we don’t trust any single selected model and decide to average them. Bayesian model averaging (BMA) weights each model by its posterior probability and blends the estimates.

The catch is that with the usual flat prior over models, posterior weights are driven by marginal likelihood — the same criterion that just betrayed our estimates. Running through all 2⁵ = 32 outcome models and applying BIC weights reveals where the probability settles.

Figure 2: BMA puts ~94% of its weight on the model that omits confounder U2. Linking model selection to the exposure model (BAC’s idea, with ω→∞) flips nearly all the weight onto fully adjusted models.

BMA assigns roughly 90 percent of its mass to the version that drops U2. Averaging does not help when the weights themselves point at the wrong target. Wang et al. observed exactly this pattern and found that BMA’s 95% intervals covered the truth less often than the nominal level would suggest [1].

The BAC Idea: Ask the Exposure Model Who the Confounders Are

A confounder, by definition, must predict the exposure as well as the outcome. We should therefore not select covariates using only the outcome model. Instead, we fit a second model for the exposure and let that model inform the outcome model’s selections.

BAC formalizes this with two linked variable-selection problems:

• an outcome model: Y as a function of X and candidate covariates,
• an exposure model: X as a function of the same candidates.

These two are joined by a dependence parameter ω, representing the prior odds that a covariate enters the outcome model given that it is already in the exposure model. At ω = 1, the models decouple and we recover plain BMA. As ω → ∞, any covariate the exposure model selects is forced into the outcome model, giving full confounding adjustment [1].

If this logic resembles propensity scores, the intuition is sound. The exposure model is a propensity model, and standard causal inference practice includes covariates related to treatment assignment when constructing an adjustment [2]. BAC encodes that intuition as a prior. The crucial difference is that BAC does not force a single hand-picked specification — it captures the uncertainty in that selection and propagates it into the final inference.

A minimal, enumerated version of BAC is shown below. Genuine BAC uses MCMC, but with only five covariates we can brute-force every model combination.

from itertools import product

def bic_weights(bics):
    b = np.array(bics)
    w = np.exp(-0.5 * (b - b.min()))
    return w / w.sum()

# All 32 outcome models
out_models, out_bics, betas = [], [], []
for inc in product([0,1], repeat=M):
    m = fit_outcome(inc)
    out_models.append(inc); out_bics.append(m.bic)
    betas.append(m.params["X"])
w_bma = bic_weights(out_bics)

# All 32 exposure models -> per-covariate inclusion probability
exp_bics, exp_models = [], []
for inc in product([0,1], repeat=M):
    feats = [c for c, i in zip(cols, inc) if i]
    Xm = sm.add_constant(df[feats]) if feats else np.ones((n,1))
    exp_bics.append(sm.OLS(df["X"], Xm).fit().bic)
    exp_models.append(inc)
w_exp = bic_weights(exp_bics)
incl_prob = [sum(w for w, inc in zip(w_exp, exp_models) if inc[j])
             for j in range(M)]

# BAC-style prior with omega -> infinity:
# outcome models must contain every covariate the exposure model demands
required = [p > 0.5 for p in incl_prob]
ok = np.array([all(inc[j] for j in range(M) if required[j])
               for inc in out_models])
w_bac = w_bma * ok
w_bac = w_bac / w_bac.sum()

print("exposure inclusion probs:", np.round(incl_prob, 3))
print("BMA estimate:", np.dot(w_bma, betas))
print("BAC estimate:", np.dot(w_bac, betas))

The printed output shows that the exposure model correctly recovers the confounder structure:

exposure inclusion probs: [1.    1.    1.    0.031 0.032]
BMA estimate: 0.169
BAC estimate: 0.144

The inclusion probabilities for U1, U2, and U3 are all approximately one, while U4 and U5 have probabilities very close to zero. When the outcome-model weights are conditioned on this knowledge, the treatment effect estimate moves from 0.169 down to 0.144. Figure 2 shows how the model weights flip decisively toward the fully adjusted specifications.

To confirm this is not a lucky draw from a single dataset, repeating the experiment 300 times validates that consulting the exposure model reliably steers variable selection toward the covariates that matter for unbiased causal estimation — a result that standard predictive criteria consistently fail to achieve.