## ----include = FALSE----------------------------------------------------------
eval_surv <- getRversion() >= "4.0.0" && requireNamespace("survival", quietly = TRUE)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6, 
  fig.height = 4,
  warning = FALSE, 
  message = FALSE
)

## ----flash_forward, echo=FALSE, results='hide', fig.width=6, fig.height=4-----
library(causaldef)
library(ggplot2)

set.seed(2026)
n <- 500
W <- rnorm(n)
prob_A <- plogis(0.5 * W)
A <- rbinom(n, 1, prob_A)
Y <- 1 + 2 * A + 1 * W + rnorm(n)
Y_nc <- 0.5 * W + rnorm(n)
df <- data.frame(W = W, A = A, Y = Y, Y_nc = Y_nc)

spec <- causal_spec(
  data = df,
  treatment = "A",
  outcome = "Y",
  covariates = "W",
  negative_control = "Y_nc"
)

results <- estimate_deficiency(
  spec,
  methods = c("unadjusted", "iptw", "aipw"),
  n_boot = 50
)

bounds <- policy_regret_bound(results, utility_range = c(0, 10), method = "aipw")
plot(bounds, type = "safety_curve")

## ----setup--------------------------------------------------------------------
library(causaldef)
library(ggplot2)

set.seed(2026)

# DAG Helper (using ggplot2 only)
plot_dag <- function(coords, edges, title = NULL) {
  edges_df <- merge(edges, coords, by.x = "from", by.y = "name")
  colnames(edges_df)[c(3,4)] <- c("x_start", "y_start")
  edges_df <- merge(edges_df, coords, by.x = "to", by.y = "name")
  colnames(edges_df)[c(5,6)] <- c("x_end", "y_end")
  
  ggplot2::ggplot(coords, ggplot2::aes(x = x, y = y)) +
    ggplot2::geom_segment(data = edges_df, ggplot2::aes(x = x_start, y = y_start, xend = x_end, yend = y_end),
                 arrow = ggplot2::arrow(length = ggplot2::unit(0.3, "cm"), type = "closed"), 
                 color = "gray40", size = 1, alpha = 0.8) +
    ggplot2::geom_point(size = 14, color = "white", fill = "#4682B4", shape = 21, stroke = 1.5) +
    ggplot2::geom_text(ggplot2::aes(label = name), fontface = "bold", size = 4, color = "white") +
    ggplot2::ggtitle(title) + 
    ggplot2::theme_void(base_size = 14) +
    ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5, face = "bold", margin = ggplot2::margin(b = 10))) +
    ggplot2::coord_fixed()
}

## ----data_gen-----------------------------------------------------------------
# Simulate data
n <- 500
W <- rnorm(n) # Severity (Confounder)

# Treatment assignment biased by severity
prob_A <- plogis(0.5 * W) 
A <- rbinom(n, 1, prob_A)

# Outcome: Treatment effect = 2, Severity effect = 1
Y <- 1 + 2 * A + 1 * W + rnorm(n)

# Negative Control: Affected by W but not A
Y_nc <- 0.5 * W + rnorm(n)

df <- data.frame(W = W, A = A, Y = Y, Y_nc = Y_nc)

# Visualize the Causal Structure
coords <- data.frame(
  name = c("W", "A", "Y", "Y_nc"),
  x = c(0, -1, 1, 0.5),
  y = c(1, 0, 0, 1) # Pushing Y_nc up/out to side
)
edges <- data.frame(
  from = c("W", "W", "W", "A"),
  to =   c("A", "Y", "Y_nc", "Y")
)
plot_dag(coords, edges, title = "Scenario: Confounding + Negative Control")

## ----spec---------------------------------------------------------------------
spec <- causal_spec(
  data = df,
  treatment = "A",
  outcome = "Y",
  covariates = "W",
  negative_control = "Y_nc"
)
print(spec)

## ----deficiency---------------------------------------------------------------
# Compare Baseline (Unadjusted) vs IPTW vs AIPW
results <- estimate_deficiency(
  spec,
  methods = c("unadjusted", "iptw", "aipw"),
  n_boot = 50 
)

print(results)
plot(results, type = "bar")

## ----diagnostics--------------------------------------------------------------
diag_res <- nc_diagnostic(spec, method = "iptw")
print(diag_res)

## ----estimation---------------------------------------------------------------
# Estimate ATE
effect <- estimate_effect(results, target_method = "iptw")
print(effect)

## ----regret-------------------------------------------------------------------
# Utility range: [0, 10] outcome units
bounds <- policy_regret_bound(results, utility_range = c(0, 10), method = "aipw")
print(bounds)
plot(bounds, type = "safety_curve")

## ----frontier, fig.width=6, fig.height=5--------------------------------------
# What if there were an unobserved confounder U?
frontier <- confounding_frontier(spec, grid_size = 30)

# Calculate the delta corresponding to an acceptable transfer-penalty limit
# Suppose our max acceptable transfer penalty is 0.2 units, and M=10
# Transfer penalty = M * delta  =>  delta = transfer penalty / M
acceptable_transfer_penalty <- 0.2
M <- 10
delta_threshold <- acceptable_transfer_penalty / M

plot(frontier, type = "heatmap", threshold = c(delta_threshold)) 

## ----survival, eval = eval_surv-----------------------------------------------
# data("hct_outcomes")
# 
# # Event: Relapse (1) vs Censored (0)
# hct_outcomes$relapse <- ifelse(hct_outcomes$event_status == "Relapse", 1, 0)
# 
# spec_surv <- causal_spec_survival(
#   data = hct_outcomes,
#   treatment = "conditioning_intensity",
#   time = "time_to_event",
#   event = "relapse",
#   covariates = c("age", "disease_status", "kps"),
#   estimand = "RMST",
#   horizon = 24
# )
# 
# # Estimate deficiency (Distance between Obs and Target Survival Expts)
# surv_results <- estimate_deficiency(spec_surv, methods = c("unadjusted", "iptw"))
# 
# # Estimate Effect (RMST Difference + Survival Curves)
# surv_effect <- estimate_effect(surv_results, target_method = "iptw")
# 
# print(surv_effect)
# plot(surv_effect)