import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
from scipy.stats import norm, ttest_1samp
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"

def rosenbaum_p(diffs, Gamma):                      # worst-case p-value of Wilcoxon signed-rank at sensitivity Gamma
    d=diffs[diffs!=0]; q=np.argsort(np.argsort(np.abs(d)))+1.0; Wobs=q[d>0].sum()
    pp=Gamma/(1+Gamma); E=pp*q.sum(); V=pp*(1-pp)*(q**2).sum()
    return 1-norm.cdf((Wobs-E)/np.sqrt(V))
def gamma_star(diffs, grid=None):
    grid=np.linspace(1,15,1400) if grid is None else grid
    for G in grid:
        if rosenbaum_p(diffs,G)>0.05: return G
    return np.inf

# Two matched studies built to separate SIGNIFICANCE from ROBUSTNESS. Significance tracks the
# standardized effect times sqrt(n); robustness tracks the standardized effect alone. So a large
# effect in a small study and a small effect in a large study can reach the SAME p-value while
# tolerating very different amounts of hidden bias.
rng=np.random.default_rng(0)
big   = rng.normal(2.0, 2.0,  30)      # large effect (d ~ 1.0),  small study
small = rng.normal(0.5, 2.0, 480)      # small effect (d ~ 0.25), large study

rows=[]
for nm,x in [("large effect, small study", big), ("small effect, large study", small)]:
    rows.append((nm, len(x), x.mean(), x.std(ddof=1), x.mean()/x.std(ddof=1),
                 ttest_1samp(x,0).pvalue, rosenbaum_p(x,1.0), gamma_star(x)))
gs, gw = rows[0][7], rows[1][7]

print(f"  {'study':27} {'n':>4} {'mean':>7} {'sd':>6} {'d':>6} {'t-test p':>11} {'p at G=1':>10} {'Gamma*':>8}")
for nm,N,mu,sd,d_,tp,r1,G in rows:
    print(f"  {nm:27} {N:>4} {mu:>7.3f} {sd:>6.2f} {d_:>6.2f} {tp:>11.2e} {r1:>10.5f} {G:>8.1f}")

