1. 程式人生 > >【腫瘤預測模型系列】如何在R軟體中求一致性指數( Harrell'concordance index:C-index)?

【腫瘤預測模型系列】如何在R軟體中求一致性指數( Harrell'concordance index:C-index)?

今晚收到一封醫生好友的郵件,諮詢如何求Harrell的C-index?曾在丁香園論壇裡遇到過朋友求助,自己也嘗試回答過,論壇裡回答的言簡意賅,針對問題論問題,沒有詳細的原理說明,今天趁回覆朋友郵件的機會,就大致總結下自己對C-index的理解以及在R軟體中的計算過程。

所謂C-index,英文名全稱concordance index,中文裡有人翻譯成一致性指數,最早是由範德堡大學(Vanderbilt University)生物統計教教授Frank E Harrell Jr 1996年提出,主要用於計算生存分析中的COX模型預測值與真實之間的區分度(discrimination);現階段用的最多的是腫瘤患者預後模型的預測精度。一般評價模型的好壞主要有兩個方面,一是模型的擬合優度(Goodness of Fit),常見的評價指標主要有R方,-2log

L,AIC,BIC等等;另外一個是模型的預測精度,主要就是模型的真實值與預測值之間的差的大小,均方誤差,相對誤差等。從臨床應用的角度來說,我們更注重後者,即統計建模主要是用於預測,而從C-index的概念大家看出它屬於模型評價指標的後者,這一指標比前面提到的幾個指標看起來更高大上,一般文獻中用的也比較多。換句話說,如果預後模型建好,效果不錯,即使不知道如何計算C-index值,報告軟體輸出結果中的預測誤差是相同效果,再新增擬合優度會更能說明效果,這樣反而更實用。

C-index本質上是估計了預測結果與實際觀察到的結果相一致的概率,即資料所有病人對子中預測結果與實際結果一致的對子所佔的比例。有點類似於ROC曲線下面積。

C-index的計算方法是:把所研究的資料中的所有研究物件隨機地兩兩組成對子。以生存分析為例,對於一對病人,如果生存時間較長的一位,其預測生存時間長於生存時間較短的一位,或預測的生存概率高的一位的生存時間長於生存概率低的另一位,則稱之為預測結果與實際結果一致。

C-index的計算步驟為:

(1)產生所有的病例配對。若有n個觀察個體,則所有的對子數應為Cn2(組合數)?

(2)排除下面兩種對子:對子中具有較小觀察時間的個體沒有達到觀察終點及對子中兩個個體都沒達到觀察終點。剩餘的為有用對子。

(3)計算有用對子中,預測結果和實際相一致的對子數,即具有較壞預測結果個體的實際觀察時間較短。

(4)計算。C=一致對子數/有用對子數。

由上述計算方法可以看出,C-index在0.5-1之間。0.5為完全不一致,說明該模型沒有預測作用,1為完全一致,說明該模型預測結果與實際完全一致。在實際應用中,很難找到完全一致的預測模型,既往研究認為,C-index在0.50-0.70為較低準確度:在0.71-0.90之間為中等準確度;而高於0.90則為高準確度。

當C-index檢驗由同一樣本建成的模型時易造成偏倚,因此再採用重抽樣技術(Bootstrap)可以幾乎無偏倚的檢驗預測模型的準確度。Bootstrap是非引數統計中一種重要的估計統計量方差進而進行區間估計的統計方法,是現代統計應用較為廣泛的一種統計方法。

Bootstrap方法核心思想和基本步驟如下:

(1)採用重抽樣技術從原始樣本中抽取一定數量的樣本,此過程允許重複抽樣。

(2)根據抽出的樣本計算給定的統計量T。

(3)重複上述N次(一般大於1000),得到N個統計量T。

(4)計算上述N個統計量T的樣木方差,得到統計量的方差。

Bootstarap方法只是對單一樣本且樣本量較小的資料,如果資料集很大可以按照不同的比例將資料集拆分,一部分用於建模一部分用於驗證。關於交叉驗證(Cross-validation),由於篇幅有限,留作下次探討。

R軟體實現:

C-index的R軟體計算實現有兩種實現方法,一種是用到Harrell本人的的R包Hmisc package;另一種是Le Kang, Weijie Chen 2014年12月18日釋出的R compareCPackage;

############################

#### Method 1.Hmisc code ###

############################

data <- read.csv("survivaldta.csv") ###讀入csv格式資料####

library(Hmisc) ###載入Hmisc包,前提是安裝了,如果沒有安裝,請百度如果安裝R package,這裡就不詳細講了!

library(survival) ###載入survival包,主要用於建立模型###

f <- cph(Surv(time,death)~x1+x2+x3,data=survivldata) ###擬合cox模型

fp <- predict(f)###模型的預測值

cindex.orig=1-rcorr.cens(fp,Surv(time,death)) [[1]]###計算出的C-index

 

###############################

#### Method 2.compareC code ###

###############################

data <- read.csv("survivaldta.csv") ###讀入csv格式資料####

library(compareC) ###載入compareC包,前提是安裝了,如果沒有安裝,請百度如果安裝R package,這裡就不詳細講了!

library(survival) ###載入survival包,主要用於建立模型###

cindex <- cindex(Surv(time,death)  ~ x1+x2+x3,data=survivldata)###計算出的C-index

###############################

#### Bootstrap code ###

###############################

bootit=200
for(i in 1:bootit){
case=noNA[group=="long",]  
control=noNA[group=="<24",]
bootindex.case=sample(1:nrow(case),replace=T)
boot.case.data=case[bootindex.case,]
bootindex.control=sample(1:nrow(control),replace=T)
boot.control.data=control[bootindex.control,]
boot.data=rbind(boot.case.data,boot.control.data)
dstr.boot=svydesign(id=~1, prob=~inv_weight, fpc=~ssize, data=boot.data)
boot.fit=svycoxph(Surv(survival,surv_cens) ~x1+x2+x3,data=boot.data,x=TRUE,design=dstr.boot)
cindex.train=1-rcorr.cens(lp.boot,Surv(boot.data$survival, boot.data$surv_cens))[[1]]
cindex.test=1-rcorr.cens(lp_=.test,Surv(noNA$survival,noNA$surv_cens))[[1]]
bias=rep(1,bootit)
bias[i]=abs(cindex.train-cindex.test) }

以上R程式碼給的略微簡潔,還請見諒!

 

參考文獻:

Harrell FE, Califf RM, Pryor DB, Lee KL, and Rosati RA. (1982) Evaluating the yield of medical tests. The Journal of the American Medical Association, 247(18), 2543–2546

Pencina MJ and D’Agostino RB. (2004) Overall C as a measure of discrimination in survival analysis: model specific population value and confidence interval estimation. Statistics in Medicine, 23(13), 2109–2123

Kang L, Chen W, Petrick NA, and Gallas BD. (2014) Comparing two correlated C indices with right-censored survival outcome: a one-shot nonparametric approach. Statistics in Medicine, 34(4), 685–703, doi: 10.1002/sim.6370

Hmisc Reference manual:http://cran.r-project.org/web/packages/Hmisc/Hmisc.pdf

compareC Reference manual: http://cran.r-project.org/web/packages/compareC/compareC.pdf

Frank.Harrell :http://biostat.mc.vanderbilt.edu/wiki/Main/FrankHarrell



 

本公眾號精彩歷史文章:

04:如何在R軟體中求一致性指數( Harrell'concordance index:C-index)?

05:Nomogram 繪製原理及R&SAS實現.

06  : Lasso方法簡要介紹及其在迴歸分析中的應用

07  : 最優模型選擇中的交叉驗證(Cross validation)方法

08  : 用R語言進行分位數迴歸(Quantile Regression)

09  : 樣本資料中異常值(Outliers)檢測方法及SPSS & R實現

10  : 原始資料中幾類缺失值(Missing Data)的SPSS及R處理方法

11  :  [Survival analysis] Kaplan-Meier法之SPSS實現

12  :  [Survival analysis] COX比例風險迴歸模型在SPSS中的實現

13  :  用R繪製地圖:以疾病流行趨勢為例

14  :  資料探勘方法:聚類分析簡要介紹 及SPSS&R實現

15  :  醫學研究中的Logistic迴歸分析及R實現

16  :  常用的非引數檢驗(Nonparametric Tests)總結

17  :  高中生都能看懂的最小二乘法原理

18  :  R語言中可實現的常用統計假設檢驗總結(側重時間序列)

19  :  如何根據樣本例數、均數、標準差進行T-Test和ANOVA

20  :  統計學中自由度的理解和應用

21  :  ROC和AUC介紹以及如何計算AUC

22  :  支援向量機SVM介紹及R實現

23  :  SPSS如何做主成分分析?

24  : Bootstrap再抽樣方法簡介

25  :  定量測量結果的一致性評價及 Bland-Altman 法的應用 

26  :  使用R繪製熱圖及網路圖  

27  :  幾種常用的雙座標軸圖形繪製 

28  :  遺失的藝術—諾謨圖(Nomogram) 

29  :  Nomogram 繪製原理及R&SAS實現(二) 

30  :  WOE:信用評分卡模型中的變數離散化方法 

31  :  結構方程模型(SEM)簡介及教程下載  

32  :  重複測量的多因素方差分析SPSS實現操作過程 

回覆文章前程式碼數字如“04”即可檢視或直接檢視歷史文章。

公眾號:survival-analysis QQ:8243033

郵箱:8243033 @ qq.com 歡迎關注!