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

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

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
    ) +
    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)
hct <- transform(
  hct_outcomes,
  event_death = as.integer(as.character(event_status) == "Death")
)

head(hct)
table(hct$event_status)

## ----dag----------------------------------------------------------------------
coords <- data.frame(
  name = c("Covariates", "Conditioning", "Death"),
  x = c(0, -1.5, 1.5),
  y = c(1, 0, 0)
)
edges <- data.frame(
  from = c("Covariates", "Covariates", "Conditioning"),
  to = c("Conditioning", "Death", "Death")
)
plot_dag(coords, edges, title = "Death-Endpoint Survival Structure")

## ----spec---------------------------------------------------------------------
spec_hct <- causal_spec_survival(
  data = hct,
  treatment = "conditioning_intensity",
  time = "time_to_event",
  event = "event_death",
  covariates = c("age", "disease_status", "kps", "donor_type"),
  estimand = "RMST",
  horizon = 24
)

print(spec_hct)

## ----deficiency---------------------------------------------------------------
results_hct <- estimate_deficiency(
  spec_hct,
  methods = c("unadjusted", "iptw"),
  n_boot = 0
)

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

## ----regret-------------------------------------------------------------------
horizon <- 24
bound <- policy_regret_bound(
  results_hct,
  utility_range = c(0, horizon),
  method = "iptw"
)

print(bound)

## ----effect, eval = eval_surv_effect------------------------------------------
effect_iptw <- estimate_effect(
  results_hct,
  target_method = "iptw",
  contrast = c("Myeloablative", "Reduced")
)

print(effect_iptw)

## ----effect-note, results='asis', eval = !eval_surv_effect--------------------
# cat("Effect-estimation is skipped on runtimes without the required survival support. The deficiency and regret-bound calculations above still provide the main diagnostic quantities for this example.")

## ----frontier-----------------------------------------------------------------
frontier <- confounding_frontier(
  spec_hct,
  alpha_range = c(-2, 2),
  gamma_range = c(-2, 2),
  grid_size = 30
)

print(frontier)
plot(frontier)

## ----competing-risks, eval = FALSE--------------------------------------------
# spec_cr <- causal_spec_competing(
#   data = hct_outcomes,
#   treatment = "conditioning_intensity",
#   time = "time_to_event",
#   event = "event_status",
#   covariates = c("age", "disease_status", "kps", "donor_type"),
#   event_of_interest = "Relapse",
#   horizon = 24
# )
# 
# def_cr <- estimate_deficiency_competing(
#   spec_cr,
#   method = "cshr",
#   n_boot = 100
# )
# 
# print(def_cr)