|
凈重分類改善指數（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年內發(fā)生結局事件的的預測概率(患病風險）Pstd.e3Y.test<-1-summary(survfit(BC.std,newdata=testset),times=365.25*3)$surv[1,] #模型BC.std對測試集中每個個體在3年內發(fā)生結局事件的的預測概率(患病風險） Pnew.e3Y.train<-1-summary(survfit(BC.new,newdata=trainset),times=365.25*3)$surv[1,] #模型BC.new對訓練集中每個個體在3年內發(fā)生結局事件的的預測概率(患病風險） Pnew.e3Y.test<-1-summary(survfit(BC.new,newdata=testset),times=365.25*3)$surv[1,] #模型BC.new對測試集中每個個體在3年內發(fā)生結局事件的的預測概率(患病風險） 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()可將因子變量設置為啞變量并保持連續(xù)變量不變來創(chuàng)建模型矩陣,，[,-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指定結局狀態(tài),，0為刪失,，1為發(fā)生結局事件；mdl.std和mdl.new分別指定標準模型和新模型的coxph或survreg對象,，在擬合Cox或參數生存模型時其x參數需要設置為x=TRUE（預測變量需要從glm對象中提?。?/span>z.std和z.new分別指定標準模型和新模型的預測變量矩陣,，不允許有因子,、字符和缺失值；p.std和p.new分別指定標準模型和新模型的預測風險向量,；t0指定評估時間點,；updown指定確定up和down的方法，默認為"category",，計算的是基于風險分類的NRI,，取值"diff"，將計算基于風險差的NRI,；cut指定一個標量或向量,，作為確定up和down的預測風險（風險類別或風險差）的界值，標量（1個值）用于二分類或風險差,，向量（多個值）可用于多分類,。連續(xù)性NRI是風險差的一種特殊類型,，可以通過指定updown="diff"且cut=0來實現；point.method指定點估計方法,，"km" 采用Kaplan-Meier (KM)法,，"ipw"采用逆概率加權 (IPW) 估計；niter自主抽樣的次數,，具體使用的是百分位自助式方法,，用于計算置信區(qū)間，默認1000次,，若取值為0則不計算置信區(qū)間,；,alpha指定α值，1-α置信區(qū)間將被計算,；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,，該模式下自助抽樣時并未執(zhí)行預測模型的擬合,，可用于外部數據或交叉驗證的驗證研究。這種模式下,，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年內發(fā)生結局事件的病例159人,，沒有發(fā)生結局事件的對照組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法估計的置信區(qū)間,。
標準模型和新模型對測試集3年發(fā)生結局事件的的分類變化情況： 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年預測結果按連續(xù)性NRI重分類變化情況： 在《模型預測能力的比較：IDI[綜合判別改善指數]》一文中，在計算兩個Cox模型IDI時結果同時輸出了連續(xù)NRI的結果,。我們發(fā)現IDI.INF::survIDINRI與nricens::nricens計算的連續(xù)NRI結果并不一致,。原因有①導入數據劃分訓練集和測試集時使用了不同的隨機種子，所以建模時采用的數據集并不相同；②survIDINRI包中,，連續(xù)性NRI定義為1/2 NRI(>0),，即普通NRI的一半。NRI一般定義為事件組NRI和對照組NRI的和,，取值一半相當于把兩個子組進行平均,，這樣的NRI的值域就定義在[-1,1]之間了。NRI(>0)表示連續(xù)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 > 《待分類》
