Cox PH - R 'rms' package - Concordance (Train, CV & Test)

Question

Please could someone verify if my methodology for calculating these three versions of concordance are correct when using the rms package? I have really struggled to find many examples for performing cross-validation & test concordance, in particular.

I am working on a dataset of ~550k rows at work, but please see the reproducible example below:

library(rms)
library(caret)

set.seed(1)

index <- createDataPartition(veteran$status, p = 0.6, list = F)
veteran_train <- veteran[index, ] # 60%
veteran_test <- veteran[-index, ] # 40%

cox <- cph(Surv(time, status) ~ karno + celltype, 
           data = veteran_train, 
           x = T, 
           y = T)

Using the relationship: $C = \frac{D_{xy} + 1}{2}$

'Train' Concordance:

train_C <- as.numeric((cox$stats[9] + 1)/2)
round(train_C, 4) # 0.7236

'Train' (10-fold CV) Concordance:

cv_version <- validate(cox, method = "crossvalidation", B = 10)
train_CV_C <- (cv_version[1, 3] + 1)/2
round(train_CV_C, 4) # 0.7167

'Test' Concordance:

actuals <- Surv(veteran_test$time, veteran_test$status)
estimates <- survest(cox, newdata = veteran_test, times = 1337)$surv # can pick an arbitrary time here, same results
test_C <- as.numeric(rcorr.cens(x = estimates, S = actuals)[1])
round(test_C, 4) # 0.6977

I investigated the rms package after attempting to find the test concordance through the survival package.

However I think (please confirm) survival::concordance() is bugged when specifying the newdata argument, only running if veteran_test has the exact same number of rows as veteran_train:

cox <- coxph(Surv(time, status) ~ karno + celltype, data = veteran_train)
# test_C_v2 <- survival::concordance(cox, newdata = veteran_test)$concordance

This returns: Error... x and y are not the same length.

If I then create a situation where veteran_train & veteran_test do have the same number of rows, I get different test concordance results to the rms package!:

veteran_train <- veteran_train[1:nrow(veteran_test), ]
nrow(veteran_train) == nrow(veteran_test) # 41


# 'rms' method:
cox_rms <- rms::cph(Surv(time, status) ~ karno + celltype, 
           data = veteran_train, 
           x = T, 
           y = T)

actuals <- Surv(veteran_test$time, veteran_test$status)
estimates <- survest(cox_rms, newdata = veteran_test, times = 1337)$surv
test_C_v1 <- as.numeric(rcorr.cens(x = estimates, S = actuals)[1])
round(test_C_v1, 4) # 0.6977


# 'survival' method
cox_survival <- survival::coxph(Surv(time, status) ~ karno + celltype, data = veteran_train)

test_C_v2 <- survival::concordance(cox_survival, newdata = veteran_test)$concordance
round(test_C_v2, 4) # 0.4333

Is one of these approaches correct? When applying this methodology to my data, there was a common theme that the train/CV/test Concordance were all very close when using the rms package.

However, when using the survival package, the test concordance would see a significant drop (e.g. 0.7 train -> 0.55 test), which doesn't seem believable for a model with just one predictor and a large volume of data.

Answer 1

Your RMS method solution agrees with the solution that Frank Harrell approved in this answer. I've found survest to be extremely slow, particularly when running predictions for >10 bootstrap re-samples of test folds. But aside from code performance, the accuracy of the calculations should be more robust. HTH.