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

Deficiency Estimation

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

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

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.