## ----include = FALSE----------------------------------------------------------
eval_surv <- getRversion() >= "4.0.0" && requireNamespace("survival", quietly = TRUE)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(causaldef)
library(ggplot2) # Ensure ggplot2 is loaded

# DAG Helper
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 = "#CD5C5C", shape = 21, stroke = 1.5) + # IndianRed
    ggplot2::geom_text(ggplot2::aes(label = name), fontface = "bold", size = 3.5, 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(hct_outcomes)

head(hct_outcomes)

## ----spec, eval = eval_surv---------------------------------------------------
# spec_surv <- causal_spec_survival(
#   data = hct_outcomes,
#   treatment = "conditioning_intensity",
#   time = "time_to_event",
#   event = "event_status",
#   # Note: The event setup needs binary 0/1 or Surv object logic
#   # We will convert it internally or preprocess
#   covariates = c("age", "disease_status", "kps"),
#   estimand = "ATE" # Difference in survival probability
# )

## ----prep, eval = eval_surv---------------------------------------------------
# # Create indicator for Primary Event (Relapse)
# hct_outcomes$relapse <- ifelse(hct_outcomes$event_status == "Relapse", 1, 0)
# 
# # Create indicator for Competing Event (Death)
# # Patients who die without relapse are removed from the risk set for relapse
# # but are NOT merely censored in the same way as loss-to-follow-up.
# hct_outcomes$death_cr <- ifelse(hct_outcomes$event_status == "Death", 1, 0)
# 
# # Check the distribution of events
# # Check the distribution of events
# table(hct_outcomes$event_status)
# 
# # Visualize Competing Risks
# coords <- data.frame(
#   name = c("Covariates", "CondIntensity", "Relapse", "Death"),
#   x = c(0, -1.5, 1.5, 1.5),
#   y = c(1, 0, 0.5, -0.5)
# )
# edges <- data.frame(
#   from = c("Covariates", "Covariates", "Covariates", "CondIntensity", "CondIntensity"),
#   to =   c("CondIntensity", "Relapse", "Death", "Relapse", "Death")
# )
# plot_dag(coords, edges, title = "Competing Risks Structure")

## ----real_spec, eval = eval_surv----------------------------------------------
# spec_hct <- causal_spec_survival(
#   data = hct_outcomes,
#   treatment = "conditioning_intensity",
#   time = "time_to_event",
#   event = "relapse",           # Primary event of interest
#   competing_event = "death_cr", # Competing event
#   covariates = c("age", "disease_status", "kps", "donor_type"),
#   estimand = "ATE"
# )
# 
# print(spec_hct)

## ----real_est, eval = eval_surv-----------------------------------------------
# # Estimate deficiency
# # Note: In survival settings, this measures how well we can recover the
# # interventional survival curves.
# results_hct <- estimate_deficiency(
#   spec_hct,
#   methods = c("unadjusted", "iptw"),
#   n_boot = 0 # Set > 0 for bootstrap confidence intervals
# )
# 
# print(results_hct)

## ----surv_model, message=FALSE, warning=FALSE, eval = eval_surv---------------
# # 1. Visualize Survival Curves
# # We can use the standard survival package to visualize the unadjusted curves first
# library(survival)
# hct_outcomes$rfs_event <- ifelse(hct_outcomes$event_status != "Censored", 1, 0)
# surv_obj <- Surv(time = hct_outcomes$time_to_event, event = hct_outcomes$rfs_event)
# km_fit_naive <- survfit(surv_obj ~ conditioning_intensity, data = hct_outcomes)
# 
# plot(km_fit_naive, col = c("blue", "red"), lwd = 2,
#      main = "Relapse-Free Survival (Unadjusted)",
#      xlab = "Time (Months)", ylab = "Survival Probability")
# legend("topright", legend = levels(hct_outcomes$conditioning_intensity),
#        col = c("blue", "red"), lwd = 2)
# 
# # 2. Estimate Causal Effect (RMST Difference)
# # We use estimate_effect() on our deficiency object to compute the RMST difference
# # using the IPTW weights we found were effective (low deficiency).
# 
# # Horizon = 24 months
# horizon <- 24
# # We reuse the deficiency object (and its learned weights) for the new estimand.
# # Since deficiency measures the distributional match, the same weights apply
# # to any outcome derived from the survival curve.
# results_hct$spec$estimand <- "RMST"
# results_hct$spec$horizon <- 24
# 
# effect_iptw <- estimate_effect(
#   results_hct,
#   target_method = "iptw",
#   contrast = c("Myeloablative", "Reduced") # Defined as Treated vs Control
# )
# 
# print(effect_iptw)
# 
# # Compare with Unadjusted Effect
# effect_naive <- estimate_effect(
#   results_hct,
#   target_method = "unadjusted",
#   contrast = c("Myeloablative", "Reduced")
# )
# 
# print(effect_naive)
# 
# cat(sprintf("\nCausal RMST Difference (IPTW): %.2f months\n", effect_iptw$estimate))
# cat(sprintf("Biased RMST Difference (Unadjusted): %.2f months\n", effect_naive$estimate))

## ----safety_floor, eval = eval_surv-------------------------------------------
# # Calculate policy regret bounds (transfer penalty + minimax floor)
# # Utility range is [0, horizon] (months of survival).
# # The bound is tied to the previously estimated IPTW deficiency proxy.
# bound <- policy_regret_bound(results_hct, utility_range = c(0, horizon), method = "iptw")
# print(bound)