R语言机器学习与临床预测模型25--净重新分类指数

本内容为【科研私家菜】R语言机器学习与临床预测模型系列课程

R小盐准备介绍R语言机器学习与预测模型的学习笔记

你想要的R语言学习资料都在这里， 快来收藏关注【科研私家菜】

---

01 什么是净重新分类指数NRI？

净重新分类指数NRI（Net Reclassification Index, NRI）评价的是在两模型采用最优诊断截点进行预测时，与旧模型相比，使用新模型使得个体的预测结果得到改善的概率，包括两个部分：

1. 事件发生组：在R中以NRI+表示
2. 事件不发生组：在R中以NRI-表示
   以预测结局为患病和不患病的研究为例，当各自选定好了最佳阈值之后，我们可以得到一个混淆矩阵。
   NRI用于我们在临床上比较新旧模型的效能，比如说：平时我们都用心电图评估心梗，现在有了新的指标肌钙蛋白，我们想了解肌酐蛋白联合心电图评估心梗和单用心电图哪个模型更好。我们知道一个模型建立后能产生一个切点，把一个人群分切阳性组和阴性组，新模型的建立后阳性组和阴性组肯定和原来模型不同，假阳性和假阴性也不同，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都是连续变量，大小介于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为数据集；cOutcome = 7表示数据及中的第七列为标签列；cutoff为模型阈值
# pred1旧模型的预测值，pred2为新模型的预测值，pred1和pred2都是连续变量，大小介于0和1之间。

---

02 R语言实现

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回归函数
mod.new<-coxph(Surv(time,status)~ .,data.frame(time, status,j2),x=TRUE) #生成新模型COX回归函数
p.std = get.risk.coxph(mod.std, t0=48)#生成预测函数
p.new = get.risk.coxph(mod.new, t0=48)#生成预测函数
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为重抽样次数，category为分类的意思，alpha为置信区间

我们要关注NRI+这个值，它是正值表示，新模型优于旧模型，最后生成图表

还按分类变量及广义线性模型来处理，要先设置一下时间变量，和结局变量，其他操作一样的

bc= bc[ bc$time > 48| (bc$time < 48 & bc$status == 1), ]#重新定义一下数据
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)，有问题请联系R小盐。让我们一起来学习 R语言机器学习与临床预测模型