Topic 13. 临床预测模型之一致性指数 (C-index)

发布时间：2023-02-07 19:00

前言

C-index，C指数即一致性指数（concordance index），用来评价模型的预测能力。C指数是指所有病人对子中预测结果与实际结果一致的对子所占的比例。它估计了预测结果与实际观察到的结果相一致的概率。c指数的计算方法是：把所研究的资料中的所有研究对象随机地两两组成对子。以生存分析为例，对于一个病人，如果生存时间较长的一位的预测生存时间也长于另一位的预测生存时间，或预测的生存概率高的一位的生存时间长于生存概率低的另一位，则称之为预测结果与实际结果一致。

C-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值，报告软件输出结果中的预测误差是相同效果，再添加拟合优度会更能说明效果，这样反而更实用。

What is a C-Statistic?

Figure 1: The concordance statistic is equal to the area under a ROC curve.

If by the “c-statisitc” you are referring to the measure of the discriminative power of the logistic equation, you can calculate it by saving the predictive probabilities from the logistic regression analysis and running a ROC curve with the predictive probabilities as the “test variable.” The c-statistic is the area under the curve value.

评价一个预测模型的表现可以从三方面来度量：

区分能力（discrimination）：指的是模型区分有病/没病，死亡/活着等结局的预测能力。简单举个例子，比如说，现有100个人，50个有病，50个健康；你用预测模型预测出46个有病，54个没病。那么这46个覆盖到50个真正有病的人的多少就直接决定了你模型预测的靠谱程度，和准确性。通常用ROC、C指数来度量，当然NRI（Net reclassification improvement）和 IDI（integrated discrimination improvement）也是度量指标之一。
一致性（Calibration）：指结局实际发生的概率和预测的概率的一致性。读起来有点费解，我们还是举上面这个例子，我们预测的100个人，不是指我们真用模型预测出来有病/没病，模型只给我们有病的概率，根据概率大于某个截断值（比如说0.5）来判断有病/没病。100个人，我们最终通过模型得到了100个概率，也就是100个0-1之间的数，我们将这100个数，按照从小到大排列，再依次将这100个人分成10组，每组10个人，实际的概率就是这10个人中发生疾病的比例，预测的概率就是每组预测得到的10个数的平均值，然后比较这两个数，一个作为横坐标，一个作为纵坐标，就得到了一致性曲线图（当然得到95%可信区间后更完整了）。当然一致性还可以通过 Hosmer-Lemeshow goodness-of-fit test 来度量。
总体上（overall）：事实上就是综合了区分能力和一致性的度量指标，比如 R2。

在很多临床文章中经常看见统计方法里面描述模型的区分能（discrimination ability）用C指数来度量，其实准确来说，这个C指数应该指明是哪个人提出来的C指数:

C-index本质上是估计了预测结果与实际观察到的结果相一致的概率，即资料所有病人对子中预测结果与实际结果一致的对子所占的比例；分类问题中，正类预测score大于负类预测score的概率即是C指数（Mann–Whitney U检验的C统计量），也称AUC，推导可以参考：AUC与Mann–Whitney U test 以及 wikipedia/Mann_Whitney_U_test。当然，生存数据也有专门的AUC计算方法，integratedAUC（ iAUC），对应的ROC一般称之为riskROC、时间依赖ROC等。integratedAUC和C-index其实差不太多，计算的结果也相差不大。

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

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

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

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

计算，C-index = 一致对子数/有用对子数。

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

这和AUC其实是类似的，如果完全预测反了，那么C-index其实应该是1。但是这种情况仅在训练集中可以这么计算，在验证集中，预测反了就是预测反了。因此看到C-index或AUC小于0.5的情况，也不要大惊小怪。

我们从pubmed 上搜索 C-index，发现有很多临床文章均使用过一致性指数，可见这种模型评估方式是临床预测模型构建时不可或缺的分析，掌握了就等于文章不愁了！

https://pubmed.ncbi.nlm.nih.gov/

Figure 1: The C-index in papers

实例解析

我们这里利用6种方法分别计算 C-index ，比较每种方法的区别，再对模型评估时选择准确的方法十分重要。

1. 安装及加载

我们选择的这几种方法需要安装及加载5个程序包，有安装过的，可以略过，安装如下：

