生存模型评价常用指标总结
生存分析：研究生存现象和响应时间数据及其统计规律的一种分析方法。
常见的生存分析方法有：
寿命表法：主要研究不同时间段内生存率的变化
单因素 Kaplan-Meier 分析：主要研究单个因素对生存率的影响
多因素的Cox分析：能支持多个影响因素进行分析
生存分析常见术语：
删失（Censoring）：删失是指在对样本进行观测或调查时，某个个体的确切生存时间未知，只知道其生存时间位于某个区间内。删失可分为右删失，左删失和区间删失。
右删失：（right censoring) 当删失时间C＞生存时间X时，绝大多数删失为右删失。
左删失：（left censoring)当个体的生存时间X＜删失时间C，称为左删失。
区间删失：（interval censoring） 个体的生存时间只能确定在某一区间之间，但不知道确切时间 （L，R）L代表左观测时间，R代表右观测时间。
截尾（Truncation)：截尾数据指只有个体的生存时间处于某区特定区间才能被观测到，在区间之外的个体不能被观测到。
打结(tied) :当不同个体拥有相同的失效时间，称之为打结。例如，i,j为两个不同的个体（i≠j）当Ti=Tj时，发生打结。
一个生存模型的建立需要经过描述，模型拟合，和预测几个阶段。当我们根据病人的生存数据拟合构建了生存模型后，需要评估模型的预测能力，考察它将研究对象正确的划分为病人和非病人的能力。模型预测能力是评价模型的优劣性的重要指标。
目前，用于生存模型风险预测效果评价的指标有以下几种：
一、Harrell’s c 指标、Somers’D 指标
二、ROC曲线法、时依ROC曲线法
三、IDI、NRI法
一、Harrell’s c 指标、Somers’D 指标
Harrell’s c 指标
**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 个对子）。
根据对子的数据特征，可进一步的区分不同类型的对子
**可用对（usable pair）：
**两个实际生存时间不一致的被试（Yi≠Yj）组成的对子称为可用对。（pencina 等人对可用对的范围进行了更严格的划分，两个病人的实际生存时间与预测的生存概率都不一致时（Yi≠Yj、Xi≠Xj），才被选为可用对。[1] [2]）
一致对（concordant pair）：在后文我们用P表示，如果一个配对中两个被试实际生存时间的趋势与模型预测的生存概率趋势一致（即当Yi>Yj时，Xi>Xj, 或 Yi＜Yj时，Xi＜Xj），这个配对为一致对。
非一致对 （discordant pair）:在后文我们用Q表示。如果一个配对中两个被试实际生存时间的趋势与模型预测的生存概率趋势不一致（即当Yi>Yj时，Xi＜Xj，或 Yi＜Yj时，Xi＞Xj），这个配对为不一致对。
不确定对（unsure pair） ：两个被试实际生存时间不同，而模型预测的生存概率相同时，该配对为不确定对。Yi=Yj, Xi≠Xj的配对，在后文我们用T表示。（在pencina 等人后续的定义中，unsure pair ，也属于不可用对[2]。）
公式：C=一致对对数数/可用对总对数
在 Harrell 等人的定义中，可用对总对数=P+Q+T，实际生存时间（Yi≠Yj）不一致的配对即为可用对[3]，而在pencina等人后续的定义中，可用对总对数=P+Q，不确定对（Xi≠Xj, Yi=Yj的配对）为不可用对[2]。
Somers’D 指标
其概念与Harrell'c 有相似之处，由 1962年 Somers 提出[4]。
公式：Somers'D(X|Y) =（ 一致对对数-不一致对对数）/可用对总数
有人认为，Somers’D 指标与Harrell's C指标之间存在关系，C=(D+1)/2，也因此认为，当D=(P-Q)/(P+Q+T)时， C=(P+0.5T)/(P+Q+T)
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
# 可以看出，两种不同的函数代码，计算出的结果是一样的（可对比变量age计算的C 及somers'D指标，完全一致哦）
参考文献
[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)，认为该病例阴性。
灵敏度（Sensitivity）: P(pi＜c | Yi=1) 最终生存结局为存活的病例中，这些病例的生存概率函数pi＞c的比例
特异度（Specificity）: P(pi≥c | Yi = 0) 最终生存结局为死亡的病例中，这些病例的生存概率函数pi≤c的比例
令c∈（-∞，+∞），获得在不同的c值下，sensitivity 和对应的specificity的值，以sensitivity 为横坐标，1-specificity 为纵坐标，绘制生存模型的ROC曲线。通过ROC曲线的AUC值，可评估生存模型的准确性。
可用survivalROC包分析生存模型的ROC曲线
代码：
library(survivalROC)
data(mayo)
nobs <- NROW(mayo)
cutoff <- 365
#使用survival包的survivalROC函数构建生存模型
Mayo4.2= survivalROC(Stime=mayo$time,
status=mayo$censor,
marker = mayo$mayoscore4,
predict.time = cutoff, method="KM")
#绘图
plot(Mayo4.2$FP, Mayo4.2$TP, type="l", xlim=c(0,1), ylim=c(0,1),
xlab=paste( "FP", "\n", "AUC = ",round(Mayo4.2$AUC,3)),
ylab="TP",main="Mayoscore 4, Method = KM \n Year = 1")
abline(0,1)
代码输出图片如下图：
时依ROC
时间依赖性受试者工作特征曲线（time dependent receiver operating characer curve, ROC），简称时依ROC曲线。经典的ROC分析中，研究对象的疾病状态是固定不变的，然而，很多疾病对患者生存的影响是时间依赖性的。因此，在ROC曲线的基础上进一步发展出了时依ROC曲线[1,2]。传统的ROC曲线，以生存结局Y=0，Y=1的二元变量来作为区分病例阴阳性的标准，计算模型在不同的阈值C的情况下的灵敏度、特异度并绘制ROC曲线。最后根据ROC曲线的AUC值对模型的预测能力进行评价。而时依ROC曲线，在传统的ROC曲线基础上加入时间节点t，当病人的生存时间T＜t判定该病例为阳性，T＞t判定该病例为阴性。这样，设置不同的时间节点t，可获得一系列ROC曲线[1,2]。
目前，有三大类时依ROC模型，它们的主要区别在于灵敏度和特异度的具体计算方法上，我们将一一介绍。
我们用Mi来表示模型预测的病人生存概率，Ti来表示病人的实际生存时间，c为阈值。
累积/动态(Cumulative/Dynamic)ROC曲线
Sensitivity(c,t):P (Mi>c | Ti ≤ t) = P (Mi ＞c | Ti ≤ t)
Specificity(c,t):P(Mi≤c | Ti ＞t) = P(Mi ≤ c | Ti ＞t)
灵敏度：生存时间小于t的人群之中，Mi大于阈值c的人群所占的比例 特异度：生存时间大于t的人群之中，Mi小于等于阈值c的人群所占的比例
时点/静态(Incident/Static)ROC曲线
Sensitivity(c,t):P (Mi>c | Ti = t) = P (Mi ＞c | Ti = t)
Specificity(c,t):P(Mi≤c | Ti ＞t*)= P(Mi ≤ c | Ti ＞t)
灵敏度：在生存时间等于t的人群中，Mi大于阈值c的人群所占的比例
特异度：生存时间大于t的人群中，Mi小于等于阈值c的人群所占的比例
时点/动态(Incident/dynamic )ROC曲线 Sensitivity(c,t):P (Mi>c | Ti ＝ t) = P (Mi ＞c | Ti＝ t) Specificity(c,t):P(Mi≤c | Ti ＞t) = P(Mi ≤ c | Ti ＞t)
灵敏度：在生存时间等于t的人群中，Mi大于阈值c的人群所占的比例
特异度：生存时间大于t的人群中，Mi小于等于阈值c的人群所占的比例
可将多个时间点t下绘制的时依ROC曲线绘制在同一张图来表示，如下图：
图片来源：[1]
时依ROC曲线，可用risksetROC包进行分析
代码：
library(risksetROC)
library(survival)
library(MASS)
data(VA) #加载数据集 VA
#> head(VA)
# stime status treat age Karn diag.time cell prior
#1 72 1 1 69 60 7 1 0
#2 411 1 1 64 70 5 1 10
#3 228 1 1 38 60 3 1 0
#4 126 1 1 63 60 9 1 10
#5 118 1 1 65 70 11 1 10
#6 10 1 1 49 20 5 1 0
survival.time=VA$stime #实际生存时间
survival.status=VA$status #事件发生=1，无事件发生=0
score <- VA$Karn
cell.type <- factor(VA$cell)
tx <- as.integer( VA$treat==1 )
age <- VA$age
survival.status[survival.time>500 ] <- 0
survival.time[survival.time>500 ] <- 500
fit0 <- coxph( Surv(survival.time,survival.status)
~ score + cell.type + tx + age, na.action=na.omit )
eta <- fit0$linear.predictor
#使用risksetROC函数， 建立incident/dynamic型时依ROC曲线
ROC.CC30=risksetROC(Stime=survival.time, status=survival.status,
marker=eta, predict.time=30, method="Cox",
main="ROC Curve", lty=2, col="red")
代码输出图片如下图：
rickset函数包进行时依ROC曲线分析，也可参考该网站的代码
参考文献：[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
三、IDI、NRI法
IDI，综合区分改善度（Integrated Discrimination Improvement, IDI）
公式：IDI =（ISnew – ISold）- （IPnew – IPold）
IS：不同分类阈值设定下的总灵敏度（Integral of sensitivity over all possible cutoff values）
IP：不同分类阈值设定下的1-特异度
其中IPnew、IPold表示在患者组中，新模型和旧模型对于每个个体预测疾病发生概率的平均值，两者相减表示预测概率提高的变化量，对于患者来说，预测患病的概率越高，模型越准确，因此差值越大则提示新模型越好。而IPnew,non-events、IPold,non-events表示在非患者组中，新模型和旧模型对于每个个体预测疾病发生概率的平均值，两者相减表示预测概率减少的量，对于非患者来说，预测患病的概率越低，模型越准确，因此差值越小则提示新模型越好。**总体来说IDI越大，则提示新模型预测能力越好。若IDI>0，则为正改善，说明新模型比旧模型的预测能力有所改善，若IDI<0，则认为相比旧模型而言，新模型预测能力下降，若IDI=0，则认为新模型没有改善。**
NRI，净重分类改善度（Net Reclassification Improvement, NRI）
公式：NRI =（ISnew- ISold）+（IPnew - IPold）=（ISnew + IPnew）-（ISold + IPold）
IS：不同分类阈值设定下的总灵敏度（Integral of sensitivity over all possible cutoff values）
IP：不同分类阈值设定下的1-特异度
若NRI>0，则为正改善，说明新模型比旧模型的预测能力有所改善；若NRI<0，则为负改善，新模型预测能力下降；若NRI=0，则认为新模型没有改善。[2]
可用survIDINRI包计算模型的IDI、NRI值
代码：
library('survIDINRI')
library('survival')
#--- sample data (pbc in survival package) ---
D=subset(pbc,select=c("time","status","age","albumin","edema","protime","bili"))
###让我们来看一下建立的这个D数据集包含的内容
#> D
# time status age albumin edema protime bili
#1 400 2 58.76523 2.60 1.0 12.2 14.5
#2 4500 0 56.44627 4.14 0.0 10.6 1.1
#3 1012 2 70.07255 3.48 0.5 12.0 1.4
#4 1925 2 54.74059 2.54 0.5 10.3 1.8
#5 1504 1 38.10541 3.53 0.0 10.9 3.4
#6 2503 2 66.25873 3.98 0.0 11.0 0.8
D$status=as.numeric(D$status==2) # 先根据status的数值是否等于2，将数据转换成布尔值TRUE或者FALSE，再利用as.numeric将TRUE转换成1，FALSE转换成0。
D=D[!is.na(apply(D,1,mean)),] ; dim(D) #！is.na去空值，apply(D,1,mean)对行取平均值，
mydata=D[1:100,] #选第1~100行
t0=365*5
indata1=mydata;
indata0=mydata[,-7] ; n=nrow(D) ; #去掉倒数第7列
covs1<-as.matrix(indata1[,c(-1,-2)]) #去掉倒数第一行，去掉倒数第二列
covs0<-as.matrix(indata0[,c(-1,-2)])
#让我们来看一下最后建立的模型,
##################################################
#> covs1
# age albumin edema protime bili
#1 58.76523 2.60 1.0 12.2 14.5
#2 56.44627 4.14 0.0 10.6 1.1
#3 70.07255 3.48 0.5 12.0 1.4
#4 54.74059 2.54 0.5 10.3 1.8
#5 38.10541 3.53 0.0 10.9 3.4
#6 66.25873 3.98 0.0 11.0 0.8
#> covs0
# age albumin edema protime
#1 58.76523 2.60 1.0 12.2
#2 56.44627 4.14 0.0 10.6
#3 70.07255 3.48 0.5 12.0
#4 54.74059 2.54 0.5 10.3
#5 38.10541 3.53 0.0 10.9
#6 66.25873 3.98 0.0 11.0
#可以看到，模型covs1有5个指标，模型covs0有4个指标
#--- inference ---
x<-IDI.INF(mydata[,1:2], covs0, covs1, t0, npert=200) ;
#详解：第一个参数输入生存时间，和生存结局组成的data.frame，第二，第三个模型指标两个模型指标，t0时间，第三个参数，时间节点，第五个参数，置换检验的次数
#--- results ---
IDI.INF.OUT(x) ; #可以看计算结果
#> IDI.INF.OUT(x)
# Est. Lower Upper p-value
#M1 0.073 0.020 0.137 0.01
#M2 0.444 0.201 0.677 0.00
#M3 0.042 0.010 0.118 0.00
#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.