C-index/C-statistic 計算的5種不同方法及比較
阿新 • • 發佈：2019-01-01
前言
宣告： 所有計算基於R軟體，如果有人問其他軟體如何實現，請自行Google。
評價一個預測模型的表現可以從三方面來度量：
- 區分能力（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指數:
- 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