if (!require(survival)) {
    install.packages("survival")
}
if (!require("rms")) {
    install.packages("rms")
}
if (!require("ROCR")) {
    install.packages("ROCR")
}
## Warning: 程辑包'ROCR'是用R版本4.1.3 来建造的
if (!require("Hmisc")) {
    install.packages("Hmisc")
}
if (!require("survcomp")) {
    BiocManager::install("survcomp")
}
library(survival)
library(rms)
library(ROCR)
library(Hmisc)
library(survcomp)

2. 数据读取

我们仍然采用软件包 survival 自带的肺癌数据库 NCCTG Lung Cancer Data 作为输入数据，如下：

Descrption

Survival in patients with advanced lung cancer from the North Central Cancer Treatment Group. Performance scores rate how well the patient can perform usual daily activities.

data(package = "survival")
# 打包数据
lung$sex = factor(lung$sex)
lung = na.omit(lung)
dd <- datadist(lung)
options(datadist = "dd")
head(dd)
## $limits
##                 inst   time status  age  sex ph.ecog ph.karno pat.karno
## Low:effect         3  174.5      1 57.0        0       70        70
## Adjust to         11  268.0      1 64.0    1       1       80        80
## High:effect       15  419.5      2 70.0        1       90        90
## Low:prediction     1   53.0      1 44.3    1       0       60        60
## High:prediction   26  737.3      2 75.0    2       2      100       100
## Low                1    5.0      1 39.0    1       0       50        30
## High              32 1022.0      2 82.0    2       3      100       100
##                 meal.cal wt.loss
## Low:effect         619.0     0.0
## Adjust to          975.0     7.0
## High:effect       1162.5    15.0
## Low:prediction     320.5    -5.0
## High:prediction   1477.5    35.4
## Low                 96.0   -24.0
## High              2600.0    68.0
## 
## $values
## $values$status
## [1] 1 2
## 
## $values$sex
## [1] "1" "2"
## 
## $values$ph.ecog
## [1] 0 1 2 3
## 
## $values$ph.karno
## [1]  50  60  70  80  90 100
## 
## $values$pat.karno
## [1]  30  40  50  60  70  80  90 100

3. 模型构建

我们根据数据情况构建 Logistic 回归模型，当然我们之前已经通过各种建模的方式进行比较，其中显著的几个变量为 sex 和 ph.ecog，次之的为 age 和 ph.karno，那么我们就通过这四个变量进行 Logistic 回归建模，如下：

cph <- cph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno, data = lung, x = TRUE,
    y = TRUE, surv = TRUE)
cph
## Cox Proportional Hazards Model
##  
##  cph(formula = Surv(time, status) ~ age + sex + ph.ecog + ph.karno, 
##      data = lung, x = TRUE, y = TRUE, surv = TRUE)
##  
##                         Model Tests    Discrimination    
##                                               Indexes    
##  Obs       167    LR chi2     23.46    R2       0.131    
##  Events    120    d.f.            4    Dxy      0.273    
##  Center 2.9977    Pr(> chi2) 0.0001    g        0.557    
##                   Score chi2  23.50    gr       1.745    
##                   Pr(> chi2) 0.0001                      
##  
##           Coef    S.E.   Wald Z Pr(>|Z|)
##  age       0.0129 0.0114  1.13  0.2586  
##  sex=2    -0.5239 0.1978 -2.65  0.0081  
##  ph.ecog   0.7426 0.2086  3.56  0.0004  
##  ph.karno  0.0205 0.0113  1.81  0.0699  
##

使用survival程序包中的 coxph 函数构造 Cox 回归模型，如下：

coxph <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno, data = lung)
coxph
## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.ecog + ph.karno, 
##     data = lung)
## 
##              coef exp(coef) se(coef)      z        p
## age       0.01290   1.01298  0.01142  1.129 0.258911
## sex2     -0.52418   0.59204  0.19777 -2.650 0.008038
## ph.ecog   0.74274   2.10169  0.20863  3.560 0.000371
## ph.karno  0.02048   1.02069  0.01130  1.813 0.069796
## 
## Likelihood ratio test=23.46  on 4 df, p=0.0001026
## n= 167, number of events= 120

4. C-index 计算

