I am using a clinical dataset (n=624) to estimate survival via the following stratified cox model to overcome the violation of PH assumption (‘Age’ and ‘Steps’ are continuous the rest are factor variables):

cox_avoid_overfit <- survival::coxph(obj_surv ~ Age + Sex + simplerECOG + strata(PrimaryRiskGroup) + LRD + identity(Steps/5000) , data = data_avoid_overfit)

Calculating concordance as one of the performance metrics I get:

• Apparent C-Index of 0.821 [via survival::concordance(cox_avoid_overfit)]
• Perplexingly a higher optimism-corrected C-Index of 0.832 (95%CI 0.816 − 0.849; via 5-times repeated 10-fold CV using tidymodels)

When using the same model except the stratification I expectedly get a lower optimism-corrected than apparent value:

• Apparent C-Index: 0.886
• Optimism-corrected C-Index: 0.876 (95%CI 0.859 − 0.892)

I tried to reproducibly emulate this behaviour using the lung data set in the survival package using the same tidymodels approach (see code and sessionInfo() below):

1. Compared to the non-stratified model where expectedly the apparent C-Index is higher (0.647) than the optimism-corrected C-Index (0.614; using 10-fold & 5 repeated CV).
2. When stratifying by ECOG performance status the apparent and optimism-corrected are the same when rounding to 3 digits (0.581). I suspect this may be a less extreme phenomenon of what I saw with my data set…
3. Could this be due to what it say says in the tidymodels/parsnip documentation (.html) in the ‘linear predictor’ paragraph? And if so, where do I set the ‘increasing’ argument.
4. If not, any pointers to other reasons +/- fixes?

Thanks for the amazing tidymodels ecosystem. I appreciate any pointers. Thanks, Thilo

library(tidyverse) 
library(tidylog)
library(survival)
library(tidymodels)
library(censored)

# Data prep and tidymodels setup ----

lung %>% glimpse()
lung_clean <- lung %>% 
  mutate(sex=factor(sex, levels=c(1,2), labels=c("M", "F")),
         ps.ecog=factor(ph.ecog, levels=c(0,1,2,3), labels=c("ECOG 0", "ECOG 1","ECOG 2","ECOG 3")),
         had_event=ifelse(status==1, 0, 1),
         time_months = interval(today(), today() + time) / months(1)) %>%
        tidylog::drop_na() %>% 
        mutate(obj_surv = Surv(time_months, had_event), .keep = "all") %>%
  select(-status, -ph.ecog, -time)
lung_clean %>% glimpse()

set.seed(321)
lung_folds <- rsample::vfold_cv(lung_clean, v = 10, repeats = 5, strata = "had_event")
lung_folds

survival_metrics <- metric_set(brier_survival_integrated, 
                               brier_survival, 
                               roc_auc_survival, 
                               concordance_survival)
evaluation_time_points <- seq(6, 24, 6) # i.e. 6, 12, 18, 24

survreg_spec <- survival_reg() %>% 
  set_engine("survival") %>% 
  set_mode("censored regression")

# NON-stratified cox model ----

cox_lung = coxph(obj_surv ~ age + sex + ps.ecog + wt.loss, data=lung_clean)
print(cox_lung)
broom::tidy(cox_lung, exponentiate = TRUE)

test.ph <- cox.zph(cox_lung)
print(test.ph)

## apparent/optimistic metrics ----

concordance(cox_lung)
broom::glance(cox_lung) %>% pivot_longer(everything(), names_to = "metric", values_to = "value")

## optimism-corrected metrics ----
# (using CV resampling via tidymodels)

rec_non_strata <- recipe(obj_surv ~ age + sex + ps.ecog + wt.loss, data = lung_clean) 

survreg_wflow_non_strata <- workflow() %>% 
  add_recipe(rec_non_strata) %>% 
  add_model(survreg_spec)

set.seed(111)
survreg_res_non_strata <- fit_resamples(
  survreg_wflow_non_strata,
  resamples = lung_folds,
  metrics = survival_metrics,
  eval_time = evaluation_time_points, 
  control = control_resamples(save_pred = TRUE)
)
survreg_res_non_strata

preds_non_strata <- collect_predictions(survreg_res_non_strata)

roc_curves_data_non_strata <- roc_curve_survival(
  preds_non_strata,
  truth = obj_surv,
  .pred
)
workflowsets::autoplot(roc_curves_data_non_strata)

C_no_steps <- collect_metrics(survreg_res_non_strata) %>% 
  filter(.metric == "concordance_survival") %>% 
  mutate(
    ci_lower = mean - 1.96 * std_err,
    ci_upper = mean + 1.96 * std_err
  ) %>%
  print()


# Stratified cox model ----

cox_lung_stratified = coxph(Surv(time_months, had_event) ~ age + sex + strata(ps.ecog) + wt.loss, data=lung_clean)
print(cox_lung_stratified)

concordance(cox_lung_stratified)
broom::glance(cox_lung_stratified) %>% pivot_longer(everything(), names_to = "metric", values_to = "value")

## optimism-corrected metrics ----

model_formula_strata <- obj_surv ~ age + sex + strata(ps.ecog) + wt.loss
preproc_formula_strata <- obj_surv ~ age + sex + ps.ecog + wt.loss

