计算 c-index 以使用 R 对 Cox PH 模型进行外部验证

机器算法验证 r 回归 cox模型 验证 有效值
2022-03-28 11:33:54

首先，我要声明我知道有很多关于 c-index 的问题。我已经搜索了这个网站和其他网站，但我还没有找到适合我情况的答案。我可以成功地validate()在rms包中使用来计算我的引导内部验证的 Dxy 和 c-index。现在我需要一个 c-index 来进行外部验证。

我刚刚使用独立数据集从外部验证了我的模型，并使用val.surv(). 不幸的是，我提交的摘要没有数字，所以我不得不报告一个 c-index 以进行外部验证。我已经搜索了这个站点和 R 帮助档案，但还没有找到关于如何计算外部验证的 c-index 的结论性答案。

rcorr.cens()我已经看到在包中提到使用Hmisc，但在我看来，您只能将它用于单个变量的一致性，而不是整个模型。到目前为止，我还找不到val.surv()用于计算 c-index 的方法。我在下面发布了一些示例代码，包括一个类似于外部验证集的测试数据集。

我非常感谢您在使用独立数据集从外部验证计算 c-index 方面的帮助。

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!
3个回答

我刚刚从一位同事那里得到了有关如何使用此功能的解释。在 的帮助页面中rcorr.cens()，它指出这x是一个“数字预测变量”。我认为这意味着它必须是一个模型变量，如年龄、阶段、转移等。我发现它x可以只是你的模型对外部数据集的生存估计的向量。因此，唯一rcorr.cens()需要的两件事是生存估计向量和一个Surv()对象。使用我上面的代码，这就是你使用它的方式：

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)

我希望这对将来有同样问题的人有所帮助！

有一个包，它是 bioconductor 的一个组件，可以帮助您计算 c-index：survcomp

如果不包括生存数据，则 survcomp 中的 cindex 与您从 ROC 曲线获得的 AUC 基本相同。

我认为@JJM 提供的代码可以进行一些更改。@Seanosapien 提出的观点可以通过以下编辑的代码来解决。

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)

我所做的更改Surv()是为测试数据创建对象Surv.obj_test。这将允许与测试数据rcorr.cens进行比较。estimatestest_dat

estimates使用该模型获得生存估计的能力cox_mod。这些估计值与test_dat“时间”和“状态”变量使用rcorr.cens.

其它你可能感兴趣的问题
上一篇BUGS/JAGS 中的细化链 下一篇解释置信区间。