如何在R软件中求一致性指数( Harrell concordance index:C-index)？

所谓C-index，英文名全称concordance index，中文里有人翻译成一致性指数，最早是由范德堡大学（Vanderbilt University）生物统计教教授Frank E Harrell Jr 1996年提出，主要用于计算生存分析中的COX模型预测值与真实之间的区分度（discrimination），也称为Harrell's concordance index ；现阶段用的最多的是肿瘤患者预后模型的预测精度。一般评价模型的好坏主要有两个方面，一是模型的拟合优度（Goodness of Fit),常见的评价指标主要有R方，-2logL,AIC,BIC等等；另外一个是模型的预测精度，主要就是模型的真实值与预测值之间的差的大小，均方误差，相对误差等。从临床应用的角度来说，我们更注重后者，即统计建模主要是用于预测，而从C-index的概念大家看出它属于模型评价指标的后者，这一指标比前面提到的几个指标看起来更高大上，一般文献中用的也比较多。换句话说，如果预后模型建好，效果不错，即使不知道如何计算C-index值，报告软件输出结果中的预测误差是相同效果，再添加拟合优度会更能说明效果，这样反而更实用。

C-index本质上是估计了预测结果与实际观察到的结果相一致的概率，即资料所有病人对子中预测结果与实际结果一致的对子所占的比例。有点类似于ROC曲线下面积。

C-index的计算方法是:把所研究的资料中的所有研究对象随机地两两组成对子。以生存分析为例,对于一对病人,如果生存时间较长的一位,其预测生存时间长于生存时间较短的一位,或预测的生存概率高的一位的生存时间长于生存概率低的另一位,则称之为预测结果与实际结果一致。

C-index的计算步骤为:

(1)产生所有的病例配对。若有n个观察个体,则所有的对子数应为Cn2(组合数)?

(2)排除下面两种对子:对子中具有较小观察时间的个体没有达到观察终点及对子中两个个体都没达到观察终点。剩余的为有用对子。

(3)计算有用对子中,预测结果和实际相一致的对子数,即具有较坏预测结果个体的实际观察时间较短。

(4)计算。C=一致对子数/有用对子数。

由上述计算方法可以看出,C-index在0.5-1之间。0.5为完全不一致,说明该模型没有预测作用,1为完全一致,说明该模型预测结果与实际完全一致。在实际应用中,很难找到完全一致的预测模型,既往研究认为,C-index在0.50-0.70为较低准确度:在0.71-0.90之间为中等准确度;而高于0.90则为高准确度。

当C-index检验由同一样本建成的模型时易造成偏倚,因此再采用重抽样技术(Bootstrap)可以几乎无偏倚的检验预测模型的准确度。Bootstrap是非参数统计中一种重要的估计统计量方差进而进行区间估计的统计方法,是现代统计应用较为广泛的一种统计方法。

Bootstrap方法核心思想和基本步骤如下:

(1)采用重抽样技术从原始样本中抽取一定数量的样本,此过程允许重复抽样。

(2)根据抽出的样本计算给定的统计量T。

(3)重复上述N次(一般大于1000),得到N个统计量T。

(4)计算上述N个统计量T的样木方差,得到统计量的方差。

Bootstarap方法只是对单一样本且样本量较小的资料，如果数据集很大可以按照不同的比例将数据集拆分，一部分用于建模一部分用于验证。关于交叉验证（Cross-validation），由于篇幅有限，留作下次探讨。

R软件实现：

C-index的R软件计算实现有两种实现方法，一种是用到Harrell本人的的R包Hmisc package；另一种是Le Kang, Weijie Chen 2014年12月18日发布的R compareC Package；

############################

#### Method 1.Hmisc code ###

############################

data <- read.csv("survivaldta.csv") ###读入csv格式数据####

library(Hmisc) ###加载Hmisc包，前提是安装了，如果没有安装，请百度如果安装R package，这里就不详细讲了！

library(survival) ###加载survival包，主要用于建立模型###

f <- cph(Surv(time,death)~x1+x2+x3，data=survivldata) ###拟合cox模型

fp <- predict(f)###模型的预测值

cindex.orig=1-rcorr.cens(fp,Surv(time,death)) [[1]]###计算出的C-index

###############################

#### Method 2.compareC code ###

###############################

data <- read.csv("survivaldta.csv") ###读入csv格式数据####

library(compareC) ###加载compareC包，前提是安装了，如果没有安装，请百度如果安装R package，这里就不详细讲了！

library(survival) ###加载survival包，主要用于建立模型###

cindex <- cindex(Surv(time,death) ~ x1+x2+x3,data=survivldata)###计算出的C-index

###############################

#### Bootstrap code ###

###############################

bootit=200

for(i in 1:bootit){

case=noNA[group=="long",] 

control=noNA[group=="<24",]

bootindex.case=sample(1:nrow(case),replace=T)

boot.case.data=case[bootindex.case,]

bootindex.control=sample(1:nrow(control),replace=T)

boot.control.data=control[bootindex.control,]

boot.data=rbind(boot.case.data,boot.control.data)

dstr.boot=svydesign(id=~1, prob=~inv_weight, fpc=~ssize, data=boot.data)

boot.fit=svycoxph(Surv(survival,surv_cens) ~x1+x2+x3,data=boot.data,x=TRUE,design=dstr.boot)

cindex.train=1-rcorr.cens(lp.boot,Surv(boot.data$survival, boot.data$surv_cens))[[1]]

cindex.test=1-rcorr.cens(lp_=.test,Surv(noNA$survival,noNA$surv_cens))[[1]]

bias=rep(1,bootit)

bias[i]=abs(cindex.train-cindex.test) }

参考文献：

Harrell FE, Califf RM, Pryor DB, Lee KL, and Rosati RA. (1982) Evaluating the yield of medical tests. The Journal of the American Medical Association, 247(18), 2543–2546

Pencina MJ and D’Agostino RB. (2004) Overall C as a measure of discrimination in survival analysis: model specific population value and confidence interval estimation. Statistics in Medicine, 23(13), 2109–2123

Kang L, Chen W, Petrick NA, and Gallas BD. (2014) Comparing two correlated C indices with right-censored survival outcome: a one-shot nonparametric approach. Statistics in Medicine, 34(4), 685–703, doi: 10.1002/sim.6370

Hmisc Reference manual：

compare C： Reference manual

Frank.Harrell