文章目录
  1. 1. 前言
  2. 2. 方法1
  3. 3. 方法2
  4. 4. 方法3
  5. 5. 方法4
  6. 6. 方法5

前言

声明: 所有计算基于R软件,如果有人问其他软件如何实现,请自行Google

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

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

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

  • Harrell’s C
  • C-statistic by Begg et al. (survAUC::BeggC)
  • C-statistic by Uno et al. (survC1::Inf.Cval; survAUC::UnoC)
  • Gonen and Heller Concordance Index for Cox models (survAUC::GHCI, CPE::phcpe, clinfun::coxphCPE)

    不同C指数的详细含义可以看这里

我们这里主要探讨Harrell C,因为文献中使用最广泛。

说点题外话,Frank E. Harrell是一位著名生物统计学家,写了很多包,其中Hmisc,H是指他的姓的首字母 ,misc指的是miscellaneous,里面有很多五花八门的函数; rms就是和他的书Regression Modeling Strategies配套的R包。

方法1

利用Hmisc包中的rcorr.cens函数
局限:

  • 只能处理一个预测变量
  • 对超过2分类的分类变量处理粗糙
# 加载包及生成数据框,这里生成数据框主要是为了方便大家理解,因为大家通常都是将Excel的数据读进R,存储为数据框格式
library(survival)
library(Hmisc)
age <- rnorm(200, 50, 10)
bp  <- rnorm(200,120, 15)
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)
#让我们看一下生成的例子数据的前6
head(sample.data)

##        age       bp       os death
## 1 33.18822 114.6965 1.106501 FALSE
## 2 41.86970 123.2265 1.365944 FALSE
## 3 50.41484 124.9522 0.867119 FALSE
## 4 45.66936 127.3237 1.155765  TRUE
## 5 39.79024 134.8846 1.257501  TRUE
## 6 31.89088 140.9382 1.125504 FALSE    

rcorr.cens的代码及结果,第一个值就是C指数,同时也有Dxy的值

rcorr.cens(sample.data$age, Surv(sample.data$os, sample.data$death))

##        C Index            Dxy           S.D.              n        missing 
##   4.528492e-01  -9.430156e-02   5.565299e-02   2.000000e+02   0.000000e+00 
##     uncensored Relevant Pairs     Concordant      Uncertain 
##   1.290000e+02   3.172800e+04   1.436800e+04   8.072000e+03

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

r <- rcorrcens(Surv(os, death) ~ age + bp,data = sample.data)
r
## Somers' Rank Correlation for Censored Data    Response variable:Surv(os, death)
## 
##         C    Dxy  aDxy    SD    Z      P   n
## age 0.453 -0.094 0.094 0.056 1.69 0.0902 200
## bp  0.498 -0.003 0.003 0.054 0.06 0.9517 200

方法2

直接从survival包的函数coxph结果中输出,需要R的版本高于2.15.

library(survival)
sum.surv <- summary(coxph(Surv(os, death) ~ age + bp,data = sample.data))
c_index <-sum.surv$concordance
c_index

## concordance.concordant            se.std(c-d) 
##             0.54469239             0.02788881

可以看出这种方法输出了C指数,也输出了标准误,那么95%可信区间就可以通过加减1.96*se得到。并且这种

方法也适用于很多指标联合。

方法3

利用函数survConcordance,这种方法和方法2类似,输出的结果相同

fit <- coxph(Surv(os, death) ~ age + bp,data = sample.data)
survConcordance(Surv(os, death) ~ predict(fit, sample.data))
## Call:
## survConcordance(formula = Surv(os, death) ~ predict(fit, sample.data))
## 
##   n= 200 
## Concordance= 0.5446924 se= 0.02788881
## concordant discordant  tied.risk  tied.time   std(c-d) 
##  8641.0000  7223.0000     0.0000     0.0000   884.8563

方法4

利用survcomp包,安装这个包我就不在这里赘述了。

library(survcomp)
fit <- coxph(Surv(os, death) ~ age + bp, data = sample.data)
cindex <- concordance.index(predict(fit),surv.time = sample.data$os, surv.event = sample.data$death,method = "noether")
cindex$c.index; cindex$lower; cindex$upper

这种方法的优点就是可以直接输出95%可信区间,不需要自己再进行计算。说实话语法有点繁琐,感觉不爽!

方法5

利用rms包中的cph函数和validate函数,可提供un-adjusted和bias adjusted C指数两种,
未校正的C指数的结果和方法4是相同的。

library(rms)
#这里设置种子,目的是为了能重复最后的结果,因为validate函数的校正结果是随机的。但是我也发现即使设置了随机数种子,这个矫正的结果也不停在变,目前还没有找到解决办法,希望知道的大侠能给与指导。
set.seed(1)
fit.cph <- cph(Surv(os, death)~ age + bp, data = sample.data, x = TRUE, y = TRUE, surv = TRUE)

# Get the Dxy
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"]

# 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.5325809
 orig_c_index
## [1] 0.5446924

