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.
survival::concordance
withnewdata
. I hope to look at that soon. Your usage ofrms
is correct. $\endgroup$