**Harrell’s c 指标有时也称为C-index,**它反映的是模型预测结果与实际情况的一致程度,即为:在一次生存模型分析中随机挑选的两个病人,模型所预测的生存时间更长的那个病人,他的实际生存时间也更长的概率。由于估算出预测的生存时间较复杂,Harrell等人指出,实际运用中,模型的预测生存时间与预测生存概率函数成对应关系,计算Harrell's c 指标时可以用模型的预测生存概率函数来代替预测生存时间[1]。
详解:假设一次生存分析中共有n个病人,病人的实际生存时间用Yi来表示(Y1,Y2, ...Yn),利用病人的生存数据拟合生存模型后,模型预测的生存时间用Ti来表示(T1,T2, …Tn),预测的生存概率函数用Xi来表示(X1, X2, …Xn)。接下来,计算Harrell's c 指标,对模型的预测能力进行评估。原理:将n个病人之间两两随机配对,对于一个配对,共获得四个观测值(Xi,Yi,Xj,Yj),i≠j。
对子(pair):在生存分析中,随机挑选两个被试,i、j,这两个被试的实际生存时间Y与建立的生存模型预测的生存概率X组成一个对子(Xi, Yi, Xj, Yj)。在一次共有n个被试的生存分析中,所有被试一共可以组成 nx(n-1) 个对子(若考虑 i, j 和 j, i 的重复情况,也可认为共组成 nx(n-1)/2 个对子)。
Hmics包的rcorrcens函数可用于计算模型的Harrell's c 指标和Somer's D指标
代码:
library(Hmisc::rcorrcens) #载入Hmisc包的rcorrcens函数 library(survival) #先使用 survival 包构建生存模型,下面建立一个虚拟模型 age <- rnorm(400, 50, 10) #构建样本量为400,均值为50,方差为10的正态分布样本。rnorm(n, mean = 0, sd = 1) bp <- rnorm(400,120, 15) bp[1] <- NA #构建样本量为400,均值为120,方差为15的正态分布样本。rnorm(n, mean = 0, sd = 1) d.time <- rexp(400) #构建指数分布函数 cens <- runif(400,.5,2) death <- d.time <= cens #很有意思的写法,两个向量,d.time和cens之间的元素一对一比较大小,将大小比较的结果转换为布尔值并赋值给death向量 d.time <- pmin(d.time, cens) #取d.time,cens中的向量的较小值赋值给d.time向量
#Hmisc包中的rcorr.cens、rcorrcens都可用来计算Harrell's、Somer's D指标,只是代码的写法有些区别,得到的结果一致 #写法一: rcorr.cens(age, Surv(d.time, death)) #年龄对生存结果的影响,建立模型并计算Harrell's、Somer's D指标 ####rcorr.cens handles one predictor variable. rcorrcens computes rank correlation measures separately by a series of predictors.
#> head(Surv(d.time, death)) #[1] 0.1081747 0.8427812+ 0.9699464+ 1.0169367+ 0.1045871 #[6] 0.3790691 #rcorr.cens函数计算模型的Harrell's、Somer's D指标 #> rcorr.cens(age, Surv(d.time, death)) # C Index Dxy S.D. # 5.224652e-01 4.493045e-02 3.691319e-02 # n missing uncensored # 4.000000e+02 0.000000e+00 2.910000e+02 #Relevant Pairs Concordant Uncertain # 1.361660e+05 7.114200e+04 2.343400e+04 #> #写法二 r <- rcorrcens(Surv(d.time, death) ~ age + bp) #年龄、血压对生存结果的影响,建立模型并分别计算Harrell's、Somer's D指标 #结果 #> r #Somers' Rank Correlation for Censored Data Response variable:Surv(d.time, death)
# C Dxy aDxy SD Z P n #age 0.522 0.045 0.045 0.037 1.22 0.2235 400 #bp 0.477 -0.045 0.045 0.038 1.17 0.2405 399
[1] Michael J Pencina, Ralph B D'Agostino. Overall C as a measure of discrimination in survival analysis: model specic population value and condence interval estimation. Stat Med. 2004 Jul 15;23(13):2109-23.
[2] R Van Oirbeek, E Lesaffre. An application of Harrell’s C-index to PH frailty models. Stat Med. 2010 Dec 30;29(30):3160-71.
[3]Harrell,Frank,E. Evaluating the Yield of Medical Tests. The Journal of the American Medical Association.
[4]Robert H. Somers. A New Asymmetric Measure of Association for Ordinal Variables. American Sociological Review , Dec., 1962, Vol. 27, No. 6. pp. 799-811
二、ROC曲线法、时依ROC曲线法
ROC曲线法
受试者工作特征曲线( receiver operating characer curve, ROC)简称ROC曲线,通过设定不同的阈值(区分阳性、阴性的临界值),并计算在不同的阈值下,计算生存模型的真阳性率,假阳性率。并以灵敏度(又叫真阳性率, sensitivity)为纵坐标,假阳性率(false positive rate, FPR)=(1-特异度)绘制而来。ROC曲线的曲线下面积(area under the ROC curve, AUC)可反映诊断实验的准确度。一般认为,AUC 在0.5~0.7之间,表示诊断准确度较低,0.7~0.9之间,准确度中等,0.9以上,准确度较高。
ROC曲线法常见概念解析:
真阳性 true positive :判定结果为阳性,实际上也为阳性
假阴性 false negative:判定结果为阴性,实际上为阳性
真阴性 true negative:判定结果为阴性,实际上也为阴性
假阴性 false positive:判定结果为阳性,实际上为阴性
灵敏度 sensitivity = TP/TP+FN
特异度 specificity= TN/TN+FN
下面详细阐述在生存模型中,如何用ROC曲线来评估模型的准确性
首先利用病人的生存数据拟合构建生存模型
Y:病人的生存结局,我们可用 Y=1 表示存活, Y=0 表示死亡
p: 该模型下不同病人的生存概率函数
C: 分类阈值(区分病例阴阳性的临界值),当生存概率p大于c时(p > c) ,认为该病例阳性,小于c时(pi ≤ c),认为该病例阴性。
参考文献:[1]Patrick J. Heagerty Yingye Zheng. Survival Model Predictive Accuracy and ROC Curves. Biometrics 61, 92–105
[2]Heagerty, P.J., Lumley, T., Pepe, M. S. Time-dependent ROC Curves for Censored Survival Data and a Diagnostic Marker. Biometrics, 56, 337 – 344
#M1:IDI #M2:NRI #M3:中值改善度(median improvement) #第一行,估计值,第二行,第三行2.5% 97.5% 置信区间,最后一行,p值 # #--- Graphical presentaion of the estimates --- IDI.INF.GRAPH(x) ;#根据计算结果输出图片
代码输出图片:
参考文献:
[1]Kathleen F Kerr, Robyn L McClelland, Elizabeth R Brown, Thomas Lumley. Evaluating the Incremental Value of New Biomarkers With Integrated Discrimination Improvement. Am J Epidemiol. 2011 Aug 1;174(3):364-74.
[2]Michael J Pencina, Ralph B D'Agostino Sr, Ralph B D'Agostino Jr, Ramachandran S Vasan. Evaluating the added predictive ability of a new marker: From area under the ROC curve to reclassification and beyond. Stat Med. 2008 Jan 30;27(2):157-72; discussion 207-12.