临床建模10-Cox回归中C-index的两种常用计算方法

临床建模10-Cox回归中C-index的两种常用计算方法

题记：今天这篇文章是预测模型系列的第5篇文章，我们主要介绍基于R语言的Cox回归中C-index的两种常用计算方法言，请各位斧正~

1. 背景介绍

近十年来，临床研究中有一类预测模型构建与验证类的文章数量逐渐增多。简言之，预测模型是通过已知参数来预测临床未知的结局，而模型本身就是一个数学公式。也就是把已知的参数通过这个所谓的模型就可以计算出未知的结局发生的可能性，即是预测。

临床预测模型的统计学本质就是回归建模分析，回归的本质就是发现因变量Y与多个自变量X之间的数学关系。临床研究中最常用的回归建模分析有三种：多重线性回归、Logistic回归与Cox回归。当我们通过训练集构建了一个回归模型，我们如何去科学评价一个回归模型预测的准确与否呢？举个例子，有两位算命的半仙儿在街角各支了一个算命的摊子，作为求姻缘的王家大小姐是选择张半仙儿算命还是选择李半仙儿呢？一个简单的选择方法：谁算命准就找谁算！临床预测模型大抵如此，基本的要求就是要“算命准”。总体来讲评价一个预测模型的优劣可以从三方面来度量：

（1）区分能力（Discrimination）：指的是回归模型区分有病/无病、有效/无效、死亡/存活等结局的预测能力。简单举个例子，比如说，现有100个人，50个确定患病，50个不患病；我们用预测模型预测出45个有病，55个没病。那么这45个覆盖到50个真正有病的人的多少就直接决定了你模型预测能力的准确程度，我们称准确性。通常用ROC、C-Statistics来度量（在Logistic回归模型中ROC曲线下面积=C-Statistics），当然NRI（Net reclassification improvement）和 IDI（integrated discrimination

improvement）也是度量指标之一。

（2）一致性（Calibration）：指结局实际发生的概率和预测的概率的一致性。看起来有点让人费解，我们还是举上面这个例子，我们预测的100个人，不是指我们真用模型预测出来有病/无病，模型只给我们有病的概率，根据概率大于某个截断值（比如说0.5）来判断有病/无病。100个人，我们最终通过模型得到了100个概率，也就是100个0-1之间的数，我们将这100个数，按照从小到大排列，再依次将这100个人分成10组，每组10个人，实际的概率其实就是这10个人中发生疾病的比例，预测的概率就是每组预测得到的10个比例的平均值，然后比较这两个数，一个作为横坐标，一个作为纵坐标，就得到了一致性曲线图（Calibration plot），也可计算这个图的95%区间范围。在Logistic回归模型中，有时一致性也可以通过拟合优度检验（Hosmer-Lemeshow goodness-of-fit test）来度量。

（3）R^2 是总体上（Overall）评估模型优劣的参数，事实上就是综合了区分能力和一致性的度量指标。模型决定系数更综合，但略显粗糙。

2. C-index计算两种方法

在很多临床文章中经常看见统计方法里面描述模型的区分能力（Discrimination ability）用C-Statistics来度量，下面我们就用R语言为大家演示这个所谓的C-Statistics如何计算？本文先讲解Cox回归中的C-Statistics（一般称为C-index）的计算，Logistic回归C-Statistics计算将在后续文章中介绍。严格说来C-index包括以下几种，我们仅介绍临床上较为常用的第一种。

（1）Harrell’s C

（2）C-statistic by Begg et al.(survAUC::BeggC)

（3）C-statistic by Uno et al.(survC1::Inf.Cval; survAUC::UnoC)

（4）Gonen and Heller Concordance Index forCox models (survAUC::GHCI, CPE::phcpe, clinfun::coxphCPE)

