---
title: "Survival Analysis with causaldef"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Survival Analysis with causaldef}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

**causaldef** supports survival analysis and a separate competing-risks workflow. This vignette focuses on a reproducible death-endpoint RMST analysis for the bundled HCT data, then shows how to switch to the dedicated competing-risks interface when relapse and death must be modeled separately.

The deficiency and regret-bound calculations below run on all supported runtimes. The effect-estimation chunk that relies on the survival runtime requires `R >= 4.0` in the current support matrix.

## Setup

```{r 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)
```

## Death-Endpoint Survival Workflow

We start with a simple survival specification where death is the binary event and the estimand is 24-month restricted mean survival time (RMST). This keeps the workflow aligned with the current survival API and with the manuscript example.

```{r 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")
```

### Specification

```{r 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 Estimation

```{r deficiency}
results_hct <- estimate_deficiency(
  spec_hct,
  methods = c("unadjusted", "iptw"),
  n_boot = 0
)

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

The PS-TV proxy is about `r formatC(results_hct$estimates["unadjusted"], digits = 3, format = "f")` before adjustment and about `r formatC(results_hct$estimates["iptw"], digits = 3, format = "f")` after IPTW. That is a meaningful reduction in the observational-to-interventional gap, but it is still large enough to remain decision-relevant on a 24-month utility scale.

In the HCT setting this matters because treatment assignment is confounded by indication: healthier patients are more likely to receive myeloablative conditioning. IPTW reduces that imbalance, but the remaining proxy shows that the observational comparison should not be treated as near-randomized.

### Regret Bounds on the Survival Scale

```{r regret}
horizon <- 24
bound <- policy_regret_bound(
  results_hct,
  utility_range = c(0, horizon),
  method = "iptw"
)

print(bound)
```

With a 24-month utility range, the current IPTW proxy implies a transfer penalty of about `r formatC(bound$transfer_penalty, digits = 2, format = "f")` months and a minimax safety floor of about `r formatC(bound$minimax_floor, digits = 2, format = "f")` months. This is a material amount of residual decision risk: even after observed confounding adjustment, the observational-to-interventional gap is not negligible relative to clinically meaningful survival differences.

### RMST Effect Estimation

```{r effect, eval = eval_surv_effect}
effect_iptw <- estimate_effect(
  results_hct,
  target_method = "iptw",
  contrast = c("Myeloablative", "Reduced")
)

print(effect_iptw)
```

```{r 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.")
```

When this chunk is available, the RMST estimate should be interpreted together with the transfer penalty. A large estimated survival benefit can still be fragile if it is of the same order as the regret bound; conversely, a modest residual proxy is more reassuring when the treatment effect is much larger than the transfer penalty.

## Sensitivity Analysis via the Confounding Frontier

```{r frontier}
frontier <- confounding_frontier(
  spec_hct,
  alpha_range = c(-2, 2),
  gamma_range = c(-2, 2),
  grid_size = 30
)

print(frontier)
plot(frontier)
```

The confounding frontier maps how the deficiency proxy changes as hypothetical unmeasured confounding becomes stronger on the treatment and outcome paths. Use it as a sensitivity map, not as proof of robustness: if the plausible region of confounding strength overlaps higher-deficiency territory, the study conclusions should be reported as sensitive to hidden bias.

## Dedicated Competing-Risks Workflow

When relapse and death must be treated as distinct event types, prefer the competing-risks interface instead of collapsing them into a single binary endpoint:

```{r 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)
```

This returns a competing-risks deficiency object for the selected event of interest and avoids treating competing events as ordinary censoring.

## Takeaway

For survival outcomes, `causal_spec_survival()` plus `estimate_deficiency()` gives a practical workflow for quantifying the observational-to-interventional gap on a clinically meaningful scale. For multi-state event processes such as relapse versus death, move to `causal_spec_competing()` and `estimate_deficiency_competing()` so the event structure is represented explicitly.