survreg_wflow_strata <- workflow() %>% 
  add_formula(preproc_formula_strata) %>% 
  add_model(survreg_spec, formula = model_formula_strata)

set.seed(222)
survreg_res_strata <- fit_resamples(
  survreg_wflow_strata,
  resamples = lung_folds,
  metrics = survival_metrics,
  eval_time = evaluation_time_points, 
  control = control_resamples(save_pred = TRUE)
)
survreg_res_strata


preds_strata <- collect_predictions(survreg_res_strata)
#preds_strata
#preds_strata$.pred[[1]]
#preds_strata %>% count(id)


roc_curves_data_strata <- roc_curve_survival(
  preds_strata,
  truth = obj_surv,
  .pred
)
workflowsets::autoplot(roc_curves_data_strata)

C_strata <- collect_metrics(survreg_res_strata) %>% 
  filter(.metric == "concordance_survival") %>% 
  mutate(
    ci_lower = mean - 1.96 * std_err,
    ci_upper = mean + 1.96 * std_err
  ) %>%
  print()

options(digits=7) # temp change to see more decimals
print.data.frame(C_strata)
#print.data.frame(collect_metrics(survreg_res_strata, summarize = FALSE) %>% filter(.metric=="concordance_survival") %>% mutate(mean_concordance = mean(.estimate)))
options(digits=3)

sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.6.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Sydney
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] censored_0.3.2     yardstick_1.3.1    workflowsets_1.1.0 workflows_1.1.4    tune_1.2.1         rsample_1.2.1      recipes_1.1.0      parsnip_1.2.1     
 [9] modeldata_1.4.0    infer_1.0.7        dials_1.3.0        scales_1.3.0       broom_1.0.7        tidymodels_1.2.0   survival_3.8-3     tidylog_1.1.0     
[17] lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
[25] ggplot2_3.5.1      tidyverse_2.0.0   

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3  naniar_1.1.0        rstudioapi_0.17.1   shape_1.4.6.1       magrittr_2.0.3      jomo_2.7-6          magick_2.8.4       
  [8] corrplot_0.95       nloptr_2.1.1        farver_2.1.2        rmarkdown_2.29      vctrs_0.6.5         bstfun_0.5.1.9003   minqa_1.2.8        
 [15] base64enc_0.1-3     blogdown_1.19       htmltools_0.5.8.1   finalfit_1.0.8      cellranger_1.1.0    Formula_1.2-5       mitml_0.4-5        
 [22] parallelly_1.38.0   plyr_1.8.9          cobalt_4.5.5        gt_0.11.1           MatchIt_4.5.5       lifecycle_1.0.4     iterators_1.0.14   
 [29] pkgconfig_2.0.3     Matrix_1.7-0        R6_2.5.1            fastmap_1.2.0       future_1.34.0       digest_0.6.37       GGally_2.2.1       
 [36] colorspace_2.1-1    lobstr_1.1.2        furrr_0.3.1         patchwork_1.2.0     rprojroot_2.0.4     gtsummary_2.0.3     labeling_0.4.3     
 [43] fansi_1.0.6         timechange_0.3.0    abind_1.4-8         compiler_4.4.1      here_1.0.1          withr_3.0.2         pander_0.6.5       
 [50] backports_1.5.0     carData_3.0-5       DBI_1.2.3           ggstats_0.7.0       performance_0.12.3  R.utils_2.12.3      pan_1.9            
 [57] MASS_7.3-60.2       lava_1.8.0          tsibble_1.1.5       ggsci_3.2.0         tools_4.4.1         visdat_0.6.0        future.apply_1.11.3
 [64] nnet_7.3-19         R.oo_1.27.0         glue_1.8.0          nlme_3.1-164        modelenv_0.2.0      grid_4.4.1          checkmate_2.3.2    
 [71] reshape2_1.4.4      generics_0.1.3      gtable_0.3.6        tzdb_0.4.0          R.methodsS3_1.8.2   class_7.3-22        data.table_1.16.2  
 [78] hms_1.1.3           car_3.1-3           xml2_1.3.6          utf8_1.2.4          foreach_1.5.2       pillar_1.9.0        splines_4.4.1      
 [85] lhs_1.2.0           pryr_0.1.6          lattice_0.22-6      tidyselect_1.2.1    knitr_1.49          gridExtra_2.3       xfun_0.49          
 [92] hardhat_1.4.0       rapportools_1.1     timeDate_4041.110   matrixStats_1.3.0   stringi_1.8.4       DiceDesign_1.10     yaml_2.3.10        
 [99] boot_1.3-30         ggokabeito_0.1.0    evaluate_1.0.1      codetools_0.2-20    tcltk_4.4.1         cli_3.6.3           rpart_4.1.23       
[106] munsell_0.5.1       Rcpp_1.0.13-1       readxl_1.4.3        globals_0.16.3      anytime_0.3.9       summarytools_1.0.1  parallel_4.4.1     
[113] ellipsis_0.3.2      gower_1.0.1         GPfit_1.0-8         lme4_1.1-35.5       listenv_0.9.1       glmnet_4.1-8        ipred_0.9-15       
[120] prodlim_2024.06.25  ggridges_0.5.6      insight_0.20.5      crayon_1.5.3        clisymbols_1.2.0    rlang_1.1.4         mice_3.16.0