# Survival Analysis with causaldef

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

    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)
    #>   id age disease_status kps  donor_type conditioning_intensity time_to_event
    #> 1  1  57          Early  89 HLA-Matched          Myeloablative         14.31
    #> 2  2  50          Early  85  Mismatched                Reduced         11.05
    #> 3  3  59   Intermediate  97  Mismatched          Myeloablative         20.23
    #> 4  4  65          Early  84   Unrelated          Myeloablative          2.90
    #> 5  5  54       Advanced  96   Unrelated                Reduced          2.22
    #> 6  6  48          Early  64  Mismatched          Myeloablative          1.25
    #>   event_status event_death
    #> 1      Relapse           0
    #> 2      Relapse           0
    #> 3      Relapse           0
    #> 4      Relapse           0
    #> 5      Relapse           0
    #> 6        Death           1
    table(hct$event_status)
    #>
    #> Censored    Death  Relapse
    #>      123      131      546

## 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.

    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

    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
    )
    #> ✔ Created survival causal specification: n=800, 131 events

    print(spec_hct)
    #>
    #> -- Survival Causal Specification ----------------------------------------
    #>
    #> * Treatment: conditioning_intensity ( binary )
    #> * Time: time_to_event (max = 86.08 )
    #> * Event: event_death ( 131 events)
    #> * Covariates: age, disease_status, kps, donor_type
    #> * Sample size: 800
    #> * Estimand: RMST
    #> * Horizon: 24

### Deficiency Estimation

    results_hct <- estimate_deficiency(
      spec_hct,
      methods = c("unadjusted", "iptw"),
      n_boot = 0
    )
    #> ℹ Inferred treatment value: Reduced
    #> ℹ Estimating deficiency: unadjusted
    #> ℹ Estimating deficiency: iptw

    print(results_hct)
    #>
    #> -- Deficiency Proxy Estimates (PS-TV) ------
    #>
    #>      Method  Delta SE            Quality
    #>  unadjusted 0.3030  - Insufficient (Red)
    #>        iptw 0.0644  -   Caution (Yellow)
    #> Note: delta is a propensity-score TV proxy (overlap/balance diagnostic).
    #>
    #> Best method: iptw (delta = 0.0644 )
    plot(results_hct, type = "bar")

The PS-TV proxy is about 0.303 before adjustment and about 0.064 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

    horizon <- 24
    bound <- policy_regret_bound(
      results_hct,
      utility_range = c(0, horizon),
      method = "iptw"
    )
    #> ℹ Transfer penalty: 1.5457 (delta = 0.0644)

    print(bound)
    #>
    #> -- Policy Regret Bounds -------------------------------------------------
    #>
    #> * Deficiency delta: 0.0644
    #> * Delta mode: point
    #> * Delta method: iptw
    #> * Delta selection: pre-specified method
    #> * Utility range: [0, 24]
    #> * Transfer penalty: 1.5457 (additive regret upper bound)
    #> * Minimax floor: 0.7728 (worst-case lower bound)
    #>
    #> Note: this is a plug-in bound using a deficiency proxy rather than an identified exact deficiency.
    #>
    #> Interpretation: Transfer penalty is 6.4 % of utility range given delta

With a 24-month utility range, the current IPTW proxy implies a transfer penalty of about 1.55 months and a minimax safety floor of about 0.77 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

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

    print(effect_iptw)
    #>
    #> -- Causal Effect Estimate ----------------------
    #> Method:    iptw
    #> Type:      RMST Difference
    #> Contrast:  Myeloablative vs Reduced
    #> Estimate:  0.9978
    #> Horizon:   24

    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

    frontier <- confounding_frontier(
      spec_hct,
      alpha_range = c(-2, 2),
      gamma_range = c(-2, 2),
      grid_size = 30
    )
    #> ℹ Computing benchmarks for observed covariates...
    #> ✔ Computed confounding frontier: 30x30 grid

    print(frontier)
    #>
    #> -- Confounding Frontier (Confounding Lower Bound) --------------------------
    #>
    #> * Grid size: 30 x 30
    #> * alpha range: [-2, 2]
    #> * gamma range: [-2, 2]
    #> * Model: gaussian
    #>
    #> -- Deficiency Summary --
    #>
    #> * Min delta: 5e-04
    #> * Max delta: 0.0525
    #> * Mean delta: 0.0156
    #>
    #> 37.3 % of grid has near-zero deficiency (delta < 0.01)
    plot(frontier)
    #> Warning: The following aesthetics were dropped during statistical transformation: fill.
    #> ℹ This can happen when ggplot fails to infer the correct grouping structure in
    #>   the data.
    #> ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
    #>   variable into a factor?

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:

    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.