## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----define-problem-----------------------------------------------------------
library(causaldef)
set.seed(123)

# Simulate a treatment decision problem
n <- 1000

# Covariates
age <- runif(n, 30, 70)
severity <- rbeta(n, 2, 5) * 10

# Confounded treatment assignment (sicker patients get treatment)
U <- rnorm(n)  # Unmeasured health status
ps_true <- plogis(-1 + 0.02 * age + 0.1 * severity + 0.5 * U)
A <- rbinom(n, 1, ps_true)

# Outcome: recovery score (0-100)
# True effect is heterogeneous
tau_true <- 10 + 0.2 * (age - 50)  # Older patients benefit more
Y <- 50 + tau_true * A - 0.3 * severity + 5 * U + rnorm(n, sd = 5)

# Clip to valid range
Y <- pmin(100, pmax(0, Y))

df <- data.frame(
  age = age,
  severity = severity,
  A = A,
  Y = Y
)

## ----estimate-deficiency------------------------------------------------------
spec <- causal_spec(
  data = df,
  treatment = "A",
  outcome = "Y",
  covariates = c("age", "severity")
)

# Estimate deficiency with multiple methods
def_results <- estimate_deficiency(
  spec,
  methods = c("unadjusted", "iptw", "aipw"),
  n_boot = 100
)

print(def_results)

## ----plot-deficiency, eval=FALSE----------------------------------------------
# plot(def_results, type = "bar")

## ----policy-bounds------------------------------------------------------------
# Define utility range (outcome is 0-100)
utility_range <- c(0, 100)

# Suppose our policy achieves 5% observed regret
obs_regret <- 5

# Compute bound
bounds <- policy_regret_bound(
  deficiency = def_results,
  utility_range = utility_range,
  obs_regret = obs_regret
)

print(bounds)

## ----plot-safety, eval=FALSE--------------------------------------------------
# # Show how safety floor varies with deficiency
# plot(bounds, type = "safety_curve")

## ----interpret----------------------------------------------------------------
cat("=== Policy Deployment Decision ===\n\n")

delta_best <- min(def_results$estimates)
M <- diff(utility_range)
transfer_penalty <- M * delta_best
minimax_floor <- 0.5 * M * delta_best

cat(sprintf("Best achievable deficiency: %.3f\n", delta_best))
cat(sprintf("Transfer penalty (M*delta): %.1f points\n", transfer_penalty))
cat(sprintf("Minimax safety floor (M/2*delta): %.1f points\n", minimax_floor))
cat(sprintf("Observed regret: %.1f points\n", obs_regret))

if (!is.null(bounds$regret_bound)) {
  cat(sprintf("Worst-case regret: %.1f points\n", bounds$regret_bound))
}

cat("\n")

# Decision thresholds
if (delta_best < 0.05) {
  cat("✓ EXCELLENT: Deficiency < 5%. High confidence in policy.\n")
} else if (delta_best < 0.10) {
  cat("⚠ MODERATE: Deficiency 5-10%. Proceed with monitoring.\n")
} else {
  cat("✗ CAUTION: Deficiency > 10%. Consider RCT before deployment.\n")
}

## ----confounding-frontier-----------------------------------------------------
# Map the confounding frontier
frontier <- confounding_frontier(
  spec,
  alpha_range = c(-2, 2),
  gamma_range = c(-2, 2),
  grid_size = 30
)

# Find the safe region
safe_region <- subset(frontier$grid, delta < 0.1)
cat(sprintf(
  "Safe operating region covers %.1f%% of confounding space\n",
  100 * nrow(safe_region) / nrow(frontier$grid)
))

## ----plot-frontier, eval=FALSE------------------------------------------------
# plot(frontier, type = "heatmap", threshold = c(0.05, 0.1, 0.2))

## ----grf-example, eval=FALSE--------------------------------------------------
# # Estimate deficiency using causal forests
# if (requireNamespace("grf", quietly = TRUE)) {
#   def_grf <- estimate_deficiency(
#     spec,
#     methods = c("aipw", "grf"),
#     n_boot = 50
#   )
# 
#   print(def_grf)
# 
#   # Get individual treatment effect predictions
#   kernel_grf <- def_grf$kernel$grf
#   if (!is.null(kernel_grf$tau_hat)) {
#     cat("\nHeterogeneous Effects Detected:\n")
#     cat(sprintf("ATE from forest: %.2f\n", kernel_grf$ate))
#     cat(sprintf("CATE range: [%.2f, %.2f]\n",
#                 min(kernel_grf$tau_hat),
#                 max(kernel_grf$tau_hat)))
#   }
# }

## ----monitoring, eval=FALSE---------------------------------------------------
# # Example: Re-estimate deficiency on new data
# new_data <- ...  # Your production data
# 
# new_spec <- causal_spec(
#   new_data,
#   treatment = "A",
#   outcome = "Y",
#   covariates = c("age", "severity")
# )
# 
# # Quick check
# def_monitor <- estimate_deficiency(
#   new_spec,
#   methods = "iptw",
#   n_boot = 50
# )
# 
# # Alert if deficiency increased
# if (def_monitor$estimates["iptw"] > 1.5 * delta_best) {
#   warning("Distribution shift detected! Deficiency increased.")
# }