Survival Analysis with causaldef

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.