|
净重分类改善指数（Net Reclassification Improvement Index,NRI）的基本概念参见《两个logistic模型预测能力的比较：NRI[净重分类改善指数]》。与之不同的是，在生存数据一般会存在删失的情况，两个生存模型NRI的计算需要预先指定一个时间点。 library(survival) GBSG<-gbsg GBSG$grade<-factor(GBSG$grade) library(ggplot2);library(lattice) set.seed(20240229) split<-createDataPartition(GBSG$status,p=0.7,list=FALSE) split<-as.vector(split) trainset<-GBSG[split,] testset<-GBSG[-split,] 【2】建立Cox模型并获得计算NRI需要的一些参数 以age、meno、size、nodes、grade五个变量建立的模型为旧模型，利用全子集回归获得4个变量size、nodes、pgr、hormon建立的模型为新模型。 BC.std<-coxph(Surv(rfstime,status)~age+meno+size+nodes+grade,data=trainset,x=TRUE) BC.new<-coxph(Surv(rfstime,status)~size+nodes+pgr+hormon,data=trainset,x=TRUE) Pstd.e3Y.train<-1-summary(survfit(BC.std,newdata=trainset),times=365.25*3)$surv[1,] #模型BC.std对训练集中每个个体在3年内发生结局事件的的预测概率(患病风险）Pstd.e3Y.test<-1-summary(survfit(BC.std,newdata=testset),times=365.25*3)$surv[1,] #模型BC.std对测试集中每个个体在3年内发生结局事件的的预测概率(患病风险） Pnew.e3Y.train<-1-summary(survfit(BC.new,newdata=trainset),times=365.25*3)$surv[1,] #模型BC.new对训练集中每个个体在3年内发生结局事件的的预测概率(患病风险） Pnew.e3Y.test<-1-summary(survfit(BC.new,newdata=testset),times=365.25*3)$surv[1,] #模型BC.new对测试集中每个个体在3年内发生结局事件的的预测概率(患病风险） library(nricens) Pstd.e3Y.train<-get.risk.coxph(BC.std, 365.25*3) Pnew.e3Y.train<-get.risk.coxph(BC.new, 365.25*3) Xstd.train<-model.matrix(~age+meno+size+grade+nodes,data=trainset)[,-1] #训练集中模型BC.std的预测变量矩阵。该矩阵中含有一个三水平分类变量grade，后面我们计算NRI的函数nricens::nricens不允许有因子变量，通过model.matrix()可将因子变量设置为哑变量并保持连续变量不变来创建模型矩阵，[,-1]是从 model.matrix 的输出中模型矩阵删除Intercept项。直接使用命令as.matrix() 只能将因子变量的赋值转换为数值格式,后面运行程序时会报错 Xstd.test<-model.matrix(~age+meno+size+grade+nodes,data=testset)[,-1] #测试集中模型BC.std预测变量矩阵 Xnew.train<-as.matrix(subset(trainset,select=c(size,nodes,pgr,hormon))) #训练集中模型BC.new的预测变量矩阵 nricens(time=NULL,event=NULL,mdl.std=NULL,mdl.new=NULL,z.std=NULL,z.new=NULL,p.std=NULL,p.new=NULL,t0=NULL,updown="category",cut=NULL,point.method="km", niter=1000,alpha=0.05,msg=TRUE) #time指定随访时间；event指定结局状态，0为删失，1为发生结局事件；mdl.std和mdl.new分别指定标准模型和新模型的coxph或survreg对象，在拟合Cox或参数生存模型时其x参数需要设置为x=TRUE（预测变量需要从glm对象中提取）；z.std和z.new分别指定标准模型和新模型的预测变量矩阵，不允许有因子、字符和缺失值；p.std和p.new分别指定标准模型和新模型的预测风险向量；t0指定评估时间点；updown指定确定up和down的方法，默认为"category"，计算的是基于风险分类的NRI，取值"diff"，将计算基于风险差的NRI；cut指定一个标量或向量，作为确定up和down的预测风险（风险类别或风险差）的界值，标量（1个值）用于二分类或风险差，向量（多个值）可用于多分类。连续性NRI是风险差的一种特殊类型，可以通过指定updown="diff"且cut=0来实现；point.method指定点估计方法，"km" 采用Kaplan-Meier (KM)法，"ipw"采用逆概率加权 (IPW) 估计；niter自主抽样的次数，具体使用的是百分位自助式方法，用于计算置信区间，默认1000次，若取值为0则不计算置信区间；,alpha指定α值，1-α置信区间将被计算；msg指定是否显示计算过程 #指定以下任何一组参数都可以计算NRI：(mdl.std,mdl.new);(time,event,z.std,z.new);(time,event,p.std,p.new)。指定参数mdl.std和mdl.new：生存时间、结局事件、标准模型和新模型的预测变量都将从模型中提取。但如果也指定参数time和event，则使用time和event指定的生存时间和结局事件而不是从coxph或survreg对象中提取；指定参数time、event、z.std和z.new，直接利用预测变量和生存时间、结局变量拟合一个标准模型和一个新的预测模型。这种模式下，z.std和z.new需要使用同一个数据集；指定参数time、event,p.std和p.new，直接使用预测风险来计算NRI，该模式下自助抽样时并未执行预测模型的拟合，可用于外部数据或交叉验证的验证研究。这种模式下，p.std和p.new也需要使用同一个数据集 set.seed(20240229) nricens(mdl.std=BC.std,mdl.new=BC.new,t0=365.25*3,updown="category",cut=c(0.2,0.5)) #0-0.2为低风险；0.2-0.5中风险；0.5-1高风险  Total、Case、Control分别为481、159和234，表示全部研究对象为481人，3年内发生结局事件的病例159人，没有发生结局事件的对照组234人。Total中除了Case、Control外，还有随访到3年时删失的人数481-159-234=88人。 三个表分别是针对所有研究对象、病例组和对照组两个模型的分类情况。因为涉及到删失数据，直接通过表格中的向上和向下数据来估计NRI并不合适。NRI点估计有Kaplan-Meier (KM)法、逆概率加权法（IPW）、半参数法（SEM）、平滑逆概率加权法（SmoothIPW）等，不同的方法计算的结果会有差异。nricens函数中提供了KM和IPW法，可通过参数point.method进行设置。  NRI点估计结果NRI、NRI+、NRI-分别对应的就是相加NRI、病例组NRI和对照组NRI。 在预测3年的生存结局上，与标准模型相比，新模型重分类正确的比例提高20.22%，其中对阳性结局重分类正确的比例提高13.74%，对阴性结局重分类正确的比例提高6.48%。 NRI+=Pr(Up|Case)-Pr(Down|Case)；NRI-=Pr(Down|Ctrl)-Pr(Up|Ctrl)。 注：Pr(Up|Case)表示利用新模型对病例组进行重分类，风险类别提高的那些病例所占的比例；Pr(Down|Case)表示利用新模型对病例组进行重分类，风险类别降低的那些病例所占的比例；Pr(Down|Ctrl)表示利用新模型对对照组进行重分类，风险类别降低的那些病例所占的比例；Pr(Up|Ctrl)表示利用新模型对对照组进行重分类，风险类别升高的那些病例所占的比例。 其后是利用Bootstrap法估计的置信区间。
标准模型和新模型对测试集3年发生结局事件的的分类变化情况： set.seed(20240229) nricens(time=testset$rfstime,event=testset$status,p.std=Pstd.e3Y.test,p.new=Pnew.e3Y.test,t0=365.25*3,updown="category",cut=c(0.2,0.5))       新模型对训练集的3年预测结果按差值15%进行二分类的重分类变化情况： set.seed(20240229) nricens(mdl.std=BC.std,mdl.new=BC.new,t0=365.25*3,updown="diff",cut=0.15) #updown设置为diff时，cut的取值只能1个，本例设置0.15表示当新模型预测的风险和旧模型相差0.15时，认为是重新分类   全部研究对象为481人，病例组159人，对照组234人；新模型比标准模型预测风险高0.15以上共35人，其中病例组14人，对照组15人；新模型比标准模型预测风险低0.15以上的共54人，其中病例组7人，对照组35人。 新模型对训练集的5年预测结果按连续性NRI重分类变化情况： 在《模型预测能力的比较：IDI[综合判别改善指数]》一文中，在计算两个Cox模型IDI时结果同时输出了连续NRI的结果。我们发现IDI.INF::survIDINRI与nricens::nricens计算的连续NRI结果并不一致。原因有①导入数据划分训练集和测试集时使用了不同的随机种子，所以建模时采用的数据集并不相同；②survIDINRI包中，连续性NRI定义为1/2 NRI(>0)，即普通NRI的一半。NRI一般定义为事件组NRI和对照组NRI的和，取值一半相当于把两个子组进行平均，这样的NRI的值域就定义在[-1,1]之间了。NRI(>0)表示连续NRI（无类别NRI），NRI(0.5)表示两类别NRI，类别划分界值为0.5，而NRI(0.3,0.8)表示三分类NRI，分类界值是0.3和0.8；③两生存模型NRI估计有不同的方法，KM、IPW、SEM、SmoothIPW，采用不同的方法结果会有出入。本例如果采用IPW法，与survIDINRI计算的结果就非常接近了（NRI：survIDINRI≈1/2nricens）： nricens(mdl.std=BC.std,mdl.new=BC.new,t0=365.25*5,updown="diff",cut=0,point.method="ipw")
— END — |
|
|
来自： Memo_Cleon > 《待分类》
