The question sensitivity analysis answers
A matched observational study (Module 6) assumes no unobserved confounding. Rosenbaum asks instead: how strong would hidden bias have to be to explain away the result?
Within a matched pair, two people look identical on observed covariates. Without hidden bias, each was equally likely to be the treated one. Rosenbaum's model lets an unobserved factor change those odds by up to a factor Γ:
1/Γ ≤ (odds that person A is treated) / (odds that person B is treated) ≤ Γ.
Γ = 1 means no hidden bias. Γ = 2 means one member of a pair could be twice as likely to be treated because of something we did not measure. For each Γ we compute the worst-case p-value. The Γ at which it crosses 0.05 measures the study's robustness.
!pip install -q numpy pandas scipy statsmodels scikit-learn matplotlib causaldata
import numpy as np, pandas as pd
from scipy.stats import norm, rankdata
import statsmodels.formula.api as smf
from sklearn.neighbors import NearestNeighbors
import matplotlib.pyplot as plt
from causaldata import nsw_mixtape, cps_mixtape
Matched pairs from the LaLonde data
We rebuild Module 6's propensity score and form one-to-one matched pairs, without replacement, so each pair has one trainee and one distinct CPS respondent.
nsw = nsw_mixtape.load_pandas().data
obs = pd.concat([nsw[nsw.treat == 1], cps_mixtape.load_pandas().data], ignore_index=True)
obs["u74"] = (obs.re74 == 0).astype(int); obs["u75"] = (obs.re75 == 0).astype(int)
ps = smf.logit("treat ~ age + I(age**2) + I(age**3) + educ + I(educ**2) + marr + nodegree + black + hisp"
" + re74 + re75 + u74 + u75 + educ:re74", data=obs).fit(disp=0).predict(obs)
obs["lps"] = np.log(ps / (1 - ps))
treated = obs[obs.treat == 1].sort_values("lps", ascending=False) # hardest-to-match first
controls = obs[obs.treat == 0].copy()
pairs = []
for _, t in treated.iterrows():
nn = NearestNeighbors(n_neighbors=1).fit(controls[["lps"]])
j = controls.index[nn.kneighbors([[t.lps]])[1][0][0]]
pairs.append((t.re78, controls.loc[j, "re78"]))
controls = controls.drop(j)
pairs = pd.DataFrame(pairs, columns=["y_treated", "y_control"])
pairs["diff"] = pairs.y_treated - pairs.y_control
print(f"{len(pairs)} pairs | mean difference ${pairs['diff'].mean():,.0f} | median ${pairs['diff'].median():,.0f}")
Rosenbaum bounds step by step
The test statistic is Wilcoxon's signed-rank T⁺: rank the absolute pair differences, then add up the ranks of the positive ones. Under hidden bias Γ, the chance that a pair's difference is positive (when the programme has no effect) is at most Γ / (1 + Γ). Using that worst case in the normal approximation gives an upper bound on the p-value.
def rosenbaum_pvalue(diffs, gamma):
d = diffs[diffs != 0]
ranks = rankdata(np.abs(d))
t_plus = ranks[d > 0].sum()
n = len(d)
p_plus = gamma / (1 + gamma)
mean = p_plus * n * (n + 1) / 2
var = p_plus * (1 - p_plus) * n * (n + 1) * (2 * n + 1) / 6
return 1 - norm.cdf((t_plus - mean) / np.sqrt(var)) # one-sided, worst case
gammas = np.round(np.arange(1.0, 2.01, 0.05), 2)
pvals = [rosenbaum_pvalue(pairs["diff"].to_numpy(), g) for g in gammas]
table = pd.DataFrame({"Gamma": gammas, "worst-case p": pvals}).round(4)
print(table.iloc[::2].to_string(index=False))
crit = next((g for g, p in zip(gammas, pvals) if p > 0.05), None)
print(f"Significance at 5% is lost at Gamma ≈ {crit}")
plt.plot(gammas, pvals, marker="o"); plt.axhline(0.05, color="red", ls="--")
plt.xlabel("Γ (hidden bias)"); plt.ylabel("Worst-case one-sided p-value"); plt.title("Rosenbaum sensitivity analysis"); plt.show()
At Γ = 1 this is the usual signed-rank test. As Γ grows, the worst-case p-value rises. For the LaLonde pairs significance disappears at Γ ≈ 1.1: a hidden factor that changed the odds of entering training by only 10% within pairs could explain the result away. The matched estimate was close to the experimental benchmark, but on its own it would be a fragile observational finding, which is exactly the lesson of the Dehejia–Wahba versus Smith–Todd debate.
Reading Γ and reporting it
Report it like this: "The result would remain significant at 5% unless an unobserved factor made otherwise similar people up to Γ times more likely to enter training." Then give the reader a benchmark. Γ around 1.1–1.3 means fragile: a weak hidden factor could explain it away. Classic robust findings, such as heavy smoking and lung cancer, survive up to Γ ≈ 5–6.
To make Γ concrete, compare it with an observed covariate: how much does omitting, say, 1975 earnings change the treatment odds within pairs? If the effect survives a hidden factor as strong as your strongest observed confounder, that is reassuring.
Simulate a matched study with a large true effect (pair differences ~ Normal(3, 2), n = 300) and a small one (Normal(0.4, 2)). Find the critical Γ for each. What determines robustness to hidden bias: statistical significance, effect size, or both?
Show a solution
rng = np.random.default_rng(3)
for mu in (3.0, 0.4):
d = rng.normal(mu, 2, 300)
g_crit = next((g for g in np.arange(1, 50, 0.05) if rosenbaum_pvalue(d, g) > 0.05), np.inf)
print(f"mean difference {mu}: p(Γ=1) = {rosenbaum_pvalue(d, 1):.2g}, critical Γ ≈ {g_crit:.2f}")
# Both: a large effect relative to noise (not just a tiny p-value from a big sample) gives robustness.
Checklist and recap
- I can explain Γ as a bound on hidden bias in treatment odds within matched pairs.
- I can compute Rosenbaum bounds with the signed-rank statistic.
- I report the critical Γ with a benchmark the reader can understand.