本內(nèi)容為【科研私家菜】R語言機(jī)器學(xué)習(xí)與臨床預(yù)測模型系列課程
R小鹽準(zhǔn)備介紹R語言機(jī)器學(xué)習(xí)與預(yù)測模型的學(xué)習(xí)筆記
你想要的R語言學(xué)習(xí)資料都在這里然遏， 快來收藏關(guān)注【科研私家菜】
01 什么是凈重新分類指數(shù)NRI？
凈重新分類指數(shù)NRI（Net Reclassification Index, NRI）評價(jià)的是在兩模型采用最優(yōu)診斷截點(diǎn)進(jìn)行預(yù)測時(shí)胚嘲，與舊模型相比畏鼓，使用新模型使得個體的預(yù)測結(jié)果得到改善的概率盹廷，包括兩個部分：
- 事件發(fā)生組：在R中以
NRI+表示 - 事件不發(fā)生組：在R中以
NRI-表示
以預(yù)測結(jié)局為患病和不患病的研究為例派昧，當(dāng)各自選定好了最佳閾值之后，我們可以得到一個混淆矩陣谬以。
NRI用于我們在臨床上比較新舊模型的效能，比如說：平時(shí)我們都用心電圖評估心梗由桌，現(xiàn)在有了新的指標(biāo)肌鈣蛋白为黎，我們想了解肌酐蛋白聯(lián)合心電圖評估心梗和單用心電圖哪個模型更好。我們知道一個模型建立后能產(chǎn)生一個切點(diǎ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為標(biāo)簽；cut為模型閾值
# p.std舊模型的預(yù)測值捌斧，p.new為新模型的預(yù)測值笛质，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ù)及中的第七列為標(biāo)簽列；cutoff為模型閾值
# pred1舊模型的預(yù)測值姓迅，pred2為新模型的預(yù)測值敲霍，pred1和pred2都是連續(xù)變量，大小介于0和1之間丁存。
02 R語言實(shí)現(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#生成時(shí)間變量
status<-bc$status#生成結(jié)局變量
j1<-as.matrix(subset(bc,select = c(age,pathsize,pr,er)))#舊模型變量肩杈，要以矩陣形式，不能有缺失值解寝，不然會報(bào)錯
j2<-as.matrix(subset(bc,select = c(age,pathsize,pr,er,ln_yesno))) #新模型變量扩然，要以矩陣形式，不能有缺失值聋伦，不然會報(bào)錯
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)#生成預(yù)測函數(shù)
p.new = get.risk.coxph(mod.new, t0=48)#生成預(yù)測函數(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分臨床切點(diǎn)夫偶，分為高中低風(fēng)險(xiǎn)，niter為重抽樣次數(shù)嘉抓，category為分類的意思索守，alpha為置信區(qū)間
我們要關(guān)注NRI+這個值晕窑，它是正值表示抑片，新模型優(yōu)于舊模型，最后生成圖表
還按分類變量及廣義線性模型來處理杨赤，要先設(shè)置一下時(shí)間變量敞斋，和結(jié)局變量截汪，其他操作一樣的
bc= bc[ bc$time > 48| (bc$time < 48 & bc$status == 1), ]#重新定義一下數(shù)據(jù)
status1= ifelse(bc$time < 48 & bc$status == 1, 1, 0)#重新定義下結(jié)局變量
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')
效果如下：
關(guān)注R小鹽，關(guān)注科研私家菜(VX_GZH: SciPrivate)植捎，有問題請聯(lián)系R小鹽衙解。讓我們一起來學(xué)習(xí) R語言機(jī)器學(xué)習(xí)與臨床預(yù)測模型
