【原】模型预测能力的比较：IDI[综合判别改善指数]

同净重分类改善指数（Net Reclassification Improvement Index，NRI）一样，综合判别改善指数（Integrated Discrimination Improvement Index，IDI）也用于两个模型预测能力的比较。不同的是NRI需要先根据设定的阈值将结局进行分类，然后比较两个模型在设定的阈值上正确分类的比例。而IDI直接比较两个模型的结局事件的预测概率，计算公式为：

IDI=(P_new|event-P_std|event)-(P_new|nonevent-P_std|nonevent)

Event是指发生结局事件的病例组人群，而Nonevent则是没有发生结局事件的对照组人群。

P_new|event和P_std|event表示发生结局事件的人群中，新模型和和标准模型（旧模型，基础模型）预测的每个个体发生结局事件的概率的均值。对患病人群来讲，两者的差值越大，模型的预测能力越好。

P_new|nonevent和P_std|nonevent表示没有发生结局事件的人群中，新模型和标准模型（旧模型，基础模型）预测的每个个体发生结局事件的概率的均值。对未患病人群来讲，两者的差值差值越小，模型的预测能力越好。

IDI就是两部分差值的差。IDI越大提示模型的预测能力越好。IDI>0，新模型的预测能力比旧模型有改善；IDI<0，新模型的预测能力有下降；IDI=0，新模型的预测能力没有改善。

两个logistic模型IDI的计算可借助PredictABEL包，两个Cox模型IDI的计算可使用survIDINRI包。

示例1：两logistic模型的IDI。示例数据同《两个logistic模型预测能力的比较：NRI[净重分类改善指数]》，该文中采用PredictABEL::reclassification计算NRI时已经同时输出了IDI结果。

【1】数据导入，并将数据随机分训练集和测试集

library(bruceR)

PE<- import("D:/Temp/PE.xlsx")

PE<-na.omit(PE)

#暴力删除缺失值

library(caret)

library(ggplot2);library(lattice)

set.seed(20240224)

split<-createDataPartition(PE$PE,p=0.7,list=FALSE)

split<-as.vector(split)

trainset<-PE[split,]

testset<-PE[-split,]

【2】利用训练集建立两个模型作为标准模型和新模型

在训练集上采用最优子集回归筛选到6个变量：age,BMI,ToS,CA153,CDU,stage。利用其中的人口学特征和血检指标建立标准模型PE.std：

PE.std<-glm(PE~age+BMI+ToS+CA153,family=binomial(link=logit),data=trainset,x=TRUE)

Pstd.train<- predict(PE.std,trainset,type ="response") #即fitted(PE.std)，获得训练集中每个观测值的预测概率(患病风险）

Pstd.test<- predict(PE.std,testset,type ="response") #测试集每个观测值的风险

在标准模型的基础上添加影像学和病理学特征建立新模型PE.new：

PE.new<-glm(PE~age+BMI+ToS+CA153+CDU+stage,family=binomial(link=logit),data=trainset,x=TRUE)

Pnew.train<-fitted(PE.new)

Pnew.test<-predict(PE.new,testset,type="response")

【3】两logistic模型的IDI计算

通过IDI的定义来计算一下并不难，以两个logistic模型在训练集IDI为例：

Event.train<-trainset[trainset$PE==1,]

Nonevent.train<-trainset[trainset$PE==0,]

P_new.event<- mean(predict(PE.new,Event.train,type="response"))

P_std.event<- mean(predict(PE.std,Event.train,type ="response"))

P_std.nonevent<- mean(predict(PE.std,Nonevent.train,type ="response"))

P_new.nonevent<- mean(predict(PE.new,Nonevent.train,type="response"))

IDI<-(P_new.event-P_std.event)+(P_std.nonevent-P_new.nonevent)

IDI

0.09049398

采用PredictABEL::reclassification可以获得更多的统计量：

library(PredictABEL)

reclassification(data=trainset,cOutcome=14,predrisk1=Pstd.train,predrisk2=Pnew.train,cutoff=c(0,0.2,0.5,1))

#data指定包含结局变量和预测变量的数据集；cOutcome结局变量所在的列是第几列；predrisk1指定标准模型的预测概率；predrisk2指定新模型的预测概率；cutoff指定风险类别的界值，第一个值必须是0，最后一个值必须是1，计算IDI用不到cutoff值，但由于该函数同时计算NRI，该参数需要指定 

测试集IDI：reclassification(data=testset,cOutcome=14,predrisk1=Pstd.test,predrisk2=Pnew.test,cutoff=c(0,0.2,0.5,1))

示例2：两Cox模型的IDI。示例数据同《绘制Cox回归模型的决策曲线》。

