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.