本內容為【科研私家菜】R語言機器學習與臨床預測模型系列課程
R小鹽準備介紹R語言機器學習與預測模型的學習筆記
你想要的R語言學習資料都在這里， 快來收藏關注【科研私家菜】
01 什么是凈重新分類指數(shù)NRI？
凈重新分類指數(shù)NRI（Net Reclassification Index, NRI）評價的是在兩模型采用最優(yōu)診斷截點進行預測時，與舊模型相比，使用新模型使得個體的預測結果得到改善的概率，包括兩個部分：
- 事件發(fā)生組：在R中以
NRI+表示 - 事件不發(fā)生組：在R中以
NRI-表示
以預測結局為患病和不患病的研究為例，當各自選定好了最佳閾值之后，我們可以得到一個混淆矩陣。
NRI用于我們在臨床上比較新舊模型的效能，比如說：平時我們都用心電圖評估心梗，現(xiàn)在有了新的指標肌鈣蛋白，我們想了解肌酐蛋白聯(lián)合心電圖評估心梗和單用心電圖哪個模型更好。我們知道一個模型建立后能產生一個切點，把一個人群分切陽性組和陰性組，新模型的建立后陽性組和陰性組肯定和原來模型不同，假陽性和假陰性也不同，NRI主要比較的是重新分布這類人群的真陽性和真陰性情況。
方法1 nricens包
# install.packages('nricens')
library('nricens')
NRIb = nribin(event = labels, p.std = pred1, p.new = pred2, updown = 'category', cut=0.5)
# event為標簽；cut為模型閾值
# p.std舊模型的預測值，p.new為新模型的預測值，p.std和p.new都是連續(xù)變量，大小介于0和1之間。
方法2 PredictABEL包
# install.packages('PredictABEL')
library('PredictABEL')
reclassification(data=data, cOutcome = 7, predrisk1 = pred1, predrisk2 = pred2, cutoff = c(0, 0.5,1))
# data為數(shù)據(jù)集；cOutcome = 7表示數(shù)據(jù)及中的第七列為標簽列；cutoff為模型閾值
# pred1舊模型的預測值，pred2為新模型的預測值，pred1和pred2都是連續(xù)變量，大小介于0和1之間。
02 R語言實現(xiàn)
library(nricens)
library(survival)
library(foreign)
bc <- read.spss("survival agec.sav",
use.value.labels=F, to.data.frame=T)
bc <- na.omit(bc)#刪除缺失值
names(bc)
time<-bc$time#生成時間變量
status<-bc$status#生成結局變量
j1<-as.matrix(subset(bc,select = c(age,pathsize,pr,er)))#舊模型變量，要以矩陣形式，不能有缺失值，不然會報錯
j2<-as.matrix(subset(bc,select = c(age,pathsize,pr,er,ln_yesno))) #新模型變量，要以矩陣形式，不能有缺失值，不然會報錯
mod.std<-coxph(Surv(time,status)~ .,data.frame(time, status,j1),x=TRUE)#生成舊模型COX回歸函數(shù)
mod.new<-coxph(Surv(time,status)~ .,data.frame(time, status,j2),x=TRUE) #生成新模型COX回歸函數(shù)
p.std = get.risk.coxph(mod.std, t0=48)#生成預測函數(shù)
p.new = get.risk.coxph(mod.new, t0=48)#生成預測函數(shù)
nricens(mdl.std = mod.std, mdl.new = mod.new, t0 = 48, cut = c(0.2, 0.4),
niter = 100,alpha = 0.05,updown = 'category')#cut分臨床切點，分為高中低風險，niter為重抽樣次數(shù)，category為分類的意思，alpha為置信區(qū)間
我們要關注NRI+這個值，它是正值表示，新模型優(yōu)于舊模型，最后生成圖表
還按分類變量及廣義線性模型來處理，要先設置一下時間變量，和結局變量，其他操作一樣的
bc= bc[ bc$time > 48| (bc$time < 48 & bc$status == 1), ]#重新定義一下數(shù)據(jù)
status1= ifelse(bc$time < 48 & bc$status == 1, 1, 0)#重新定義下結局變量
j3<-as.matrix(subset(bc,select = c(age,pathsize,pr,er)))
j4<-as.matrix(subset(bc,select = c(age,pathsize,pr,er,ln_yesno)))
mstd1= glm(status1 ~ ., binomial(logit), data.frame(status1, j3), x=TRUE)
mnew1= glm(status1 ~ ., binomial(logit), data.frame(status1, j4), x=TRUE)
nribin(mdl.std=mstd1, mdl.new = mnew1, cut = c(0.2, 0.4),
niter = 1000, updown = 'category')
效果如下：
關注R小鹽，關注科研私家菜(VX_GZH: SciPrivate)，有問題請聯(lián)系R小鹽。讓我們一起來學習 R語言機器學習與臨床預測模型