【1】数据导入，数据集拆分

library(survival)

GBSG<-gbsg

GBSG$grade<-factor(GBSG$grade)

library(caret)

library(ggplot2);library(lattice)

set.seed(20240224)

split<-createDataPartition(GBSG$status,p=0.7,list=FALSE)

split<-as.vector(split)

trainset<-GBSG[split,]

testset<-GBSG[-split,]

trainset$Subset<-0

testset$Subset<-1

【2】两个Cox回归模型的IDI计算

以age、meno、size、nodes、grade五个变量建立的模型为旧模型，利用全子集回归获得4个变量size、nodes、pgr、hormon建立的模型为新模型。

TE<-as.matrix(trainset[,c("rfstime","status")]) #生存时间和结局的数据集，此处也建议转换为矩阵格式。第一列是生存时间，第二列是结局状态（1=结局事件，0=删失）

Covs.std<-model.matrix(~age+meno+size+grade+nodes,data=trainset)[,-1] #基础模型的预测变量矩阵。grade是一个三分类变量，后面我们计算IDI的函数IDI.INF::survIDINRI不允许有因子变量，通过model.matrix()可将因子变量设置为哑变量并保持连续变量不变来创建模型矩阵，[,-1]是从 model.matrix 的输出中模型矩阵删除Intercept项。如果直接使用命令as.matrix(trainset[,c(2,3,4,5,6)]) 会将因子变量的赋值转换为数值格式（第2,3,4,5,6列分别是age、meno、size、grade、nodes）,后面运行程序时会报错

Covs.new<-as.matrix(trainset[,c(4,6,7,9)]) #新模型的预测变量矩阵。第4,6,7,9分别是size、nodes、pgr、hormon。并转换为矩阵格式

library(survIDINRI)

IDI.train.5Y<-IDI.INF(indata=TE,covs0=Covs.std,covs1=Covs.new,t0=365.25*5,npert=500,seed1=20240224)

#indata指定生存时间和结局的数据集，要求两列，第一列是生存时间，第二列是结局状态（1=结局事件，0=删失）；covs0和covs1分别制定基础模型和新模型的预测变量，数据结构需要是矩阵格式，不允许有因子、字符和缺失值，因子变量可以通过model.matrix()进行哑变量编码；t0指定评估时间点，模型会计算t0时的风险评分作为该时点的结局事件发生概率；npert指定perturbation-resampling的迭代次数；npert.rand指定扰动重抽样随机数；seed1指定生成扰动重抽样随机数的种子，因为置信区间的计算通过重抽样来完成，每次运行命令得到的置信区间会有差异，为保持结果的可重现性，可设置随机种子；npert.rand；alpha指定需要计算的置信区间(1-alpha/2)，默认95%CI

IDI.INF.OUT(IDI.train.5Y)

M1即IDI，包括点估计值为0.029、95%置信区间[-0.034,0.079]及P值=0.371。M2是连续NRI的点估计和区间估计值；M3是风险得分的中位值改善值。本例新模型对旧模型的预测能力有改善，但这种改善没有统计学意义。

实际上在在IDI.INF的返回值中除了M1、M2和M3，还包含了更多的一些其他信息m1.est、m2.est、m3.est及IDI.INF.GRAPH()作图用到的数据点。

IDI.INF.GRAPH(IDI.train.5Y)

以图形的形式显示IDI.INF的结果。

D表示预先指定时间点上新旧模型预测风险的差值，其取值应该在[-1,1]之间，如果新模型预测能力优于旧模型，则病例组倾向于正值而对照组倾向于负值。横坐标s表示D具体差值，纵坐标表示小于等于这个差值的概率。按照参考文献[Stat Med. 2013,32(14): 2430-2442]，粗实线表示发生结局事件的病例组的经验分布函数，细实线表示没有发生结束事件的对照组的经验分布函数。但从IDI的计算公式[(P_new|event-P_std|event)-(P_new|nonevent-P_std|nonevent)]上来看，细实线更像是病例组的经验分布函数，而粗实线则是对照组的经验分布函数。不过从作图数据（IDI.train.5Y结果中的point$FX、point$GX、point$cc）上看，作图用到的似乎是生存率而不是事件发生率，从这个角度上说粗实线确实是病例组的经验分布函数，细实线则是对照组的经验分布函数。抛开背后这些理论层次上的东西，结果就是：两条曲线下曲线下面积之差就是IDI，其值等于红色面积减去蓝色面积的差，即M1；两个黑点之间的垂直距离就是连续NRI，即M2；两个灰点之间的距离则是两个分布函数的中位数差值，即M3。

— END —

@_@