Cookbook: Survival Outcome with Censoring, End to End

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
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[@]dogado.de.

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.
Health stats visible at Monitor.