Causal Inference III(a) — Confounding as a Parameter¶
Turning a sensitivity threshold into a posterior, calibrated against the confounders we can see¶
The sensitivity notebook ended on a number that is hard to act on. Nearest-neighbour propensity matching recovered the experimental benchmark almost exactly — USD 1,792 against a randomized truth of USD 1,794 — and Rosenbaum bounds then showed the result is destroyed by a confounder shifting within-pair treatment odds by a factor of only $\Gamma = 1.21$.
That is a tipping point. It says where the conclusion breaks and nothing about whether we are likely to be past it. A referee asking "is 1.21 a lot?" gets no answer from the number itself, because $\Gamma$ is a worst case over all confounders of that strength, and worst cases are rarely where reality sits.
This notebook asks the question the other way round. Instead of finding the confounder strength at which the result dies, put a prior on how strong the unobserved confounder plausibly is and integrate it out. The result is a posterior on the treatment effect that already accounts for hidden confounding — a distribution rather than a breakpoint.
The move that makes this honest rather than arbitrary is benchmarking: the prior on the unobserved confounder is calibrated against the confounders we did measure. "Plausible" then means no stronger than the eight covariates already in the model, which is a defensible claim rather than a guess. Imbens (2003) made exactly this argument on exactly this dataset.
Python/PyMC lead. First in the Bayesian selection-on-observables group.
1. The estimate, and the bias an omitted confounder would cause¶
Start from regression adjustment on the observational sample: $Y$ on treatment $W$ and the eight covariates. If a confounder $U$ is omitted, the estimated effect is biased by a quantity with a clean structure,
$$\text{bias} \;=\; \gamma_U \cdot \delta_U,$$
where $\gamma_U$ is how much $U$ moves the outcome (dollars per standard deviation of $U$, holding the covariates fixed) and $\delta_U$ is how imbalanced $U$ is between treated and control after adjusting for the covariates (in standard deviations of $U$). A confounder matters only if it does both: shift the outcome and differ across arms.
Both quantities are measurable for the covariates we do observe, which is what makes calibration possible. Before trusting the formula, the notebook checks it: for each observed covariate, drop it, refit, and compare the actual change in the estimate against $\gamma_j\delta_j$.
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, statsmodels.api as sm
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)
Xr=obs[cov].values.astype(float)
Xs=(Xr-Xr.mean(0))/Xr.std(0) # standardized, so gamma is per SD
def ols(y, cols):
A=sm.add_constant(np.column_stack(cols)); m=sm.OLS(y,A).fit()
return m
full=ols(Y,[W,Xs]); tau_hat=full.params[1]; tau_se=full.bse[1]
print(f"randomized benchmark (the answer key) ${bench:,.0f}")
print(f"regression-adjusted estimate on the obs sample ${tau_hat:,.0f} (SE {tau_se:,.0f})")
print(f" naive difference in means ${Y[W==1].mean()-Y[W==0].mean():,.0f}")
print(f" n = {len(Y)} ({int(W.sum())} treated, {int((1-W).sum())} control)")
g++ not available, if using conda: `conda install gxx`
randomized benchmark (the answer key) $1,794 regression-adjusted estimate on the obs sample $1,548 (SE 781) naive difference in means $-635 n = 614 (185 treated, 429 control)
# benchmark every observed covariate: outcome strength, imbalance, and the bias it alone would cause
rows=[]
for j,nm in enumerate(cov):
g = full.params[2+j] # dollars per SD of X_j, covariates held fixed
others=[Xs[:,k] for k in range(len(cov)) if k!=j]
d = ols(Xs[:,j], [W]+others).params[1] # SDs of X_j, treated vs control, adjusted
drop = ols(Y, [W]+others).params[1] - tau_hat # what actually happens if X_j is omitted
rows.append((nm, g, d, g*d, drop))
bm=pd.DataFrame(rows, columns=["covariate","gamma (per SD)","delta (SD)","predicted bias","actual bias if dropped"])
print(bm.to_string(index=False, formatters={"gamma (per SD)":"{:,.0f}".format,
"delta (SD)":"{:.2f}".format,
"predicted bias":"{:,.0f}".format,
"actual bias if dropped":"{:,.0f}".format}))
err=np.abs(bm["predicted bias"]-bm["actual bias if dropped"]).max()
print(f"\nlargest gap between the bias formula and the refit: ${err:,.0f}")
print("The formula is not an approximation being taken on trust -- dropping each covariate reproduces")
print("gamma x delta to the dollar, so the same arithmetic can be applied to a confounder we cannot see.")
print()
bm["|bias|"]=bm["predicted bias"].abs()
top=bm.sort_values("|bias|",ascending=False)
print("Which covariate is the worst confounder is not the one you would guess:")
for _,r in top.head(3).iterrows():
print(f" {r['covariate']:9} gamma ${r['gamma (per SD)']:>8,.0f}/SD delta {r['delta (SD)']:>6.2f} SD -> bias ${r['predicted bias']:>7,.0f}")
print()
print(f"{top.iloc[0]['covariate'].upper()} dominates, and NOT because it matters most for earnings. re74 has by far the largest")
print(f"outcome coefficient (${bm.loc[bm.covariate=='re74','gamma (per SD)'].iloc[0]:,.0f} per SD) yet contributes less bias, because after adjusting for the")
print(f"other covariates it is only {abs(bm.loc[bm.covariate=='re74','delta (SD)'].iloc[0]):.2f} SD out of balance. Black is barely half as predictive of earnings")
print(f"but sits {abs(bm.loc[bm.covariate=='black','delta (SD)'].iloc[0]):.2f} SD apart between the arms -- 84% of NSW trainees against 20% of CPS controls.")
print()
print("Confounding is a PRODUCT, so a variable can be strongly prognostic and harmless, or weakly")
print("prognostic and devastating. Ranking candidate confounders by how well they predict the outcome --")
print("the usual instinct -- gets this dataset wrong.")
covariate gamma (per SD) delta (SD) predicted bias actual bias if dropped
age 128 0.13 16 16
educ 1,061 0.19 198 198
black -607 1.12 -678 -678
hispan 161 0.18 30 30
married 200 -0.24 -48 -48
nodegree 125 0.17 22 22
re74 1,918 -0.17 -327 -327
re75 762 0.06 48 48
largest gap between the bias formula and the refit: $0
The formula is not an approximation being taken on trust -- dropping each covariate reproduces
gamma x delta to the dollar, so the same arithmetic can be applied to a confounder we cannot see.
Which covariate is the worst confounder is not the one you would guess:
black gamma $ -607/SD delta 1.12 SD -> bias $ -678
re74 gamma $ 1,918/SD delta -0.17 SD -> bias $ -327
educ gamma $ 1,061/SD delta 0.19 SD -> bias $ 198
BLACK dominates, and NOT because it matters most for earnings. re74 has by far the largest
outcome coefficient ($1,918 per SD) yet contributes less bias, because after adjusting for the
other covariates it is only 0.17 SD out of balance. Black is barely half as predictive of earnings
but sits 1.12 SD apart between the arms -- 84% of NSW trainees against 20% of CPS controls.
Confounding is a PRODUCT, so a variable can be strongly prognostic and harmless, or weakly
prognostic and devastating. Ranking candidate confounders by how well they predict the outcome --
the usual instinct -- gets this dataset wrong.
2. A prior for a confounder we cannot see¶
Now suppose an unobserved $U$ remains. We need priors on $(\gamma_U, \delta_U)$, and pulling them from thin air would make the whole exercise decorative. Instead they are calibrated against the observed covariates measured above: the analyst is asserting only that the confounder they missed is not stronger than the ones they caught.
Three calibrations are run, from generous to severe:
- as strong as a typical observed covariate — prior scales at the median of $|\gamma_j|$ and $|\delta_j|$;
- as strong as the strongest observed covariate — scales at the maximum, which for this dataset means a confounder rivalling pre-programme earnings;
- twice the strongest — a deliberately pessimistic case, asserting the analyst missed something worse than anything they measured.
Half-normal priors are used for the magnitudes, with the sign of the bias left free — a confounder is as free to hide a real effect as to manufacture a false one. Sampling uncertainty in $\hat\tau$ enters at the same time, so the posterior mixes both sources of doubt rather than reporting them separately.
g_abs=np.abs(bm["gamma (per SD)"].values); d_abs=np.abs(bm["delta (SD)"].values)
CAL=[("typical observed covariate", np.median(g_abs), np.median(d_abs)),
("strongest observed covariate", g_abs.max(), d_abs.max()),
("twice the strongest", 2*g_abs.max(), 2*d_abs.max())]
def sens_posterior(sg, sd, seed=0):
with pm.Model():
tau_obs = pm.Normal("tau_obs", tau_hat, tau_se) # sampling uncertainty in the estimate
gU = pm.HalfNormal("gU", sg) # outcome strength (a magnitude)
dU = pm.Normal("dU", 0.0, sd) # imbalance across arms, signed either way
bias= pm.Deterministic("bias", gU*dU)
pm.Deterministic("tau_true", tau_obs - bias)
with contextlib.redirect_stderr(io.StringIO()):
return pm.sample(4000, tune=1500, chains=4, cores=1, progressbar=False,
random_seed=seed, target_accept=0.9)
worst_obs=np.abs(bm["predicted bias"].values).max()
print(f" {'confounder assumed':30} {'implied |bias|':>22} {'posterior effect':>17} "
f"{'95% interval':>22} {'P(tau>0)':>9}")
print(f" {'':30} {'median / 95th pct':>22}")
res={}; diag=[]
for k,(nm,sg,sd) in enumerate(CAL):
it=sens_posterior(sg,sd,seed=k)
t=it.posterior["tau_true"].values.ravel()
b=np.abs(it.posterior["bias"].values.ravel())
rh=float(pd.to_numeric(az.summary(it,var_names=["tau_true"])["r_hat"],errors="coerce").max())
dv=int(it.sample_stats["diverging"].values.sum())
diag.append((nm,rh,dv))
res[nm]=t
impl=f"${np.median(b):,.0f} / ${np.percentile(b,95):,.0f}"
lo,hi=np.percentile(t,[2.5,97.5]); iv=f"[{lo:,.0f}, {hi:,.0f}]"
print(f" {nm:30} {impl:>22} {t.mean():>17,.0f} {iv:>22} {(t>0).mean():>9.3f}")
from scipy.stats import norm
iv0=f"[{tau_hat-1.96*tau_se:,.0f}, {tau_hat+1.96*tau_se:,.0f}]"
print(f"\n {'no confounding assumed':30} {'--':>22} {tau_hat:>17,.0f} {iv0:>22} "
f"{1-norm.cdf(0,tau_hat,tau_se):>9.3f}")
print(f" randomized benchmark ${bench:,.0f}")
print()
print(f"The 'implied |bias|' column is the calibration made auditable, against a worst OBSERVED bias of")
print(f"${worst_obs:,.0f}. Note that a half-normal with scale set to the largest observed value still places")
print("plenty of mass above it, so these labels are generous rather than strict -- and the third row is")
print("frankly beyond anything the data support, included to show how fast the analysis goes uninformative.")
print()
print(" convergence: " + "; ".join(f"{nm.split()[0]} r-hat {rh:.3f}, {dv} divergences" for nm,rh,dv in diag))
print()
print("Now the part that is easy to misread. The posterior MEAN barely moves across the three rows --")
print(f"{res[CAL[0][0]].mean():,.0f}, {res[CAL[1][0]].mean():,.0f}, {res[CAL[2][0]].mean():,.0f} -- while the interval goes from about $3,000 wide to $38,000.")
print("That is because the confounder is as free to hide a real effect as to manufacture a false one, so")
print("symmetric confounding uncertainty WIDENS the posterior without SHIFTING it. Sensitivity analysis")
print("is not a bias correction. It does not tell you the answer was really smaller; it tells you how much")
print("less you know than the unadjusted interval claimed.")
confounder assumed implied |bias| posterior effect 95% interval P(tau>0)
median / 95th pct
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.
typical observed covariate $27 / $158 1,559 [29, 3,112] 0.978
strongest observed covariate $765 / $4,628 1,555 [-3,218, 6,315] 0.827
twice the strongest $3,098 / $19,029 1,699 [-17,146, 20,909] 0.663 no confounding assumed -- 1,548 [17, 3,080] 0.976 randomized benchmark $1,794 The 'implied |bias|' column is the calibration made auditable, against a worst OBSERVED bias of $678. Note that a half-normal with scale set to the largest observed value still places plenty of mass above it, so these labels are generous rather than strict -- and the third row is frankly beyond anything the data support, included to show how fast the analysis goes uninformative. convergence: typical r-hat 1.000, 1 divergences; strongest r-hat 1.000, 0 divergences; twice r-hat 1.000, 0 divergences Now the part that is easy to misread. The posterior MEAN barely moves across the three rows -- 1,559, 1,555, 1,699 -- while the interval goes from about $3,000 wide to $38,000. That is because the confounder is as free to hide a real effect as to manufacture a false one, so symmetric confounding uncertainty WIDENS the posterior without SHIFTING it. Sensitivity analysis is not a bias correction. It does not tell you the answer was really smaller; it tells you how much less you know than the unadjusted interval claimed.
3. What this says that $\Gamma^\star$ could not¶
Rosenbaum's $\Gamma^\star = 1.21$ is a statement about a worst case: there exists a confounder of that strength which would overturn the result. It attaches no probability to that confounder existing, and it is silent on what happens for confounders that are weaker, stronger, or differently shaped.
The posterior answers the question a decision-maker actually has. The last section also lets us ask it in reverse — how strong would the hidden confounder have to be to explain the estimate away entirely — and put that on the scale of the covariates we measured, which is the only scale anyone has intuition for.
need = tau_hat # bias required to zero the estimate
strongest = (g_abs*d_abs).max(); which = cov[int(np.argmax(g_abs*d_abs))]
print(f"To wipe out the estimate entirely the hidden confounder must supply ${need:,.0f} of bias.")
print(f"The strongest OBSERVED confounder, {which}, supplies ${strongest:,.0f}.")
print(f"So U would have to be about {need/strongest:.1f}x as damaging as the worst covariate we measured.")
print()
for nm,t in res.items():
print(f" under a confounder {nm:30} P(effect is positive) = {(t>0).mean():.3f}, "
f"P(effect exceeds $1,000) = {(t>1000).mean():.3f}")
print()
print("This is the reporting a sensitivity threshold cannot produce. Gamma* = 1.21 says a confounder of")
print("that strength exists which would overturn the finding; it says nothing about how probable that is,")
print("and nothing at all about confounders of any OTHER strength. The posterior prices the whole range.")
To wipe out the estimate entirely the hidden confounder must supply $1,548 of bias. The strongest OBSERVED confounder, black, supplies $678. So U would have to be about 2.3x as damaging as the worst covariate we measured. under a confounder typical observed covariate P(effect is positive) = 0.978, P(effect exceeds $1,000) = 0.761 under a confounder strongest observed covariate P(effect is positive) = 0.827, P(effect exceeds $1,000) = 0.646 under a confounder twice the strongest P(effect is positive) = 0.663, P(effect exceeds $1,000) = 0.568 This is the reporting a sensitivity threshold cannot produce. Gamma* = 1.21 says a confounder of that strength exists which would overturn the finding; it says nothing about how probable that is, and nothing at all about confounders of any OTHER strength. The posterior prices the whole range.
fig,ax=plt.subplots(1,2,figsize=(13,4.6))
XLO,XHI=-9000,11000 # readable window; the widest calibration runs far beyond it
bins=np.linspace(XLO,XHI,90)
for (nm,t),c in zip(res.items(),[GREEN,ORANGE,RED]):
frac=np.mean((t>=XLO)&(t<=XHI))
ax[0].hist(np.clip(t,XLO,XHI),bins=bins,alpha=.5,color=c,density=True,
label=f"{nm} ({frac:.0%} shown)")
ax[0].set_xlim(XLO,XHI)
ax[0].axvline(0,color="k",lw=1)
ax[0].axvline(tau_hat,color=BLUE,ls="--",lw=1.6,label=f"no confounding ${tau_hat:,.0f}")
ax[0].axvline(bench,color=PURP,ls=":",lw=1.8,label=f"randomized truth ${bench:,.0f}")
ax[0].set_xlabel("treatment effect (USD), after integrating over hidden confounding")
ax[0].set_ylabel("posterior density")
ax[0].set_title("Confounding priced in: the mean holds, the spread does not")
ax[0].legend(fontsize=7)
gg=np.linspace(0,2.2*g_abs.max(),160); dd=np.linspace(0,2.2*d_abs.max(),160)
G,Dm=np.meshgrid(gg,dd); Z=tau_hat-G*Dm
cs=ax[1].contourf(G,Dm,Z,levels=14,cmap="RdYlGn"); plt.colorbar(cs,ax=ax[1],label="effect after bias (USD)")
ax[1].contour(G,Dm,Z,levels=[0],colors="k",linewidths=2.2)
ax[1].scatter(g_abs,d_abs,s=48,color="k",zorder=4)
OFF={"age":(7,-11),"educ":(7,4),"black":(9,2),"hispan":(7,5),
"married":(7,12),"nodegree":(-9,-11),"re74":(7,4),"re75":(7,-11)}
for j,nm in enumerate(cov):
dx,dy=OFF[nm]
ax[1].annotate(nm,(g_abs[j],d_abs[j]),fontsize=7.5,xytext=(dx,dy),
textcoords="offset points",ha="right" if dx<0 else "left")
ax[1].set_xlabel("gamma: outcome strength (USD per SD)"); ax[1].set_ylabel("delta: imbalance across arms (SD)")
ax[1].set_title("Black line = confounders that would zero the effect")
plt.tight_layout(); plt.show()
4. Summary¶
The bias formula is verified, not assumed. Dropping each observed covariate in turn reproduces $\gamma_j\delta_j$ to the dollar (largest discrepancy: $0). That is what licenses applying the same arithmetic to a confounder nobody can measure.
The worst observed confounder is not the one that predicts earnings best.
re74has by far the largest outcome coefficient — $1,918 per standard deviation — yet contributes only $327 of bias, because after adjustment it is barely out of balance.blackis about half as predictive but sits 1.1 standard deviations apart between arms (84% of NSW trainees against 20% of CPS respondents) and contributes $678, the largest of the eight. Confounding is a product: a variable can be strongly prognostic and harmless, or weakly prognostic and devastating. The usual instinct — rank candidate confounders by how well they predict the outcome — gets this dataset wrong, which is worth knowing before deciding which unmeasured variable to worry about.The result survives a confounder like the ones we measured, and not much more. Against a regression-adjusted estimate of $1,548, wiping the effect out entirely requires $1,548 of bias, or about 2.3 times the damage done by the worst covariate in the model. Under a confounder calibrated to a typical observed covariate, $\Pr(\text{effect} > 0) = 0.978$; calibrated to the strongest, it falls to 0.827.
Sensitivity analysis is not a bias correction, and the output makes that visible. The posterior mean barely moves across the three calibrations — 1,559, 1,555, 1,699 — while the 95% interval widens from roughly $3,000 to $38,000. Because a hidden confounder is as free to hide a real effect as to manufacture a false one, symmetric confounding uncertainty widens the posterior without shifting it. The analysis does not say the answer was really smaller; it says how much less is known than the unadjusted interval claimed.
What this adds over the threshold. Rosenbaum's $\Gamma^\star = 1.21$ states that a confounder of that strength exists which would overturn the finding. It attaches no probability to that confounder being real, and says nothing about confounders of any other strength. The posterior prices the whole range, on a scale calibrated to the covariates actually in the model — so the reportable claim becomes "the effect is positive with probability 0.83 even against a confounder as strong as the worst one we measured" rather than "it breaks at 1.21."
The limit. None of this is a test. If the true confounder is genuinely unlike anything observed — a fourth-generation selection process, say, or something the CPS never asked about — the calibration has nothing to say, and the third row of the table shows how quickly the analysis becomes uninformative once that is allowed. Benchmarking makes the prior defensible; it does not make it true.