#!/usr/bin/env python3
"""
Script 06: Sensitivity Analysis with E-Value (Python)
Atlas: Causality Atlas — Ch 12 (Sensitivity Analysis)
"""
import numpy as np
from scipy import stats

def e_value(rr):
    """Compute E-value from observed risk ratio."""
    if rr < 1:
        rr = 1 / rr
    return rr + np.sqrt(rr * (rr - 1))

def e_value_ci(rr, lo, hi):
    """E-value for point estimate and CI bound."""
    e_point = e_value(rr)
    if lo > 1:
        e_ci = e_value(lo)
    elif hi < 1:
        e_ci = e_value(hi)
    else:
        e_ci = 1.0
    return e_point, e_ci

def odds_to_rr(or_val, p0):
    """Convert odds ratio to risk ratio given baseline risk p0."""
    return or_val / ((1 - p0) + p0 * or_val)

def bias_factor_binary(P_U1_A1, P_U1_A0, RR_UY):
    """Compute bias factor for binary confounder with binary outcome."""
    num = (P_U1_A1 / P_U1_A0) * (RR_UY - 1) / RR_UY + 1
    den = (P_U1_A0 / P_U1_A1) * (RR_UY - 1) / RR_UY + 1
    return num / den

print("=== Sensitivity Analysis: E-Value ===\n")

# ---- 1. Binary outcome (RR) ----
print("--- Binary Outcome ---")
e_point, e_ci = e_value_ci(2.5, 1.8, 3.5)
print(f"Observed RR = 2.5 [95% CI: 1.8, 3.5]")
print(f"E-value (point) = {e_point:.2f}")
print(f"E-value (CI lower) = {e_ci:.2f}")
print(f"Interpretation: Unmeasured confounder would need RR ≥ {max(e_point, e_ci):.2f} with both exposure and outcome.")

# ---- 2. Continuous outcome ----
print("\n--- Continuous Outcome ---")
# Cohen's d to approximate RR conversion
d = 0.45
se = 0.12
e_point, _ = e_value_ci(np.exp(d), np.exp(d - 1.96*se), np.exp(d + 1.96*se))
print(f"Observed SMD = {d:.2f} (SE = {se})")
print(f"E-value = {e_point:.2f}")

# ---- 3. Multiple confounder scenarios ----
print("\n--- Bias Factor Scenario Analysis ---")
scenarios = [
    (0.4, 0.2, 2.0, "Weak"),
    (0.5, 0.2, 3.0, "Moderate"),
    (0.6, 0.15, 4.0, "Strong"),
]
rr_obs = 3.0
for pu1_a1, pu1_a0, rr_uy, label in scenarios:
    bf = bias_factor_binary(pu1_a1, pu1_a0, rr_uy)
    print(f"  {label}: P(U|A=1)={pu1_a1}, P(U|A=0)={pu1_a0}, RR_UY={rr_uy}")
    print(f"         Bias factor = {bf:.3f}, Corrected RR = {rr_obs/bf:.3f}")

# ---- 4. Statistical significance sensitivity ----
print("\n--- Gamma-like p-value testing ---")
np.random.seed(42)
diff = np.random.normal(0.5, 1, 100)
for gamma in [1.0, 1.25, 1.5, 1.75, 2.0]:
    p = 2 * (1 - stats.norm.cdf(np.mean(diff) / (np.std(diff) / np.sqrt(len(diff))) / gamma))
    print(f"  Γ = {gamma:.2f}: adjusted p = {p:.4f}")

print("\n=== Script complete ===")