Causal Inference III(b) — Why the Propensity Score Resists Bayes¶
The assignment model is ancillary, a coherent posterior ignores it, and that is not obviously safe¶
Everything in the selection-on-observables group is organised around the propensity score. It is the balancing score that collapses eight covariates to one, the thing matching matches on, the thing weighting weights by, and the object whose estimation occupies most of the applied literature.
A likelihood-based Bayesian analysis throws it away.
The reason is a factorisation, not an opinion. Write the joint model for the data with parameters $\alpha$ governing treatment assignment and $\beta$ governing the outcome:
$$p(Y, W \mid X, \alpha, \beta) \;=\; \underbrace{p(W \mid X, \alpha)}_{\text{assignment}} \cdot \underbrace{p(Y \mid W, X, \beta)}_{\text{outcome}}.$$
If $\alpha$ and $\beta$ carry independent priors, the posterior for $\beta$ — and therefore for the treatment effect — does not involve $\alpha$ at all. The assignment model is ancillary. A Bayesian who writes down the full likelihood and turns the handle will get the same answer whether or not the propensity model is in the room.
This notebook does three things: demonstrates the ancillarity on the real data, shows the case where ignoring the design is genuinely dangerous, and explains why that danger is what pushed Bayesian causal inference toward response surfaces rather than propensity models.
Python/PyMC lead. Second in the Bayesian selection-on-observables group.
1. The demonstration on LaLonde¶
Two models on the same observational sample. The first is an outcome model alone: earnings on treatment and the eight covariates. The second is that plus a complete logistic model of treatment assignment — the propensity model the entire frequentist group is built around — fitted jointly in the same PyMC model.
If the factorisation above is right, the treatment-effect posterior should be unmoved. Not similar: unmoved, to within Monte Carlo error.
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings, logging, contextlib, io
warnings.filterwarnings("ignore")
import pymc as pm, arviz as az
logging.getLogger("pymc").setLevel(logging.ERROR)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
cov=["age","educ","black","hispan","married","nodegree","re74","re75"]
obs=pd.read_csv("lalonde_obs.csv"); exp=pd.read_csv("lalonde_exp.csv")
bench=exp.re78[exp.treat==1].mean()-exp.re78[exp.treat==0].mean()
W=obs.treat.values.astype(float); Y=obs.re78.values.astype(float)
X=obs[cov].values.astype(float); Xs=(X-X.mean(0))/X.std(0)
def fit(with_assignment, seed=1):
with pm.Model():
a = pm.Normal("a", 0, 10000)
tau = pm.Normal("tau", 0, 10000)
b = pm.Normal("b", 0, 5000, shape=Xs.shape[1])
s = pm.HalfNormal("s", 10000)
pm.Normal("y", a + tau*W + pm.math.dot(Xs, b), s, observed=Y)
if with_assignment: # the propensity model, fitted jointly
a0 = pm.Normal("a0", 0, 2)
al = pm.Normal("al", 0, 2, shape=Xs.shape[1])
pm.Bernoulli("w_obs", p=pm.math.invlogit(a0 + pm.math.dot(Xs, al)), observed=W)
with contextlib.redirect_stderr(io.StringIO()):
return pm.sample(3000, tune=1500, chains=4, cores=1, progressbar=False,
random_seed=seed, target_accept=0.9)
id_out, id_joint = fit(False), fit(True)
t0=id_out.posterior["tau"].values.ravel(); t1=id_joint.posterior["tau"].values.ravel()
mcse0=float(az.mcse(id_out, var_names=["tau"])["tau"].values)
mcse1=float(az.mcse(id_joint, var_names=["tau"])["tau"].values)
bench=exp.re78[exp.treat==1].mean()-exp.re78[exp.treat==0].mean()
print(f" {'model':38} {'effect':>9} {'sd':>8} {'95% interval':>22}")
for nm,t in (("outcome model only", t0), ("+ full propensity model, fitted jointly", t1)):
lo,hi=np.percentile(t,[2.5,97.5])
print(f" {nm:38} {t.mean():>9,.0f} {t.std():>8,.0f} {f'[{lo:,.0f}, {hi:,.0f}]':>22}")
print(f"\n difference in posterior mean: ${abs(t0.mean()-t1.mean()):,.0f}")
print(f" Monte Carlo error of each mean alone: ${mcse0:,.0f} and ${mcse1:,.0f}")
print(f" difference in posterior sd: ${abs(t0.std()-t1.std()):,.0f}")
print(f" randomized benchmark from the experimental sample: ${bench:,.0f}")
print()
print("The propensity model is not approximately irrelevant here -- it is exactly irrelevant, and the")
print("gap between the two answers is smaller than the noise in either one. Adding a complete model of")
print("who got treated, the object the whole frequentist group is organised around, changes nothing.")
g++ not available, if using conda: `conda install gxx`
model effect sd 95% interval outcome model only 1,542 775 [25, 3,070] + full propensity model, fitted jointly 1,559 783 [42, 3,083] difference in posterior mean: $18 Monte Carlo error of each mean alone: $8 and $7 difference in posterior sd: $8 randomized benchmark from the experimental sample: $1,794 The propensity model is not approximately irrelevant here -- it is exactly irrelevant, and the gap between the two answers is smaller than the noise in either one. Adding a complete model of who got treated, the object the whole frequentist group is organised around, changes nothing.
2. Why that should worry you¶
It is tempting to read section 1 as a convenience: one fewer model to fit. It is better read as a warning, because the propensity score was never decoration — it buys robustness to getting the outcome model wrong.
Inverse-probability weighting is consistent whenever the propensity model is right, no matter how badly the outcome is modelled. A pure outcome-model analysis has no such protection: if the regression is misspecified, the estimate inherits the misspecification. A coherent Bayesian, having established that the assignment model cannot enter, has thrown away exactly the information that would have provided the insurance.
Robins & Ritov (1997) turned this into a formal problem: in high dimensions there are settings where any procedure ignoring the known assignment mechanism performs badly, while one using it does fine. The construction below is a small, concrete version — the assignment mechanism is known by design, so there is no excuse for the estimator that cannot use it.
rng=np.random.default_rng(0)
N, P, TRUE_TAU, REPS = 600, 10, 2.0, 300
def one(seed):
r=np.random.default_rng(seed)
Z=r.normal(size=(N,P))
e=1/(1+np.exp(-(1.2*Z[:,0]-1.2*Z[:,1]+0.6*Z[:,2]))) # KNOWN assignment mechanism
Wd=r.binomial(1,e).astype(float)
g=4*np.sin(1.5*Z[:,0])+2*Z[:,1]**2-2*Z[:,2] # prognostic, and NOT linear
Yd=TRUE_TAU*Wd+g+r.normal(size=N)
# (a) outcome-model analysis, linear in Z -- misspecified. Posterior mean under vague priors = OLS.
A=np.column_stack([np.ones(N),Wd,Z])
tau_out=np.linalg.lstsq(A,Yd,rcond=None)[0][1]
# (b) weighting by the TRUE propensity score -- uses the design, ignores the outcome shape
w=np.where(Wd==1,1/e,1/(1-e))
tau_ipw=(np.sum(w*Wd*Yd)/np.sum(w*Wd))-(np.sum(w*(1-Wd)*Yd)/np.sum(w*(1-Wd)))
return tau_out,tau_ipw
res=np.array([one(s) for s in range(REPS)])
print(f" {REPS} replicates, n={N}, {P} covariates, true effect {TRUE_TAU:.1f}")
print(f" {'estimator':44} {'mean':>7} {'bias':>7} {'sd':>7} {'RMSE':>7}")
for k,nm in ((0,"outcome model, linear (misspecified)"),(1,"weighting by the TRUE propensity score")):
v=res[:,k]; bias=v.mean()-TRUE_TAU
print(f" {nm:44} {v.mean():>7.2f} {bias:>+7.2f} {v.std():>7.2f} "
f"{np.sqrt(bias**2+v.var()):>7.2f}")
print()
b_out=res[:,0].mean()-TRUE_TAU; b_ipw=res[:,1].mean()-TRUE_TAU
print("The bias is real and it goes the way the theory says. The outcome model is wrong in a way that is")
print("CORRELATED WITH TREATMENT -- assignment depends on Z0 and Z1 and the outcome bends in Z0 and Z1 --")
print(f"so its error does not average out: bias {b_out:+.2f}. Weighting by the known propensity score never")
print(f"models the outcome at all, and is essentially unbiased: {b_ipw:+.2f}.")
print()
print("But read the whole row before declaring a winner, because the honest result is a trade and not a")
print(f"victory. Weighting pays for its unbiasedness with {res[:,1].std()/res[:,0].std():.1f}x the spread, and on RMSE the")
print(f"MISSPECIFIED REGRESSION WINS -- {np.sqrt(b_out**2+res[:,0].var()):.2f} against {np.sqrt(b_ipw**2+res[:,1].var()):.2f}. An analyst minimising expected squared error")
print("would take the biased estimator here, and would be right to.")
print()
print("So the design information is worth something, and what it is worth is unbiasedness rather than")
print("accuracy. That is the precise form of the loss: the coherent Bayesian of section 1 cannot reach")
print("for this fix, because the factorisation says the assignment model carries no information about")
print("the effect. The information is real; the likelihood principle says it is not in the likelihood.")
300 replicates, n=600, 10 covariates, true effect 2.0 estimator mean bias sd RMSE outcome model, linear (misspecified) 2.23 +0.23 0.32 0.39 weighting by the TRUE propensity score 1.97 -0.03 1.02 1.02 The bias is real and it goes the way the theory says. The outcome model is wrong in a way that is CORRELATED WITH TREATMENT -- assignment depends on Z0 and Z1 and the outcome bends in Z0 and Z1 -- so its error does not average out: bias +0.23. Weighting by the known propensity score never models the outcome at all, and is essentially unbiased: -0.03. But read the whole row before declaring a winner, because the honest result is a trade and not a victory. Weighting pays for its unbiasedness with 3.2x the spread, and on RMSE the MISSPECIFIED REGRESSION WINS -- 0.39 against 1.02. An analyst minimising expected squared error would take the biased estimator here, and would be right to. So the design information is worth something, and what it is worth is unbiasedness rather than accuracy. That is the precise form of the loss: the coherent Bayesian of section 1 cannot reach for this fix, because the factorisation says the assignment model carries no information about the effect. The information is real; the likelihood principle says it is not in the likelihood.
3. What the Bayesian answer actually is¶
The resolution is not to bolt the propensity score onto a Bayesian analysis. Two-step recipes that draw $\alpha$ from its posterior and then match or weight are incoherent as a single model — they are pragmatic devices, not posteriors, and they are usually justified by frequentist properties rather than Bayesian ones.
The resolution is to notice what actually failed in section 2. It was not the absence of a propensity model. It was a misspecified outcome model: the linear regression could not represent a function that bends, and the bending was correlated with treatment.
That suggests an obvious repair — make the outcome model flexible enough that misspecification stops being the binding constraint — and it is the route the Bayesian causal literature took toward response surfaces rather than propensity scores.
The obvious repair, applied naively, does not work. It is worth seeing fail before seeing it done properly.
# does flexibility fix what the propensity score was insuring against?
from sklearn.ensemble import RandomForestRegressor
def one_flex(seed):
r=np.random.default_rng(seed)
Z=r.normal(size=(N,P))
e=1/(1+np.exp(-(1.2*Z[:,0]-1.2*Z[:,1]+0.6*Z[:,2])))
Wd=r.binomial(1,e).astype(float)
g=4*np.sin(1.5*Z[:,0])+2*Z[:,1]**2-2*Z[:,2]
Yd=TRUE_TAU*Wd+g+r.normal(size=N)
A=np.column_stack([np.ones(N),Wd,Z]); tau_lin=np.linalg.lstsq(A,Yd,rcond=None)[0][1]
# a FLEXIBLE outcome model, still with no propensity score anywhere
f=RandomForestRegressor(n_estimators=300,min_samples_leaf=5,random_state=seed,n_jobs=-1)
f.fit(np.column_stack([Wd,Z]),Yd)
y1=f.predict(np.column_stack([np.ones(N),Z])); y0=f.predict(np.column_stack([np.zeros(N),Z]))
return tau_lin,(y1-y0).mean()
rf=np.array([one_flex(s) for s in range(120)])
print(f" {'estimator':44} {'mean':>7} {'bias':>7} {'sd':>7}")
for k,nm in ((0,"outcome model, linear (misspecified)"),(1,"outcome model, flexible -- no propensity score")):
v=rf[:,k]; print(f" {nm:44} {v.mean():>7.2f} {v.mean()-TRUE_TAU:>+7.2f} {v.std():>7.2f}")
print()
bl=rf[:,0].mean()-TRUE_TAU; bf=rf[:,1].mean()-TRUE_TAU
print("That is not the result the section was set up to produce, and it is the more interesting one.")
print(f"Flexibility does NOT rescue the estimate. It makes it far worse: bias {bf:+.2f} against the linear")
print(f"model's {bl:+.2f}, roughly {abs(bf/bl):.0f} times larger and in the opposite direction. The flexible model")
print(f"attenuates the effect from {TRUE_TAU:.1f} toward zero, landing at {rf[:,1].mean():.2f}.")
print()
print("The mechanism is worth being precise about. The prognostic signal here is large -- g(Z) swings")
print("across roughly ten units -- while the treatment effect is 2. A regularized learner spends its")
print("splits where the variance is, which is Z, and the treatment indicator is left to soak up what")
print("little remains. The shrinkage that makes the model good at PREDICTING is what destroys it as an")
print("estimator of a causal contrast.")
print()
print("This is REGULARIZATION-INDUCED CONFOUNDING, and it is not an argument against flexible models --")
print("it is the argument for a particular kind of them. A model that separates the prognostic part from")
print("the treatment part, and lets the propensity score back in as a covariate to absorb the selection,")
print("fixes exactly this. That is Bayesian Causal Forests, and it is the third example in this group.")
estimator mean bias sd outcome model, linear (misspecified) 2.18 +0.18 0.31 outcome model, flexible -- no propensity score 0.48 -1.52 0.13 That is not the result the section was set up to produce, and it is the more interesting one. Flexibility does NOT rescue the estimate. It makes it far worse: bias -1.52 against the linear model's +0.18, roughly 8 times larger and in the opposite direction. The flexible model attenuates the effect from 2.0 toward zero, landing at 0.48. The mechanism is worth being precise about. The prognostic signal here is large -- g(Z) swings across roughly ten units -- while the treatment effect is 2. A regularized learner spends its splits where the variance is, which is Z, and the treatment indicator is left to soak up what little remains. The shrinkage that makes the model good at PREDICTING is what destroys it as an estimator of a causal contrast. This is REGULARIZATION-INDUCED CONFOUNDING, and it is not an argument against flexible models -- it is the argument for a particular kind of them. A model that separates the prognostic part from the treatment part, and lets the propensity score back in as a covariate to absorb the selection, fixes exactly this. That is Bayesian Causal Forests, and it is the third example in this group.
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].hist(t0,bins=60,alpha=.55,color=BLUE,density=True,label="outcome model only")
ax[0].hist(t1,bins=60,alpha=.55,color=ORANGE,density=True,label="+ propensity model, joint")
ax[0].axvline(bench,color=PURP,ls=":",lw=1.8,label=f"randomized truth ${bench:,.0f}")
ax[0].set_xlabel("treatment effect (USD)"); ax[0].set_ylabel("posterior density")
ax[0].set_title("Adding the propensity model changes nothing"); ax[0].legend(fontsize=8)
lab=["linear outcome\n(misspecified)","weighting by the\nTRUE propensity score","flexible outcome\n(no propensity score)"]
bp=ax[1].boxplot([res[:,0],res[:,1],rf[:,1]],labels=lab,showfliers=False,patch_artist=True)
for patch,c in zip(bp["boxes"],[RED,GREEN,PURP]): patch.set_facecolor(c); patch.set_alpha(.55)
ax[1].axhline(TRUE_TAU,color="k",ls="--",lw=1.5,label=f"true effect {TRUE_TAU:.1f}")
ax[1].set_ylabel("estimated effect")
ax[1].set_title("Three ways to be wrong about the same effect"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
4. Summary¶
The propensity score is exactly, not approximately, irrelevant to a coherent Bayesian analysis. Fitting a complete logistic model of treatment assignment jointly with the outcome model moves the effect posterior by USD 18 — against a Monte Carlo error of USD 8 in each mean, and a posterior standard deviation of about USD 780. The object the entire frequentist group is organised around cannot enter the effect posterior, because the likelihood factorises and the priors are independent.
What that costs is unbiasedness, and the size of the loss is a trade rather than a rout. In a constructed setting where the assignment mechanism is known and the outcome model is misspecified in a treatment-correlated way, the outcome-model estimate carries a bias of +0.23 against a true effect of 2.0, while weighting by the known propensity score is essentially unbiased at −0.03. But weighting pays with 3.2 times the spread, and on RMSE the misspecified regression wins — 0.39 against 1.02. An analyst minimising expected squared error would take the biased estimator, and would be right to. The design information buys unbiasedness, not accuracy.
The obvious repair fails, and fails informatively. Replacing the linear outcome model with a flexible one — no propensity score anywhere — does not recover the estimate. It attenuates the effect from 2.0 to 0.48, a bias of −1.52: roughly eight times the linear model's error and in the opposite direction. The prognostic signal swings across about ten units while the treatment effect is 2, so a regularized learner spends its splits on the covariates and leaves the treatment indicator to soak up the remainder. The shrinkage that makes a model good at prediction is what destroys it as an estimator of a causal contrast.
That last result has a name — regularization-induced confounding — and it is the reason this group needs a third example rather than stopping at "use a flexible outcome model". A model that separates the prognostic part from the treatment part, and admits the propensity score back in as a covariate to absorb the selection, repairs it. That is Bayesian Causal Forests.
The shape of the argument. A coherent Bayesian cannot use the assignment model, so robustness has to come from the outcome model instead; making the outcome model flexible introduces a new bias; fixing that bias requires letting the propensity score back in through the only door still open to it — as a predictor rather than as a weight. The propensity score returns, but not in the role the frequentist group gave it.