Survival measures behavior in tasks with 100% vs 0% censoring

[bblodfon]

Based on #311, I did a bit of an investigation: what if a dataset/task has all observations censored or the complete opposite? Would the results we get from some of the measures implemented in this package make sense? See table of results and discussion points at the end:

library(mlr3verse)
#> Loading required package: mlr3
library(mlr3proba)
library(tibble)

set.seed(0)
# use the veteran dataset in two versions
vet_obs = survival::veteran
vet_obs$status = 1 # all dead

vet_cen = survival::veteran
vet_cen$status = 0 # all censored

# take of missing data
task_cen = as_task_surv(x = vet_cen, time = 'time', event = 'status')
task_obs = as_task_surv(x = vet_obs, time = 'time', event = 'status')
poe = po('encode')
task_cen = poe$train(list(task_cen))[[1L]]
task_obs = poe$train(list(task_obs))[[1L]]

# RSF by default computes both ranking scores and distr predictions
learner = lrn('surv.ranger', verbose = FALSE,
  id = 'SurvivalForestLogRank',
  label = 'Random Forest (Logrank splitrule)',
  fallback = lrn('surv.kaplan'),
  num.threads = 4,
  splitrule = 'logrank',
  num.trees = 50,
  mtry.ratio = 0.8,
  min.node.size = 3,
  importance = 'permutation'
)

# some measures of interest
rcll = msr('surv.rcll')
ibrier_proper = msr('surv.graf', proper = TRUE)
ibrier_improper = msr('surv.graf')
dcal = msr('surv.dcalib', B = 10, chisq = FALSE)
uno_c = msr('surv.cindex', weight_meth = 'G2')
harrell_c = msr('surv.cindex')

grid = benchmark_grid(tasks = list(task_cen, task_obs), learners = learner,
  resamplings = rsmp('cv', folds = 6))
bm = benchmark(design = grid, store_models = TRUE)
#> INFO  [01:48:33.871] [mlr3] Running benchmark with 12 resampling iterations
#> INFO  [01:48:33.934] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_cen' (iter 1/6)
#> INFO  [01:48:34.014] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_cen' (iter 2/6)
#> INFO  [01:48:34.075] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_cen' (iter 3/6)
#> INFO  [01:48:34.129] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_cen' (iter 4/6)
#> INFO  [01:48:34.180] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_cen' (iter 5/6)
#> INFO  [01:48:34.233] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_cen' (iter 6/6)
#> INFO  [01:48:34.302] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_obs' (iter 1/6)
#> INFO  [01:48:34.368] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_obs' (iter 2/6)
#> INFO  [01:48:34.442] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_obs' (iter 3/6)
#> INFO  [01:48:34.515] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_obs' (iter 4/6)
#> INFO  [01:48:34.643] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_obs' (iter 5/6)
#> INFO  [01:48:34.710] [mlr3] Applying learner 'SurvivalForestLogRank' on task 'vet_obs' (iter 6/6)
#> INFO  [01:48:34.784] [mlr3] Finished benchmark
dt = as.data.table(bm)

data_list = list()
index = 1
for(p in dt$prediction) {
  tsk_id = dt$task[[index]]$id
  task = dt$task[[index]]
 
  cen_num = sum(task$truth()[,2] == 0)
  train_set = dt$resampling[[index]]$train_set(dt$iteration[index])

  #print(p$score(rcll)) # due to issue we skip it
  ibrier_proper_score = p$score(ibrier_proper, task = task, train_set = train_set)
  ibrier_proper_test_score = p$score(ibrier_proper)
  ibrier_improper_score = p$score(ibrier_improper, task = task, train_set = train_set)
  ibrier_improper_test_score = p$score(ibrier_improper)

  dcal_score = p$score(dcal)
  unoc_score = p$score(uno_c, task = task, train_set = train_set)
  harc_score = p$score(harrell_c)

  data_list[[index]] = tibble(tsk_id = tsk_id, cen_num = cen_num,
    ibrier_proper_score = ibrier_proper_score,
    ibrier_proper_test_score = ibrier_proper_test_score,
    ibrier_improper_score = ibrier_improper_score,
    ibrier_improper_test_score = ibrier_improper_test_score,
    dcal_score = dcal_score, unoc_score = unoc_score, harc_score = harc_score)

  index = index + 1
}