Gg=np.linspace(1,8,60)
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].hist(big,bins=15,alpha=.6,color=GREEN,density=True,label=f"large effect, n=30  (d={rows[0][4]:.2f})")
ax[0].hist(small,bins=40,alpha=.6,color=ORANGE,density=True,label=f"small effect, n=480  (d={rows[1][4]:.2f})")
ax[0].axvline(0,color="k",lw=.7); ax[0].set_xlabel("matched-pair outcome difference")
ax[0].set_title("Same p-value, very different effect-to-noise ratio"); ax[0].legend(fontsize=8)
ax[1].plot(Gg,[rosenbaum_p(big,G) for G in Gg],color=GREEN,lw=2.5,label=f"large effect, small study (Γ*={gs:.1f})")
ax[1].plot(Gg,[rosenbaum_p(small,G) for G in Gg],color=ORANGE,lw=2.5,label=f"small effect, large study (Γ*={gw:.1f})")
ax[1].axhline(0.05,color=RED,ls="--",label="p = 0.05"); ax[1].axvline(gs,color=GREEN,ls=":"); ax[1].axvline(gw,color=ORANGE,ls=":")
ax[1].set_xlabel("Γ (hidden-bias sensitivity parameter)"); ax[1].set_ylabel("worst-case p-value")
ax[1].set_title("Γ* = where the finding loses significance"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()

print()
print(f"Both studies are overwhelmingly significant, and the SMALL effect is the more significant of the two:")
print(f"p = {rows[1][5]:.1e} on 480 pairs against {rows[0][5]:.1e} on 30. A referee reading p-values alone would rank")
print(f"the 480-pair study as the stronger evidence.")
print()
print(f"Rosenbaum bounds reverse that. The large effect tolerates hidden bias up to Γ* = {gs:.1f}; the small one")
print(f"breaks at Γ* = {gw:.1f} -- a confounder shifting within-pair treatment odds by barely {100*(gw-1):.0f}% would erase it.")
print("The reason is that Γ measures how far the OBSERVED differences could be reshuffled before the sign")
print("pattern stops looking systematic, and that depends on the effect relative to the noise -- not on the")
print("sample size that turned it into a small p-value. Significance can always be bought with n; robustness")
print("to hidden bias cannot.")

  study                          n    mean     sd      d    t-test p   p at G=1   Gamma*
  large effect, small study     30   1.757   1.65   1.07    2.43e-06    0.00001      4.2
  small effect, large study    480   0.465   2.05   0.23    9.25e-07    0.00000      1.4

Both studies are overwhelmingly significant, and the SMALL effect is the more significant of the two:
p = 9.3e-07 on 480 pairs against 2.4e-06 on 30. A referee reading p-values alone would rank
the 480-pair study as the stronger evidence.

Rosenbaum bounds reverse that. The large effect tolerates hidden bias up to Γ* = 4.2; the small one
breaks at Γ* = 1.4 -- a confounder shifting within-pair treatment odds by barely 43% would erase it.
The reason is that Γ measures how far the OBSERVED differences could be reshuffled before the sign
pattern stops looking systematic, and that depends on the effect relative to the noise -- not on the
sample size that turned it into a small p-value. Significance can always be bought with n; robustness
to hidden bias cannot.

from sklearn.linear_model import LogisticRegression
obs=pd.read_csv("lalonde_obs.csv"); cov=["age","educ","black","hispan","married","nodegree","re74","re75"]
W=obs.treat.values; Y=obs.re78.values; X=obs[cov].values.astype(float)
Xs=(X-X.mean(0))/X.std(0); ps=LogisticRegression(penalty=None,max_iter=5000).fit(Xs,W).predict_proba(Xs)[:,1]
ti=np.where(W==1)[0]; ci=np.where(W==0)[0]; m=ci[np.abs(ps[ti][:,None]-ps[ci][None,:]).argmin(1)]
diffs=Y[ti]-Y[m]
gL=gamma_star(diffs, np.linspace(1,3,600))
print(f"LaLonde matched pairs: n={len(diffs)}, mean earnings difference = ${diffs.mean():,.0f} (matches experimental ~$1,794)")
print(f"worst-case p at Gamma = 1 (no hidden bias) = {rosenbaum_p(diffs,1.0):.4f}  "
      f"-- significant to begin with, so Gamma* means what it appears to")
print(f"Rosenbaum Gamma* = {gL:.2f}")
print(f"  -> a hidden confounder raising enrollment odds only {gL:.2f}x would overturn the result: FRAGILE")
Gg=np.linspace(1,2.5,60)
fig,ax=plt.subplots(figsize=(8,4))
ax.plot(Gg,[rosenbaum_p(diffs,G) for G in Gg],color=BLUE,lw=2.5)
ax.axhline(0.05,color=RED,ls="--",label="p = 0.05"); ax.axvline(gL,color=GREEN,ls=":",label=f"Γ* = {gL:.2f}")
ax.set_xlabel("Γ (hidden-bias sensitivity)"); ax.set_ylabel("worst-case p-value"); ax.set_title("LaLonde matching estimate: robust to hidden bias only up to Γ ≈ 1.2")
ax.legend(); plt.tight_layout(); plt.show()
print("The estimate recovered the experimental benchmark yet is fragile to hidden bias -- 'balance is testable, unconfoundedness")
print("is not' (previous notebook), made precise. A referee who asks 'what about unobserved confounders?' now has a number.")

LaLonde matched pairs: n=185, mean earnings difference = $1,792 (matches experimental ~$1,794)
worst-case p at Gamma = 1 (no hidden bias) = 0.0033  -- significant to begin with, so Gamma* means what it appears to
Rosenbaum Gamma* = 1.21
  -> a hidden confounder raising enrollment odds only 1.21x would overturn the result: FRAGILE

The estimate recovered the experimental benchmark yet is fragile to hidden bias -- 'balance is testable, unconfoundedness
is not' (previous notebook), made precise. A referee who asks 'what about unobserved confounders?' now has a number.

def evalue(rr):
    rr=rr if rr>=1 else 1.0/rr
    return rr+np.sqrt(rr*(rr-1))
examples=[("strong (smoking-scale)",3.9,3.0),("weak",1.3,1.05)]
for nm,rr,lo in examples:
    print(f"{nm:22s} RR={rr:.1f} (CI low {lo:.2f}):  E-value(point) = {evalue(rr):.2f},  E-value(CI) = {evalue(lo):.2f}")
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
rrs=np.linspace(1.05,6,100); ax[0].plot(rrs,[evalue(r) for r in rrs],color=PURP,lw=2.5)
for nm,rr,lo in examples: ax[0].scatter([rr],[evalue(rr)],s=70,zorder=5); ax[0].annotate(f"{nm}\nE={evalue(rr):.1f}",(rr,evalue(rr)),textcoords="offset points",xytext=(8,-4),fontsize=8)
ax[0].set_xlabel("observed risk ratio"); ax[0].set_ylabel("E-value"); ax[0].set_title("Bigger effects need stronger confounders to explain away")
# bias-factor frontier for the strong effect: RR_EU * RR_UD / (RR_EU+RR_UD-1) = RR  -> the curve of just-explaining confounders
RR=3.9; g=np.linspace(evalue(RR)*0.5,12,200); frontier=RR*(1-1.0/g)/((1)-RR/g)   # solve RR_UD given RR_EU=g
g2=np.linspace(evalue(RR),12,200); RR_UD=[(RR*(x-1))/(x-RR) if x>RR else np.nan for x in g2]
ax[1].plot(g2,RR_UD,color=RED,lw=2.5); ax[1].scatter([evalue(RR)],[evalue(RR)],s=80,color="k",zorder=5,label=f"E-value = {evalue(RR):.1f}")
ax[1].plot([1,12],[1,12],ls=":",color=GREY); ax[1].set_xlim(1,12); ax[1].set_ylim(1,12)
ax[1].set_xlabel("confounder–treatment association (RR)"); ax[1].set_ylabel("confounder–outcome association (RR)")
ax[1].set_title(f"Frontier of confounders that could explain RR={RR}"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"Smoking-scale RR {RR}: E-value {evalue(RR):.1f} -- a confounder would need ~{evalue(RR):.0f}-fold associations with BOTH")
print(f"smoking and cancer to explain it; none is known, so the finding is robust. Weak RR 1.3: E-value {evalue(1.3):.2f} -- a")
print("modest confounder could explain it. The E-value turns 'could confounding explain this?' into a number anyone can judge.")

strong (smoking-scale) RR=3.9 (CI low 3.00):  E-value(point) = 7.26,  E-value(CI) = 5.45
weak                   RR=1.3 (CI low 1.05):  E-value(point) = 1.92,  E-value(CI) = 1.28

Smoking-scale RR 3.9: E-value 7.3 -- a confounder would need ~7-fold associations with BOTH
smoking and cancer to explain it; none is known, so the finding is robust. Weak RR 1.3: E-value 1.92 -- a
modest confounder could explain it. The E-value turns 'could confounding explain this?' into a number anyone can judge.

Causal Inference II(b) — Sensitivity Analysis for Unobserved Confounding¶

Rosenbaum bounds and the E-value: how strong must a hidden confounder be to overturn the result?¶

1. Rosenbaum bounds — how much hidden bias would it take?¶

2. A humbling real example — the LaLonde matching estimate¶

3. The E-value — a scale-free confounding threshold¶

4. Summary¶