A question that sounds easy but isn’t

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

However, 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. It is the uncertainty about which variables to adjust for, and they proposed a fix called Bayesian Adjustment for Confounding (BAC) [1]. It also has a famous frequentist twin that most people don’t realise is the same insight: post-double-selection and double machine learning [3][4].

In this post, I will demonstrate three things. I will show you the failure mode in runnable Python, build the fix, and then finally show you 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 needs?

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.

Now, 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. Also, 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 “effect” is now actually 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.

Now consider the entire argument reduced to two t-statistics. On the simulated data below, the variable U2 shows a t-statistic of 1.2 inside the outcome model. That yields p equal to 0.23. It appears indistinguishable from noise. Inside the exposure model, the same variable shows a t-statistic of 31.5. Same variable, same data set. One model calls it irrelevant. The other calls it central.

Every routine we normally invoke examines only the first column. This outcome follows from arithmetic. 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. That 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. Two quantities. Fit rules register only one.BIC remains selection consistent. With sufficient observations, it would retain U2. The present difficulty is finite sample. The asymptotic disappearance of bias offers little solace when the data set at hand contains roughly 1,000 rows. The Code below follows the pattern already described.

The failure mode, in Python

Let’s simulate the setup in Figure 1 and fit two outcome models, i.e., the fully adjusted one and one that drops U2.

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
Here, three things happened:

  • The model missing a true confounder has the better BIC.By the usual logic, it is the model we are going to choose.
  • Its confidence interval excludes the true value.We would report a significant effect of ~0.17 with confidence and actually be completely wrong.
  • 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. However, 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.

Here, the catch is that with the usual flat prior over models, the posterior weights are driven by marginal likelihood. Again, it is the same criterion that just betrayed our estimates. Run through all 2⁵ = 32 outcome models, apply the usual BIC weights, and watch where the probability settles.

BMA assigns roughly 90 percent of its mass to the version that drops U2. Hence, we can conclude that averaging does not help when the weights themselves point at the wrong target. Wang et al. observed exactly the same pattern in their paper 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. This means we should not select covariates using only the outcome model; instead, we should fit a second model for the exposure and let that model veto 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 ω, which represents 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 that the exposure model wants is forced into the outcome model, which gives us the full confounding adjustment [1].

If this logic reminds you of propensity scores, that is absolutely intuitive. The exposure model is a propensity model, and the textbook advice from the causal inference literature includes covariates related to treatment assignment when constructing our adjustment [2]. This is exactly the intuition BAC encodes as a prior. The crucial difference is that BAC does not force us to commit to a single hand‑picked specification. Instead, 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 possible 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.

A single dataset could be a lucky draw, so I repeated the whole experiment 300 times.

This bias is not random noise. It is structural. The prediction‑driven weights systematically favour under‑adjusted models, and the final estimates inherit that systematic slant.

The catch: an instrument looks exactly like a confounder

This is the part that most write‑ups of BAC skip, and it explains why the ω knob actually exists. The rule “adjust for anything that predicts the exposure” is not a very safe advice. Let’s consider a variable Z that predicts X strongly but has no path to Y whatsoever except through X. That variable is an instrument, and an exposure model cannot tell it apart from a genuine confounder.

Now let’s see what happens when we obey the exposure model in a setting where UC is an unmeasured confounder (as there always is in real data) and Z is a pure instrument:

Uc = rng.normal(size=n) # unmeasured confounder Z = rng.normal(size=n) # instrument: affects X only X = Uc + Z + rng.normal(size=n) Y = 0.1*X + Uc + rng.normal(size=n) # Z does NOT appear here
Over 500 replications with a true effect of 0.1, the results are striking:

no adjustment β̂ = 0.434 (bias +0.334) adjust for Z β̂ = 0.601 (bias +0.501)
Adjusting for Z made things 50% worse. The estimate did not become less efficient but it became more biased. This phenomenon is known as bias amplification or Z‑bias [5][6]. The intuition works as follows: conditioning on Z strips out the part of X’s variation that was clean (the instrument‑driven part). This leaves behind a residual X that is more dominated by the unmeasured confounder than before. So, we have basically thrown away our good variation and kept the contaminated variation, so the residual confounding got concentrated.

At , ignoring leaves a bias of , whereas adjusting for Z leaves . And we cannot detect this distinction from an exposure model alone. In a simulation where Z is a pure instrument and C is a genuine confounder, both entering X identically, we obtain the following:

exposure model t-statistics: Z = 32.2 C = 33.6These numbers are statistically indistinguishable. However, causally they are opposite. No amount of data will separate them, because the difference lies in the DAG, not in the covariance matrix. This is precisely why BAC’s is a dial and not a switch. The setting is the aggressive version that forces in every exposure predictor. A finite softens the rule for exposure. The choice of is a causal assumption we are making, not a parameter we can estimate from the data. If you take one thing from this post, take this sentence.

The same idea, three ways:

BAC dates from 2012, but the core insight was independently rebuilt in econometrics. It eventually became Double Machine Learning, and is now used under various different names

  • Post‑double‑selection (Belloni, Chernozhukov & Hansen, 2014 [3]: The procedure is quite simple. We run a lasso of Y on the covariates and run a lasso of X on the covariates. Then take the union of the selected variables, and refit OLS on that union. The lasso‑on‑X step exists for exactly the reason BAC’s exposure model exists i.e., to rescue confounders that are weak on Y but strong on X. The impact is real. Using the theoretically motivated plug‑in penalty from the BCH paper on the same simulation, the outcome lasso drops- U2 in 100 out of 300 datasets (33% of the time). Single‑selection (lasso on Y only) yields a mean estimate of 0.112 with bias , while double‑selection (the union) yields a mean of 0.101 with bias and 95% confidence interval coverage of 0.94. A third of the time, the outcome lasso throws away a true confounder but the exposure lasso always catches it, and the interval coverage returns to nominal.
  • Double/debiased machine learning (Chernozhukov et al., 2018) [4]: This is the same two‑model structure generalised. We predict Y from the covariates and X from the covariates. Then we take the residuals of both, and regress residual‑on‑residual. Cross‑fitting is added to eliminate the overfitting bias. Now the two nuisance models can be gradient boosting or random forests instead of lassos, and the effect estimate remains consistent. A sanity check on the same data‑generating process confirms that DML lands at 0.101 i.e., the same place that BAC and double selection reach.
  • Outcome‑adaptive lasso (Shortreed & Ertefaie, 2017) [7]: Here is the counter‑punch, which exists precisely because of the bias‑amplification problem described above. Outcome‑adaptive lasso builds the propensity model but penalises each covariate by how weakly it relates to the outcome. It deliberately down-weights variables that predict exposure and nothing else. This is the exact opposite instinct from BAC with . It is a bet that pure instruments are the bigger danger; BAC’s bet is that dropped confounders are. Both bets are defensible.

So what is the conclusion?

  • We need to stop letting predictive criteria pick our adjustment set. BIC, AIC, lasso paths, and stepwise p-values all score models on outcome fit, and none of them knows what a confounder is. A model can win on all of these metrics and still hand us a biased effect.
  • Second, we need to fit the exposure model even if we never run BAC. It costs us four lines of code, and it is the only diagnostic in our toolkit that can see the right‑hand column of that t-statistic table. When a covariate is very prominent in the exposure model but silent in the outcome model, we have just found the variable that our selection procedure is about to throw away.
  • But, also we should not blindly obey the exposure model. Before we include a strong exposure‑predictor, ask the question the data cannot answer: could this variable affect Y through any path other than X? If the answer is yes, it is a confounder and we should include it. If the answer is genuinely no, it is an instrument, and adjusting for it will amplify whatever unmeasured confounding we have left. This is the step where we have to think, and no estimator will do it for us.
  • The asymmetry between under‑adjustment and over‑adjustment should shape our decisions. Over-adjustment for variables that are merely prognostic usually costs efficiency, whereas adjusting for instruments can actually increase bias. [1].

Finally, we need to treat adjustment uncertainty as real uncertainty. If our confidence interval is conditional on one hand‑selected model, it understates what we do not know.

The uncomfortable summary is that the tools most of us reach for first were built for a different job. Prediction asks, “Which model explains Y?” Causal inference asks, “Which model isolates the effect of X on Y?” Those questions have different answers. Knowing that is most of the battle, and knowing that the second question carries an assumption that no amount of data will resolve for us.

References

[1] C. Wang, G. Parmigiani and F. Dominici, Bayesian Effect Estimation Accounting for Adjustment Uncertainty (2012), Biometrics, 68(3), 661–671.

[2] A. Gelman, J. Hill and A. Vehtari, Regression and Other Stories (2020), Cambridge University Press, Ch. 20.

[3] A. Belloni, V. Chernozhukov and C. Hansen, Inference on Treatment Effects after Selection among High-Dimensional Controls (2014), Review of Economic Studies, 81(2), 608–650.

[4] V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey and J. Robins, Double/Debiased Machine Learning for Treatment and Structural Parameters (2018), The Econometrics Journal, 21(1), C1–C68.

[5] J. Pearl, Invited Commentary: Understanding Bias Amplification (2011), American Journal of Epidemiology, 174(11), 1223–1227.

[6] J. A. Myers, J. A. Rassen, J. J. Gagne, K. F. Huybrechts, S. Schneeweiss, K. J. Rothman, M. M. Joffe and R. J. Glynn, Effects of Adjusting for Instrumental Variables on Bias and Precision of Effect Estimates (2011), American Journal of Epidemiology, 174(11), 1213–1222.

[7] S. M. Shortreed and A. Ertefaie, Outcome-Adaptive Lasso: Variable Selection for Causal Inference (2017), Biometrics, 73(4), 1111–1122.