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:



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.

  • $\begingroup$ As an aside, split-sample validation is not reliable unless N > 20,000. It's too dependent on the luck of the split, so is volatile. $\endgroup$ Oct 29, 2021 at 13:11
  • $\begingroup$ I have not had experience with survival::concordance with newdata. I hope to look at that soon. Your usage of rms is correct. $\endgroup$ Mar 4, 2022 at 11:54

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.