Cookbook: Survival Outcome with Censoring, End to End

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.

One complete, runnable script for a time-to-event outcome with right-censoring. Same four steps as every cookbook (see vignette("cookbook-continuous")); what changes is how a response is recorded — a survival response is either an exact event time or a censoring interval — and the estimand, here a log hazard ratio from the InferenceSurvival* family (Cox, Weibull AFT, log-rank, RMST, …).

Setup

EDI is not on CRAN yet, so install.packages("EDI") fails — install from R-universe (fallback: GitHub, subdir = "R/EDI"). Not evaluated here.

install.packages("EDI", repos = c("https://kapelner.r-universe.dev", "https://cloud.r-project.org"))
# or: remotes::install_github("kapelner/EDI", subdir = "R/EDI")

library(EDI)
set.seed(20260916)

n = 100
X = data.frame(
  age   = round(rnorm(n, 60, 10)),
  stage = sample(1:3, n, replace = TRUE)
)
true_log_hr = -0.5   # treatment lowers the hazard

Fixed design, recording events and censoring

Event times come from an exponential model; each subject is independently right-censored (lost to follow-up) with probability 0.3. EDI records a survival response as an interval (y_L, y_R]: an exact event at time t is y = t; right-censoring at t is y_L = t, y_R = Inf (“known event-free through t”). Left- and interval-censoring use the same representation (see Design$add_one_subject_response()); this cookbook uses the per-subject method so each case is explicit.

des = DesignFixedBernoulli$new(n = n, response_type = "survival", verbose = FALSE)
des$add_all_subjects_to_experiment(X)
des$assign_w_to_all_subjects()
w = des$get_w()

rate       = exp(-2 + true_log_hr * w + 0.02 * (X$age - 60) + 0.3 * (X$stage - 2))
event_time = rexp(n, rate)
censored   = rbinom(n, 1, 0.3) == 1
follow_up  = pmin(event_time, runif(n, 0, 2 * median(event_time)))  # observed time

for (i in seq_len(n)) {
  if (censored[i]) {
    des$add_one_subject_response(i, y_L = follow_up[i], y_R = Inf)   # right-censored at follow_up
  } else {
    des$add_one_subject_response(i, y = event_time[i])               # exact event
  }
}
table(censored = censored)
#> censored
#> FALSE  TRUE 
#>    68    32

Inference: Cox proportional hazards

inf = InferenceSurvivalCoxPHRegr$new(des, verbose = FALSE)
inf$num_cores = 1L
inf$compute_estimate()                         # log hazard ratio for treatment
#> [1] -1.109091
inf$compute_asymp_confidence_interval(alpha = 0.05)
#>       2.5%      97.5% 
#> -1.6429557 -0.5752272
inf$compute_asymp_two_sided_pval()
#> [1] 4.665464e-05

The randomization test replays the design’s assignment mechanism; the bootstrap resamples subjects carrying their (w, time, censoring) along. Note that a randomization confidence interval is deliberately not offered for the Cox-family (log-hazard-ratio) classes — the generic randomization CI inverts an accelerated-failure-time null on a log-time scale, which is not the same axis as a log hazard ratio (see NEWS.md, 1.0.1). The randomization p-value and the bootstrap CI are the right tools here.

inf$set_seed(1)
inf$compute_rand_two_sided_pval(r = 200, show_progress = FALSE)
#> [1] 0.01
inf$set_seed(1)
inf$compute_bootstrap_confidence_interval(alpha = 0.05, B = 200, show_progress = FALSE)
#>       2.5%      97.5% 
#> -1.6269088 -0.5115814

Everything at once

suite = InferenceSuite$new(des)
res = suite$run_all_inference(screen = TRUE, plots = FALSE, num_cores = 1L,
                              methods = c("wald", "score", "lik_ratio"), max_secs_per_class = 15)
