#!/usr/bin/env Rscript
# -----------------------------------------------------------------------------
# Script 06: Sensitivity Analysis with E-Value (R)
# Atlas: Causality Atlas — Ch 12 (Sensitivity Analysis)
# -----------------------------------------------------------------------------

library(EValue)

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

# ---- 1. Binary outcome (RR) ----
cat("--- Binary Outcome ---\n")
e_rr <- e_value(est = 2.5, lo = 1.8, hi = 3.5)
print(e_rr)

# ---- 2. Continuous outcome (SMD) ----
cat("\n--- Continuous Outcome (SMD) ---\n")
e_smd <- e_value(est = 0.45, se = 0.12)
cat(sprintf("Observed effect: d = 0.45 (SE = 0.12)\n"))
print(e_smd)

# ---- 3. Odds ratio (rare outcome) ----
cat("\n--- Odds Ratio (rare outcome) ---\n")
e_or <- e_value(est = 1.8, lo = 1.3, hi = 2.5, true = 1)
print(e_or)

# ---- 4. E-value for CI bound ----
cat("\n--- CI Bound Analysis ---\n")
e_ci <- e_value(est = 2.0, lo = 1.05, hi = 3.8)
cat("Point estimate E-value and CI bound E-value:\n")
print(e_ci)

# ---- 5. Rosenbaum-type sensitivity for matched study ----
cat("\n--- Rosenbaum Sensitivity (Matched Pairs) ---\n")
# Simulate matched pairs
set.seed(42)
n_pairs <- 100
Y_diff <- rnorm(n_pairs, 0.5, 1)

# Gamma sensitivity
for (gamma in c(1, 1.25, 1.5, 1.75, 2.0, 2.5)) {
  p_upper <- psignrank(sum(Y_diff > 0), n_pairs, gamma = gamma)
  cat(sprintf("Γ = %.2f: upper bound p-value = %.4f\n", gamma, 2 * min(p_upper, 1 - p_upper)))
}

# ---- 6. Bias factor for specific confounder scenario ----
cat("\n--- Bias Factor Calculation ---\n")
# Scenario: observed RR = 3.0, U has P(U=1|A=1)=0.5, P(U=1|A=0)=0.2, RR_UY = 2.5
bf <- bias_factor(0.5, 0.2, 2.5)
cat(sprintf("Observed RR = 3.0\n"))
cat(sprintf("Bias factor = %.3f\n", bf))
cat(sprintf("Corrected RR = %.3f\n", 3.0 / bf))
cat(sprintf("95%% CI (if observed CI = [2.0, 4.5]): [%.2f, %.2f]\n",
            2.0 / bf, 4.5 / bf))

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