本文首发于公众号：医学和生信笔记
“
医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
NRI，net reclassification index，净重新分类指数，是用来比较模型准确度的，这个概念有点难理解，但是非常重要，在临床研究中非常常见，是评价模型的一大利器！
在R语言中有很多包可以计算NRI，但是能同时计算logistic回归和cox回归的只有nricens包，PredictABEL可以计算logistic模型的净重分类指数，survNRI可以计算cox模型的净重分类指数。
- logistic的NRI
- nricens包
- PredictABEL包
- 生存分析的NRI
- nricens包
- survNRI包
logistic的NRI
nricens包
#install.packages("nricens")#安装R包
library(nricens)
##Loadingrequiredpackage:survival
使用survival包中的pbc数据集用于演示，这是一份关于原发性硬化性胆管炎的数据，其实是一份用于生存分析的数据，是有时间变量的，但是这里我们用于演示logistic回归，只要不使用time这一列就可以了。
library(survival)
#只使用部分数据
dat=pbc[1:312,]
dat=dat[dat$time>2000|(dat$time<2000&dat$status==2),]
str(dat)#数据长这样
##'data.frame':232obs.of20variables:
##$id:int1234689101112...
##$time:int400450010121925250324662400513762304...
##$status:int2022222222...
##$trt:int1111221222...
##$age:num58.856.470.154.766.3...
##$sex:Factorw/2levels"m","f":2212222222...
##$ascites:int1000000100...
##$hepato:int1101100010...
##$spiders:int1101001111...
##$edema:num100.50.5000100...
##$bili:num14.51.11.41.80.80.33.212.61.43.6...
##$chol:int261302176244248280562200259236...
##$albumin:num2.64.143.482.543.9843.082.744.163.52...
##$copper:int15654210645052791404694...
##$alk.phos:num171873955166122944...
##$ast:num137.9113.596.160.693...
##$trig:int17288559263189881437995...
##$platelet:int190221151183NA37325130225871...
##$protime:num12.210.61210.311111111.51213.6...
##$stage:int4344332444...
dim(dat)#23220
##[1]23220
然后就是准备计算NRI所需要的各个参数。
#定义结局事件，0是存活，1是死亡
event=ifelse(dat$time<2000&dat$status==2,1,0)
#两个只由预测变量组成的矩阵
z.std=as.matrix(subset(dat,select=c(age,bili,albumin)))
z.new=as.matrix(subset(dat,select=c(age,bili,albumin,protime)))
#建立2个模型
mstd=glm(event~age+bili+albumin,family=binomial(),data=dat,x=TRUE)
mnew=glm(event~age+bili+albumin+protime,family=binomial(),data=dat,x=TRUE)
#取出模型预测概率
p.std=mstd$fitted.values
p.new=mnew$fitted.values
然后就是计算NRI，对于二分类变量，使用nribin()函数，这个函数提供了3种参数使用组合，任选一种都可以计算出来（结果一样），以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。
#这3种方法算出来都是一样的结果
#两个模型
nribin(mdl.std=mstd,mdl.new=mnew,
cut=c(0.3,0.7),
niter=500,
updown='category')
#结果变量+两个只有预测变量的矩阵
nribin(event=event,z.std=z.std,z.new=z.new,
cut=c(0.3,0.7),
niter=500,
updown='category')
##结果变量+两个模型得到的预测概率
nribin(event=event,p.std=p.std,p.new=p.new,
cut=c(0.3,0.7),
niter=500,
updown='category')
其中，cut是判断风险高低的阈值，我们使用了0.3,0.7，代表0-0.3是低风险，0.3-0.7是中风险，0.7-1是高风险，这个阈值是自己设置的，大家根据经验或者文献设置即可。
niter是使用bootstrap法进行重抽样的次数，默认是1000，大家可以自己设置。
updown参数，当设置为category时，表示低、中、高风险这种方式；当设置为diff时，此时cut的取值只能设置1个，比如设置0.2，即表示当新模型预测的风险和旧模型相差20%时，认为是重新分类。
上面的代码运行后结果是这样的：
UPandDOWNcalculation:
#oftotal,case,andcontrolsubjectsatt0:23288144
ReclassificationTableforallsubjects:
New
Standard<0.3<0.7>=0.7
<0.313540
<0.71314
>=0.70255
ReclassificationTableforcase:
New
Standard<0.3<0.7>=0.7
<0.31400
<0.70183
>=0.70152
ReclassificationTableforcontrol:
New
Standard<0.3<0.7>=0.7
<0.312140
<0.71131
>=0.7013
NRIestimation:
Pointestimates:
Estimate
NRI0.001893939
NRI+0.022727273
NRI--0.020833333
Pr(Up|Case)0.034090909
Pr(Down|Case)0.011363636
Pr(Down|Ctrl)0.013888889
Pr(Up|Ctrl)0.034722222
Nowinbootstrap..
Point&Intervalestimates:
EstimateStd.ErrorLowerUpper
NRI0.0018939390.027816095-0.0539955130.055354449
NRI+0.0227272730.021564394-0.0198019800.065789474
NRI--0.0208333330.017312438-0.0588235290.007518797
Pr(Up|Case)0.0340909090.0190076290.0000000000.072164948
Pr(Down|Case)0.0113636360.0109242710.0000000000.039603960
Pr(Down|Ctrl)0.0138888890.0093346850.0000000000.035211268
Pr(Up|Ctrl)0.0347222220.0147160460.0069930070.066176471
首先是3个混淆矩阵，第一个是全体的，第2个是case（结局为1）组的，第3个是control（结局为2）组的，有了这3个矩阵，我们可以自己计算净重分类指数。
看case组：
净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273
再看control组：
净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333
相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 – 3/144 ≈ 0.000315657
再往下是不做bootstrap时得到的估计值，其中NRI就是绝对净重分类指数，NRI+是case组的净重分类指数，NRI-是control组的净重分类指数（和我们计算的一样哦），最后是做了500次bootstrap后得到的估计值，并且有标准误和可信区间。
最后还会得到一张图： 
这张图中的虚线对应的坐标，就是我们在cut中设置的阈值，这张图对应的是上面结果中的第一个混淆矩阵，反应的是总体的情况，case是结果为1的组，也就是发生结局的组，control是结果为0的组，也就是未发生结局的组。
P值没有直接给出，但是可以自己计算。
#计算P值
z<-abs(0.001893939/0.027816095)
p<-(1-pnorm(z))*2
p
##[1]0.9457157
PredictABEL包
#install.packages("PredictABEL")#安装R包
library(PredictABEL)
#取出模型预测概率，这个包只能用预测概率计算
p.std=mstd$fitted.values
p.new=mnew$fitted.values
然后就是计算NRI：
dat$event<-event
reclassification(data=dat,
cOutcome=21,#结果变量在哪一列
predrisk1=p.std,
predrisk2=p.new,
cutoff=c(0,0.3,0.7,1)
)
##_________________________________________
##
##Reclassificationtable
##_________________________________________
##
##Outcome:absent
##
##UpdatedModel
##InitialModel[0,0.3)[0.3,0.7)[0.7,1]%reclassified
##[0,0.3)121403
##[0.3,0.7)113113
##[0.7,1]01325
##
##
##Outcome:present
##
##UpdatedModel
##InitialModel[0,0.3)[0.3,0.7)[0.7,1]%reclassified
##[0,0.3)14000
##[0.3,0.7)018314
##[0.7,1]01522
##
##
##CombinedData
##
##UpdatedModel
##InitialModel[0,0.3)[0.3,0.7)[0.7,1]%reclassified
##[0,0.3)135403
##[0.3,0.7)131414
##[0.7,1]02554
##_________________________________________
##
##NRI(Categorical)[95%CI]:0.0019[-0.0551-0.0589];p-value:0.94806
##NRI(Continuous)[95%CI]:0.0391[-0.2238-0.3021];p-value:0.77048
##IDI[95%CI]:0.0044[-0.0037-0.0126];p-value:0.28396
结果得到的是相加净重分类指数，还给出了IDI和P值。两个包算是各有优劣吧，大家可以自由选择。
生存分析的NRI
还是使用survival包中的pbc数据集用于演示，这次要构建cox回归模型，因此我们要使用time这一列了。
nricens包
library(nricens)
library(survival)
dat<-pbc[1:312,]
dat$status<-ifelse(dat$status==2,1,0)#0表示活着，1表示死亡
然后准备所需参数：
#两个只由预测变量组成的矩阵
z.std=as.matrix(subset(dat,select=c(age,bili,albumin)))
z.new=as.matrix(subset(dat,select=c(age,bili,albumin,protime)))
#建立2个cox模型
mstd<-coxph(Surv(time,status)~age+bili+albumin,data=dat,x=TRUE)
mnew<-coxph(Surv(time,status)~age+bili+albumin+protime,data=dat,x=TRUE)
#计算在2000天的模型预测概率，这一步不要也行，看你使用哪些参数
p.std<-get.risk.coxph(mstd,t0=2000)
p.new<-get.risk.coxph(mnew,t0=2000)
计算NRI：
nricens(mdl.std=mstd,mdl.new=mnew,
t0=2000,
cut=c(0.3,0.7),
niter=1000,
updown='category')
UPandDOWNcalculation:
#oftotal,case,andcontrolsubjectsatt0:31288144
ReclassificationTableforallsubjects:
New
Standard<0.3<0.7>=0.7
<0.320270
<0.713536
>=0.70031
ReclassificationTableforcase:
New
Standard<0.3<0.7>=0.7
<0.31930
<0.73324
>=0.70027
ReclassificationTableforcontrol:
New
Standard<0.3<0.7>=0.7
<0.312630
<0.7572
>=0.7001
NRIestimationbyKMestimator:
Pointestimates:
Estimate
NRI0.05377635
NRI+0.03748660
NRI-0.01628974
Pr(Up|Case)0.07708938
Pr(Down|Case)0.03960278
Pr(Down|Ctrl)0.04256352
Pr(Up|Ctrl)0.02627378
Nowinbootstrap..
Point&Intervalestimates:
EstimateLowerUpper
NRI0.05377635-0.0822303810.16058172
NRI+0.03748660-0.0842451970.13231776
NRI-0.01628974-0.0308612130.06753616
Pr(Up|Case)0.077089380.0000000000.19102291
Pr(Down|Case)0.039602780.0000000000.15236016
Pr(Down|Ctrl)0.042563520.0046715350.09863170
Pr(Up|Ctrl)0.026273780.0064004630.05998424
结果的解读和logistic的一模一样。
survNRI包
#安装R包
devtools::install_github("mdbrown/survNRI")
加载R包并使用，还是用上面的pbc数据集。
library(survNRI)
##Loadingrequiredpackage:MASS
library(survival)
#使用部分数据
dat<-pbc[1:312,]
dat$status<-ifelse(dat$status==2,1,0)#0表示活着，1表示死亡
res<-survNRI(time="time",event="status",
model1=c("age","bili","albumin"),#模型1的自变量
model2=c("age","bili","albumin","protime"),#模型2的自变量
data=dat,
predict.time=2000,#预测的时间点
method="all",
bootMethod="normal",
bootstraps=500,
alpha=.05)
查看结果，$estimates给出了不同组的NRI以及总的NRI，包括了使用不同方法（KM/IPW/SmoothIPW/SEM/Combined）得到的结果；$CI给出了可信区间。
res
##$estimates
##NRI.eventNRI.noneventNRI
##KM0.204454220.31874080.5231951
##IPW0.224244340.32735440.5515987
##SmoothIPW0.196450060.31442630.5108763
##SEM0.074786110.26321270.3379988
##Combined0.196338670.31437940.5107181
##
##$CI
##$CI$NRI.event
##KMIPWSmoothIPWSEMCombined
##lowerbound-0.03915924-0.02185068-0.04724202-0.1162587-0.0473723
##upperbound0.448067680.470339360.440142140.26583090.4400496
##
##$CI$NRI.nonevent
##KMIPWSmoothIPWSEMCombined
##lowerbound0.13171080.13963150.12866850.086389330.1286426
##upperbound0.71022510.73932160.69663410.514822120.6964549
##
##$CI$NRI
##KMIPWSmoothIPWSEMCombined
##lowerbound-0.05112533-0.04569046-0.05439863-0.04132364-0.05443409
##upperbound0.893061220.924643590.879701250.642535100.87953153
##
##
##$bootMethod
##[1]"normal"
##
##$predict.time
##[1]2000
##
##$alpha
##[1]0.05
##
##attr(,"class")
##[1]"survNRI"
OK，这就是NRI的计算，除此之外，随机森林、决策树、lasso回归、SVM等，这些模型，都是可以计算的NRI的，后面会继续介绍。大家如果有问题欢迎在评论区留言。
本文首发于公众号：医学和生信笔记
“
医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
本文由 mdnice 多平台发布