这里我从网络和统计书籍中共整理总结了6种方法可以很快的计算出 C-index，希望能够满足大家的需求。

1. validate {rms}

利用rms包中的cph函数和validate函数，可提供un-adjusted和bias adjusted C指数两种，这种方法我觉得最大的优势就是给出了校正的C指数，但是都没有95%可信区间。如下：

# Get the Dxy
v <- validate(cph, dxy = TRUE, B = 1000)
Dxy = v[rownames(v) == "Dxy", colnames(v) == "index.corrected"]
orig_Dxy = v[rownames(v) == "Dxy", colnames(v) == "index.orig"]

# The c-statistic according to Dxy=2(c-0.5)
bias_corrected_c_index <- abs(Dxy)/2 + 0.5
orig_c_index <- abs(orig_Dxy)/2 + 0.5

bias_corrected_c_index
## [1] 0.6195689
orig_c_index
## [1] 0.6364067

2. summary$concordance {survival}

直接从survival包的函数coxph结果中输出，可以看出这种方法输出了C指数，也输出了标准误，那么95%可信区间就可以通过加减1.96*se得到。并且这种方法也适用于很多指标联合。

sum.surv <- summary(coxph)
c_index <- sum.surv$concordance
c_index
##          C      se(C) 
## 0.63640666 0.03046212

3. survConcordance {survival}

利用survival 软件包中的函数survConcordance，这种方法和方法2类似，如下：

concordance(Surv(time, status) ~ predict(coxph), lung)
## Call:
## concordance.formula(object = Surv(time, status) ~ predict(coxph), 
##     data = lung)
## 
## n= 167 
## Concordance= 0.3636 se= 0.03046
## concordant discordant     tied.x     tied.y    tied.xy 
##       3830       6712         22         11          0

4. concordance.index {survcomp}

这种方法的优点就是可以直接输出95%可信区间，不需要自己再进行计算，如下：

cindex <- concordance.index(predict(coxph), surv.time = lung$time, surv.event = lung$status)
cindex$c.index
## [1] 0.4335472
cindex$lower
## [1] 0.2439767
cindex$upper
## [1] 0.6447896

5. performance {ROCR}

利用软件包 ROCR 中的函数 prediction来求出预测值，在进行 performance 计算假阳和假阴性，求出 AUC 即为 C-index，如下：

Figure 3: The C-index is equal to the area under a ROC curve.

6. rcorr.cens {Hmisc}

rcorr.cens的代码及结果，第一个值就是C指数，同时也有Dxy的值，注意rcorrcens的写法是写成formula（公式）的形式，较为方便；而rcorr.cens的写法是只能在前面写上一个自变量,并且不支持data = ...的写法，有点繁琐。较为遗憾的是这种方法得到的C指数的标准误需要通过S.D./2间接得到。

r <- rcorrcens(Surv(time, status) ~ age + sex + ph.ecog + ph.karno, data = lung)
r
## 
## Somers' Rank Correlation for Censored Data    Response variable:Surv(time, status)
## 
##              C    Dxy  aDxy    SD    Z      P   n
## age      0.441 -0.118 0.118 0.060 1.98 0.0475 167
## sex      0.567  0.135 0.135 0.049 2.73 0.0063 167
## ph.ecog  0.385 -0.229 0.229 0.056 4.08 0.0000 167
## ph.karno 0.598  0.195 0.195 0.062 3.17 0.0015 167

这期总结了6种方法可以很快地计算 C-index，方便大家调用，在使用的过程中选一种即可。

References:

Harrell FE Jr, Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Stat Med. 1996;15(4):361-387.
Loprinzi CL. Laurie JA. Wieand HS. Krook JE. Novotny PJ. Kugler JW. Bartel J. Law M. Bateman M. Klatt NE. et al. Prospective evaluation of prognostic variables from patient-completed questionnaires. North Central Cancer Treatment Group. Journal of Clinical Oncology. 12(3):601-7, 1994.
Mao S, Yu X, Shan Y, Fan R, Wu S, Lu C. Albumin-Bilirubin (ALBI) and Monocyte to Lymphocyte Ratio (MLR)-Based Nomogram Model to Predict Tumor Recurrence of AFP-Negative Hepatocellular Carcinoma. J Hepatocell Carcinoma. 2021;8:1355-1365.