这种方法我觉得最大的优势就是给出了校正的C指数,但是都没有95%可信区间。并且最大的缺点就是代码比较多。为了简化,我自己写了一个函数:

cindex.boot <- function(fit.cph) {
set.seed(1234)
validate <- rms::validate(fit.cph, dxy = TRUE, B = 1000)
cindex <- (validate["Dxy", c("index.orig","training","test","optimism","index.corrected")])/2 + 0.5
n <- validate["Dxy", c("n")]
res <- rbind(validate, C_index = c(cindex, n))
res["C_index","optimism"] <- res["C_index","optimism"] - 0.5
res
}

代码简化为:

fit.cph <- cph(Surv(os, death)~ age + bp, data= sample.data, x = TRUE, y = TRUE, surv = TRUE)
#结果请看最后一行
cindex.boot(fit.cph)
##           index.orig     training         test     optimism
## Dxy      0.089384771  0.103815344 0.0805584972  0.023256847
## R2       0.026114094  0.034118373 0.0214282183  0.012690154
## Slope    1.000000000  1.000000000 0.9201702953  0.079829705
## D        0.003561857  0.005008749 0.0027677656  0.002240983
## U       -0.001664800 -0.001667190 0.0008882841 -0.002555474
## Q        0.005226657  0.006675938 0.0018794815  0.004796457
## g        0.227956865  0.248434659 0.2036042516  0.044830408
## C_index  0.544692385  0.551907672 0.5402792486  0.011628423
##         index.corrected    n
## Dxy        0.0661279238 1000
## R2         0.0134239399 1000
## Slope      0.9201702953 1000
## D          0.0013208739 1000
## U          0.0008906738 1000
## Q          0.0004302001 1000
## g          0.1831264574 1000
## C_index    0.5330639619 1000

细心的读者可以看出,方法2、3、4、5的结果都是相同的。但不代表他们之间没有差别。本质上,方法2、3是相同的;4、5是相同的。这两类的区别就在于处理tied risk上,也就是当两个观测拥有相同的生存时间和相同的自变量X时,方法4和5忽略tied risk,而方法2和3则考虑了tied risk。

方法4和5的计算:Concordance = #all concordant pairs/#total pairs ignoring ties.

方法2和3的计算:Concordance = (#all concordant pairs + #tied pairs/2)/(#total pairs including ties)

说了那么多方法,唯一不同是否在计算时考虑tied risk,其他只是实现方法和函数不同罢了。那么我们能不能不要这么复杂,只需要二个函数来解决C指数和可信区间的事呢?当然!!

我写了如下函数,随心所欲!

cindex <- function(time = sample.data$os, event =sample.data$death, variable = c("age","bp"), data = sample.data, ties=TRUE,adj = FALSE){
    require(rms)
    surv <- Surv(time,event)
    form <- as.formula(paste("surv~",paste(variable,collapse=" + ")))
    fit.coxph <- coxph(form,data)
    fit.cph <- cph(form, data = data, x = TRUE, y = TRUE, surv = TRUE)

    if (ties==FALSE){
      require(survcomp)
      coxPredict <- predict(fit.coxph, data = data, type="risk")
      c_index <-concordance.index(x=coxPredict, surv.time=time, surv.event=event, method="noether")
      res <- paste(c_index$c.index, " (", c_index$lower, " - ", c_index$upper,")", sep = "")
  }

  else if (ties==TRUE) {
    sum.surv <- summary(fit.coxph)
    c_index <- sum.surv$concordance
    res <- paste(c_index[1], " (", c_index[1]-1.96*c_index[2], " - ", c_index[1]+1.96*c_index[2],")", sep = "")
  }
  if(adj == FALSE){
      bias_corrected_c_index <- NA
  }
  else if (adj==TRUE) { 
    set.seed(1234)
      v <- rms::validate(fit.cph, dxy = TRUE, B = 1000)
      Dxy  <-  v[rownames(v)=="Dxy", colnames(v)=="index.corrected"]
      bias_corrected_c_index  <- abs(Dxy)/2+0.5
      bias_corrected_c_index
  }
    final <- list()
    final["C-index and 95%CI"] <- res
    final["Bias corrected C-index"] <- bias_corrected_c_index
    final
}

最后的最后,我们用自编的函数来求解试试:

默认计算节点的情况

cindex(sample.data$os,event = sample.data$death,variable=c("age","bp"),data = sample.data,adj = TRUE)
## $`C-index and 95%CI`
## [1] "0.544692385274836 (0.490030309427598 - 0.599354461122074)"
## 
## $`Bias corrected C-index`
## [1] 0.533064

忽略节点的情况

cindex(sample.data$os,event = sample.data$death,variable=c("age","bp"), ties=FALSE,data = sample.data,adj = TRUE)

## $`C-index and 95%CI`
## [1] "0.544692385274836 (0.492539122789943 - 0.596845647759729)"
## 
## $`Bias corrected C-index`
## [1] 0.533064
文章目录
  1. 1. 前言
  2. 2. 方法1
  3. 3. 方法2
  4. 4. 方法3
  5. 5. 方法4
  6. 6. 方法5