#> inference       cov      estimand    est       se        pval       pval method     status 
#> class           mod                                                                        
#> ===========================================================================================
#> Classes 0/17  [                0%                 ] Status: Estimating...[KAvg Δ                    mean Δ      5.67      1.93      4.30e-03   wald            ok     
#> Classes 1/17  [=              5%                ] Estimated Time Left: 0s[KCox PH Regr     ~.       log hazar…  -1.11     0.272     4.67e-05   wald            ok     
#> Classes 2/17  [===            11%               ] Estimated Time Left: 0s[KCox PH Regr     ~.       log hazar…  -1.11     0.272     2.25e-05   score           ok     
#> Classes 3/17  [=====          17%               ] Estimated Time Left: 0s[KCox PH Regr     ~.       log hazar…  -1.11     0.272     3.57e-05   lik_ratio       ok     
#> Classes 4/17  [=======        23%               ] Estimated Time Left: 0s[KDep Cens Tran…  ~.       log time …  1.20      0.386     1.89e-03   wald            ok     
#> Classes 5/17  [=========      29%               ] Estimated Time Left: 0s[KDep Cens Tran…  ~.       log time …  1.20      0.386     6.09e-04   score           ok     
#> Classes 6/17  [===========    35%               ] Estimated Time Left: 0s[KDep Cens Tran…  ~.       log time …  1.20      0.386     2.29e-03   lik_ratio       ok     
#> Classes 7/17  [=============  41%               ] Estimated Time Left: 0s[KGehan Wilcox             gehan wil…  -0.258    0.0756    2.70e-04   wald            ok     
#> Classes 8/17  [============== 47%               ] Estimated Time Left: 0s[KKaplan-Meier Δ           survival …  7.34      3.62      4.26e-02   wald            ok     
#> Classes 9/17  [============== 52%               ] Estimated Time Left: 0s[KLog Rank                 log rank …  -0.539    0.153     4.09e-04   wald            ok     
#> Classes 10/17 [============== 58%               ] Estimated Time Left: 0s[KRestricted Av…           restr mea…  8.72      2.57      7.08e-04   wald            ok     
#> Classes 11/17 [============== 64% ==            ] Estimated Time Left: 0s[KStrat Cox PH …  ~.       log hazar…  -1.05     0.289     2.76e-04   wald            ok     
#> Classes 12/17 [============== 70% ====          ] Estimated Time Left: 0s[KStrat Cox PH …  ~.       log hazar…  -1.05     0.289     1.60e-04   score           ok     
#> Classes 13/17 [============== 76% ======        ] Estimated Time Left: 0s[KStrat Cox PH …  ~.       log hazar…  -1.05     0.289     1.95e-04   lik_ratio       ok     
#> Classes 14/17 [============== 82% ========      ] Estimated Time Left: 0s[KWeibull Regr    ~.       log time …  0.779     0.201     1.10e-04   wald            ok     
#> Classes 15/17 [============== 88% ==========    ] Estimated Time Left: 0s[KWeibull Regr    ~.       log time …  0.779     0.201     3.60e-06   score           ok     
#> Classes 16/17 [============== 94% ============  ] Estimated Time Left: 0s[KWeibull Regr    ~.       log time …  0.779     0.201     2.10e-04   lik_ratio       ok     
#> Classes 17/17 [============= 100% ==============] Estimated Time Left: 0s[K-------------------------------------------------------------------------------------------
#> Status: Completed in 1s.
#> 
#>   Estimand: gehan wilcoxon statistic (1 inferences)  : p =       NA
#>   Estimand: log hazard ratio (6 inferences)          : p = 0.000133
#>   Estimand: log rank martingale Δ (1 inferences)     : p =       NA
#>   Estimand: log time ratio (6 inferences)            : p = 0.000227
#>   Estimand: mean Δ (1 inferences)                    : p =       NA
#>   Estimand: restr mean survival time Δ (1 inferences): p =       NA
#>   Estimand: survival median Δ (1 inferences)         : p =       NA
#> 
#> Combined evidence against the sharp null across 7 estimands
#> (17 inferences, weighting = uniform within estimand):
#> p = 0.000355

Sequential design

Recording is identical per subject; the design decides each arrival’s treatment on the covariates before the outcome is known — as in a real trial, where events accrue after enrolment.

des_seq = DesignSeqOneByOneKK14$new(n = n, response_type = "survival", verbose = FALSE)
for (i in seq_len(n)) {
  w_i  = des_seq$add_one_subject_to_experiment_and_assign(X[i, , drop = FALSE])
  t_i  = rexp(1, exp(-2 + true_log_hr * w_i + 0.02 * (X$age[i] - 60) + 0.3 * (X$stage[i] - 2)))
  if (rbinom(1, 1, 0.3) == 1) {
    des_seq$add_one_subject_response(i, y_L = min(t_i, 3), y_R = Inf)
  } else {
    des_seq$add_one_subject_response(i, y = t_i)
  }
}

inf_seq = InferenceSurvivalCoxPHRegr$new(des_seq, verbose = FALSE)
inf_seq$num_cores = 1L
inf_seq$compute_estimate()
#> [1] -0.2158
inf_seq$set_seed(1)
inf_seq$compute_rand_two_sided_pval(r = 200, show_progress = FALSE)
#> [1] 0.49

Where to go next

Accelerated-failure-time alternative with a time-ratio estimand and a randomization CI: InferenceSurvivalWeibullRegr.
Nonparametric: InferenceSurvivalLogRank, InferenceSurvivalGehanWilcox, InferenceSurvivalRestrictedMeanDiff.
Interval-censored data and which classes accept it: Design$add_one_subject_response()’s documentation.
vignette("validation-evidence"): checked against survival::coxph.

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.