方法1：直接从survival包的函数coxph结果中输出，需要R的版本高于2.15.需要提前安装survival包可以看出这种方法输出了C-index (对应模型参数C)，也输出了标准误，95%可信区间就可以通过C加减1.96*se得到。并且这种方法也适用于很多指标联合。

方法2：利用rms包中的cph函数和validate函数，可提供un-adjusted和bias adjusted C指数两种。

3. R语言代码及解读

模拟一组数据并设置为数据框结构，并展示数据框sample.data的前6行：

age <- rnorm(200,50,5) bp <- rnorm(200,120,10) d.time <- rexp(200) cens <- runif(200,.5,2)

death <- d.time <= cens os <- pmin(d.time, cens) sample.data <- data.frame(age =age,bp = bp,os = os,death = death) head(sample.data) ## age bp os death ## 1 51.17440 122.43753 1.43170878 FALSE ## 2 55.31156 117.43701 0.56842712 FALSE ## 3 54.81181 130.29056 0.07273923 TRUE ## 4 39.92118 107.40071 0.08936517 TRUE ## 5 46.05266 119.96172 1.46515740 FALSE ## 6 50.19414 99.39102 0.89704829 TRUE 方法1. {survival}包，载入survival包，coxph函数拟合cox回归模型，summary函数展示模型结果并赋值给对象sum.surv，所展示得模型参数concordance，即是c-index

library(survival) fit <- coxph(Surv(os, death) ~ age + bp,data = sample.data) sum.surv<- summary(fit) c_index <- sum.surv$concordance c_index ## C se(C) ## 0.50747443 0.02838136 方法2. {rms}包，模型参数 Dxy*0.5+0.5 即是c-index

library(rms) ## Loading required package: Hmisc ## Loading required package: lattice ## Loading required package: Formula ## Loading required package: ggplot2

## Attaching package: 'Hmisc' ## The following objects are masked from 'package:base': ## ## format.pval, units ## Loading required package: SparseM ## Attaching package: 'SparseM' ## The following object is masked from 'package:base': ## ## backsolve set.seed(1) #这里设置种子，目的是为了能重复最后的结果，因为validate函数的校正结果是随机的。 dd<- datadist(sample.data) options(datadist='dd') fit.cph <- cph(Surv(os, death)~ age + bp, data = sample.data, x = TRUE, y = T

RUE, surv = TRUE) fit.cph ## Cox Proportional Hazards Model ## ## cph(formula = Surv(os, death) ~ age + bp, data = sample.data, ## x = TRUE, y = TRUE, surv = TRUE) ## ## Model Tests Discrimination ## Indexes ## Obs 200 LR chi2 0.15 R2 0.001 ## Events 131 d.f. 2 Dxy 0.015 ## Center -0.5451 Pr(> chi2) 0.9277 g 0.038 ## Score chi2 0.15 gr 1.039

## Pr(> chi2) 0.9279 ## ## Coef S.E. Wald Z Pr(>|Z|) ## age -0.0051 0.0176 -0.29 0.7735 ## bp -0.0024 0.0101 -0.24 0.8088 计算校正c-index与未校正c-index

v <- validate(fit.cph, dxy=TRUE, B=1000) Dxy = v[rownames(v)==\"Dxy\", colnames(v)==\"index.corrected\"] orig_Dxy = v[rownames(v)==\"Dxy\", colnames(v)==\"index.orig\"] bias_corrected_c_index <- abs(Dxy)/2+0.5 #计算校正c-index orig_c_index <- abs(orig_Dxy)/2+0.5 #计算未校正c-index bias_corrected_c_index ## [1] 0.521968

orig_c_index ## [1] 0.5074744 4. 总结与讨论

C-index是Cox回归模型评价中最为重要的参数，它反映了模型预测效果的优劣，但这个参数在SPSS中无法计算，我们今天介绍的用R语言计算的两种方法，能掌握其中一种即可，推荐第一种，因其同时报告了C-index标准误，可以很方便的计算出C-index的可信区间。