题记：从今天开始，本公众号将陆续发布笔者本人的R语言学习笔记，从一个医生的视角，与您一起学习R语言~
一、背景介绍
近十年来，临床研究中有一类预测模型构建与验证类的文章数量逐渐增多。简言之，预测模型是通过已知参数来预测临床未知的结局，而模型本身就是一个数学公式。也就是把已知的参数通过这个所谓的模型就可以计算出未知的结局发生的可能性，即是预测。
临床预测模型的统计学本质就是回归建模分析，回归的本质就是发现因变量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. R2 是总体上（Overall）评估模型优劣的参数，事实上就是综合了区分能力和一致性的度量指标。模型决定系数更综合，但略显粗糙。
二、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指数两种。
代码及代码解读，结果解读如下：
> # 模拟一组数据并设置为数据框结构
> age <->
> bp <->
> d.time <->
> cens <->
> death <- d.time=""><=>
> os <- pmin(d.time,="">
> sample.data <- data.frame(age="age,bp" =="" bp,os="os,death" =="">
> head(sample.data) # 展示数据框sample.data的前6行
> # 方法1. {survival}包
> library(survival) # 载入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.53318557 0.02741619
> #方法2. {rms}包
> library(rms)
> set.seed(1) # 这里设置种子，目的是为了能重复最后的结果，因为validate函数的校正结果是随机的。
> dd<->
> options(datadist='dd')
> fit.cph <- cph(surv(os,="" death)~="" age="" +="" bp,="" data="sample.data," x="TRUE," y="TRUE," surv="">
> fit.cph # 模型参数 Dxy*0.5+0.5 即是c-index
> # Get the Dxy
> v <- validate(fit.cph,="" dxy="TRUE," b="">
> 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.5152632
> orig_c_index
[1] 0.5331856
联系客服
