6
$\begingroup$

First off, I'll state that I'm aware many questions get asked about the c-index. I've searched this site and others, and I haven't found an answer for my situation. I can successfully use validate() in the rms package to calculate the Dxy and c-index for my boot-strapped internal validation. Now I need a c-index for my external validation.

I just externally validated my model with an independent data set, and produced the desired plot of actual vs. predicted probabilities using val.surv(). Unfortunately the abstract I'm submitting cannot have figures, so I feel compelled to report a c-index for the external validation. I have searched this site and the R help archive and haven't found a conclusive answer for how to calculate the c-index on external validation.

I've seen mentions of using rcorr.cens() in the Hmisc package, but it seems to me that you can only use it for a single variable's concordance, not an entire model. As of yet I can't find a way to use val.surv() for calculating a c-index as well. I posted some sample code below, including a test data set to act like an external validation set.

I really appreciate your help in calculating a c-index from an external validation with an independent data set.

library(rms)
library(Hmisc)
data(veteran)

##Create a Cox PH model for the training data.
survmod=with(veteran,Surv(time,status))
cox.mod=cph(survmod~celltype+karno,data=veteran,x=T,y=T,surv=TRUE,time.inc=5*365)

##Here is the test data set that is the external "independent" data.
test_dat=data.frame(trt=replicate(500,NA), celltype=replicate(500,NA), time=replicate(500,NA), status=replicate(500,NA), karno=replicate(500,NA), diagtime=replicate(500,NA), age=replicate(500,NA), prior=replicate(500,NA))
for(i in seq(8)){
test_dat[,i]=sample(veteran[,i],500,replace=T)
}

##Validate the model with the test data
test_surv=with(test_dat,Surv(time,status))
validated=val.surv(cox.mod,newdata=test_dat,S=test_surv)

##Now what I need is to take the external validation and compute the 
#c-index.  This is where I'm stuck.  
#I've seen people mention `rcorr.cens()`, but I can't figure out a way to use 
#`rcorr.cens()` with a Cox model of several variables.  I appreciate your help!
Improve this question
$\endgroup$

3 Answers 3

Reset to default
9
$\begingroup$

I just got an explanation from a colleague about how to use this feature. In the help page for rcorr.cens(), it states that x is a "numeric predictor variable". I thought that this meant it had to be a model variable like Age, Stage, Metastasis, etc. What I found out is that x can just be a vector of your model's survival estimates for an external data set. Therefore the only two things rcorr.cens() needs is that vector of survival estimates and a Surv() object. Using my code from above, this is how you use it:

library(rms)
surv.obj=with(veteran,Surv(time,status))   ####This will be used for rcorr.cens
cox.mod=cph(surv.obj~celltype+karno,data=veteran,x=T,y=T,surv=TRUE,time.inc=5*365)

##Here is the test data set that is the external "independent" data.
test_dat=data.frame(trt=replicate(500,NA), celltype=replicate(500,NA), time=replicate(500,NA), status=replicate(500,NA), karno=replicate(500,NA), diagtime=replicate(500,NA), age=replicate(500,NA), prior=replicate(500,NA))
for(i in seq(8)){
test_dat[,i]=sample(veteran[,i],500,replace=T)
}

###Create your survival estimates
estimates=survest(cox.mod,newdata=test_dat,times=5*365)$surv


###Determine concordance
rcorr.cens(x=estimates,S=surv.obj)

I hope this helps anyone in the future who has the same question!

Improve this answer
$\endgroup$
2
  • 7
    $\begingroup$ This is correct. Or just use the linear predictor in the call to rcorr.cens if you do not have any strata. Take $D_{xy}$ from the output in that case, negate it, and use the relationship $D_{xy} = 2\times (C - \frac{1}{2})$. $\endgroup$
    – Frank Harrell
    Jan 23, 2013 at 20:48
  • 1
    $\begingroup$ The final line of code doesn't work. I think you need to use 137 instead of 500 so as to match initial n in veteran data. $\endgroup$
    – Seanosapien
    Nov 2, 2017 at 17:00
2
$\begingroup$

There is a package, which is a component of bioconductor that can help you calculate the c-index : survcomp

If you don't include survival data, than cindex in survcomp is basically the same as the AUC you get from the ROC curve.

Improve this answer
$\endgroup$
1
$\begingroup$

I think code provided by @JJM could work with some changes. The point raised by @Seanosapien could be addressed by following edited code.

library(rms)
surv.obj=with(veteran,Surv(time,status))   ####This will NOT be used for rcorr.cens
cox.mod=cph(surv.obj~celltype+karno,data=veteran,x=T,y=T,surv=TRUE,time.inc=5*365)

##Here is the test data set that is the external "independent" data.
test_dat=data.frame(trt=replicate(500,NA), celltype=replicate(500,NA), 
time=replicate(500,NA), status=replicate(500,NA), karno=replicate(500,NA), 
diagtime=replicate(500,NA), age=replicate(500,NA), prior=replicate(500,NA))
for(i in seq(8)){
   test_dat[,i]=sample(veteran[,i],500,replace=T)
   }

Surv.obj_test=with(test_dat,Surv(time,status)) #This will be used for rcorr.cens

###Create your survival estimates
estimates=survest(cox.mod,newdata=test_dat,times=5*365)$surv


###Determine concordance
rcorr.cens(x=estimates,S=Surv.obj_test)

The change I have made is that Surv() object is created for the test data Surv.obj_test. This will allow rcorr.cens to compare the estimates with the test data test_dat.

The estimates enable in getting survival estimates using the model cox_mod. These estimates are compared with the test_dat "time" and "status" variable using the rcorr.cens.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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