res = dplyr::bind_rows(data_list)
res
#> # A tibble: 12 × 9
#>    tsk_id  cen_num ibrier_pro…¹ ibrie…² ibrie…³ ibrie…⁴ dcal_s…⁵ unoc_…⁶ harc_…⁷
#>    <chr>     <int>        <dbl>   <dbl>   <dbl>   <dbl>    <dbl>   <dbl>   <dbl>
#>  1 vet_cen     137       0       0       0       0      3.43e-30   0.5     0.5  
#>  2 vet_cen     137       0       0       0       0      3.43e-30   0.5     0.5  
#>  3 vet_cen     137       0       0       0       0      3.43e-30   0.5     0.5  
#>  4 vet_cen     137       0       0       0       0      3.43e-30   0.5     0.5  
#>  5 vet_cen     137       0       0       0       0      3.43e-30   0.5     0.5  
#>  6 vet_cen     137       0       0       0       0      8.96e-31   0.5     0.5  
#>  7 vet_obs       0       0.134   0.134   0.134   0.134  2.87e+ 1   0.739   0.739
#>  8 vet_obs       0       0.140   0.140   0.140   0.140  8.74e+ 0   0.669   0.669
#>  9 vet_obs       0       0.0904  0.0904  0.0904  0.0904 1.48e+ 1   0.671   0.671
#> 10 vet_obs       0       0.0871  0.0871  0.0871  0.0871 7.87e+ 0   0.698   0.698
#> 11 vet_obs       0       0.126   0.126   0.126   0.126  9.61e+ 0   0.671   0.671
#> 12 vet_obs       0       0.146   0.146   0.146   0.146  1.53e+ 1   0.716   0.716
#> # … with abbreviated variable names ¹​ibrier_proper_score,
#> #   ²​ibrier_proper_test_score, ³​ibrier_improper_score,
#> #   ⁴​ibrier_improper_test_score, ⁵​dcal_score, ⁶​unoc_score, ⁷​harc_score

^{Created on 2023-02-01 with reprex v2.0.2}

• For a fully censored dataset, ibrier and dcal scores give perfect prediction score (the ranger learner returns survival equal to 1 for all timepoints) - is that sensible though? Uno's and Harrell's C-indexes are random (0.5)
• For a dataset with no censoring, Uno's and Harrell's C-indexes are equivalent (as they should be I guess)
• No difference between ibrier score with proper = TRUE vs proper = FALSE

[RaphaelS1]

See response in previous chat. If everyone is censored then it's impossible to train a model, in fact the whole survival paradigm collapses. If no one is censored then it's a (probabilistic) regression problem and survival analysis doesn't really make sense

[bblodfon]

I always try to think in terms of spectrums: e.g. here going from fully observed events to fully censored. In the former case (all events observed, a regression problem as you said), I would have expected using survival analysis/models (which provide a probabilistic/distribution prediction) to be more informative versus simple regression models (where prediction is like "this patient will die in X days"). But maybe it is the case that survival analysis models are not constructed with this type of thinking in mind (e.g. a Cox model is like linear regression but not exactly, so it cannot "regress" to a simple linear regression when there is absolutely no censoring). All I am trying to say is that it might be beneficial to make such models and think in those terms.

Either way, it was nice to see that the above measures actually provide meaningful results even in these extreme case scenarios.

[RaphaelS1]

simple regression models

To be pedantic, you're actually referring to deterministic regression, which returns a single point prediction, but you can also use probabilistic regression. Now this raises some questions about the status quo of modelling, which is that we've got to a point where deterministic case is more common despite probabilistic being more informative. Yet in classification it's the reverse (most people know that $p(Y = 1) = x\%$ is better than $p(Y) = 1$. So the issue here goes beyond survival. Yes you could use survival here but it shouldn't be optimal -- so tbh a deterministic regr model could yield better overall predictions (even if less informative).

All I am trying to say is that it might be beneficial to make such models and think in those terms.

prob regr ;)

And we have this in mlr3proba, e.g.:

library(mlr3verse)
library(mlr3proba)
library(ggplot2)

# create our pipeline
l = as_learner(ppl("probregrcompositor",
  learner = lrn("regr.ranger"),
  learner_se = lrn("regr.featureless"),
  dist = "Normal")
)

# train and predict
task = tsk("mtcars")
split = partition(task, ratio = 0.9)
l$train(task, split$train)
p = l$predict(task, split$test)

[bblodfon]

I see, that's super nice! Is there a metric then to assess the calibration of the composed distr output in this case?

[RaphaelS1]

Yes but need to be fixed. Also we need to add more sophisticated learners

[RaphaelS1]

Update: now done first part babf866