净重新分类指数NRI的计算
本文首发于公众号:医学和生信笔记
“
医学和生信笔记,专注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)
## Loading required package: 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': 232 obs. of 20 variables:## $ id : int 1 2 3 4 6 8 9 10 11 12 ...## $ time : int 400 4500 1012 1925 2503 2466 2400 51 3762 304 ...## $ status : int 2 0 2 2 2 2 2 2 2 2 ...## $ trt : int 1 1 1 1 2 2 1 2 2 2 ...## $ age : num 58.8 56.4 70.1 54.7 66.3 ...## $ sex : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...## $ ascites : int 1 0 0 0 0 0 0 1 0 0 ...## $ hepato : int 1 1 0 1 1 0 0 0 1 0 ...## $ spiders : int 1 1 0 1 0 0 1 1 1 1 ...## $ edema : num 1 0 0.5 0.5 0 0 0 1 0 0 ...## $ bili : num 14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...## $ chol : int 261 302 176 244 248 280 562 200 259 236 ...## $ albumin : num 2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...## $ copper : int 156 54 210 64 50 52 79 140 46 94 ...## $ alk.phos: num 1718 7395 516 6122 944 ...## $ ast : num 137.9 113.5 96.1 60.6 93 ...## $ trig : int 172 88 55 92 63 189 88 143 79 95 ...## $ platelet: int 190 221 151 183 NA 373 251 302 258 71 ...## $ protime : num 12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...## $ stage : int 4 3 4 4 3 3 2 4 4 4 ...
dim(dat) # 232 20
## [1] 232 20
然后就是准备计算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.valuesp.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%时,认为是重新分类。
上面的代码运行后结果是这样的:
UP and DOWN calculation: #of total, case, and control subjects at t0: 232 88 144
Reclassification Table for all subjects: NewStandard < 0.3 < 0.7 >= 0.7 < 0.3 135 4 0 < 0.7 1 31 4 >= 0.7 0 2 55
Reclassification Table for case: NewStandard < 0.3 < 0.7 >= 0.7 < 0.3 14 0 0 < 0.7 0 18 3 >= 0.7 0 1 52
Reclassification Table for control: NewStandard < 0.3 < 0.7 >= 0.7 < 0.3 121 4 0 < 0.7 1 13 1 >= 0.7 0 1 3
NRI estimation:Point estimates: EstimateNRI 0.001893939NRI+ 0.022727273NRI- -0.020833333Pr(Up|Case) 0.034090909Pr(Down|Case) 0.011363636Pr(Down|Ctrl) 0.013888889Pr(Up|Ctrl) 0.034722222
Now in bootstrap..
Point & Interval estimates: Estimate Std.Error Lower UpperNRI 0.001893939 0.027816095 -0.053995513 0.055354449NRI+ 0.022727273 0.021564394 -0.019801980 0.065789474NRI- -0.020833333 0.017312438 -0.058823529 0.007518797Pr(Up|Case) 0.034090909 0.019007629 0.000000000 0.072164948Pr(Down|Case) 0.011363636 0.010924271 0.000000000 0.039603960Pr(Down|Ctrl) 0.013888889 0.009334685 0.000000000 0.035211268Pr(Up|Ctrl) 0.034722222 0.014716046 0.006993007 0.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))*2p
## [1] 0.9457157
PredictABEL包
#install.packages("PredictABEL") #安装R包library(PredictABEL)
# 取出模型预测概率,这个包只能用预测概率计算p.std = mstd$fitted.valuesp.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) )
## _________________________________________## ## Reclassification table ## _________________________________________## ## Outcome: absent ## ## Updated Model## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified## [0,0.3) 121 4 0 3## [0.3,0.7) 1 13 1 13## [0.7,1] 0 1 3 25## ## ## Outcome: present ## ## Updated Model## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified## [0,0.3) 14 0 0 0## [0.3,0.7) 0 18 3 14## [0.7,1] 0 1 52 2## ## ## Combined Data ## ## Updated Model## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified## [0,0.3) 135 4 0 3## [0.3,0.7) 1 31 4 14## [0.7,1] 0 2 55 4## _________________________________________## ## 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')
UP and DOWN calculation: #of total, case, and control subjects at t0: 312 88 144
Reclassification Table for all subjects: NewStandard < 0.3 < 0.7 >= 0.7 < 0.3 202 7 0 < 0.7 13 53 6 >= 0.7 0 0 31
Reclassification Table for case: NewStandard < 0.3 < 0.7 >= 0.7 < 0.3 19 3 0 < 0.7 3 32 4 >= 0.7 0 0 27
Reclassification Table for control: NewStandard < 0.3 < 0.7 >= 0.7 < 0.3 126 3 0 < 0.7 5 7 2 >= 0.7 0 0 1
NRI estimation by KM estimator:
Point estimates: EstimateNRI 0.05377635NRI+ 0.03748660NRI- 0.01628974Pr(Up|Case) 0.07708938Pr(Down|Case) 0.03960278Pr(Down|Ctrl) 0.04256352Pr(Up|Ctrl) 0.02627378
Now in bootstrap..
Point & Interval estimates: Estimate Lower UpperNRI 0.05377635 -0.082230381 0.16058172NRI+ 0.03748660 -0.084245197 0.13231776NRI- 0.01628974 -0.030861213 0.06753616Pr(Up|Case) 0.07708938 0.000000000 0.19102291Pr(Down|Case) 0.03960278 0.000000000 0.15236016Pr(Down|Ctrl) 0.04256352 0.004671535 0.09863170Pr(Up|Ctrl) 0.02627378 0.006400463 0.05998424
结果的解读和logistic的一模一样。
survNRI包
# 安装R包devtools::install_github("mdbrown/survNRI")
加载R包并使用,还是用上面的pbc
数据集。
library(survNRI)
## Loading required package: 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.event NRI.nonevent NRI## KM 0.20445422 0.3187408 0.5231951## IPW 0.22424434 0.3273544 0.5515987## SmoothIPW 0.19645006 0.3144263 0.5108763## SEM 0.07478611 0.2632127 0.3379988## Combined 0.19633867 0.3143794 0.5107181## ## $CI## $CI$NRI.event## KM IPW SmoothIPW SEM Combined## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723## upperbound 0.44806768 0.47033936 0.44014214 0.2658309 0.4400496## ## $CI$NRI.nonevent## KM IPW SmoothIPW SEM Combined## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.1286426## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.6964549## ## $CI$NRI## KM IPW SmoothIPW SEM Combined## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409## upperbound 0.89306122 0.92464359 0.87970125 0.64253510 0.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 多平台发布
净重新分类指数NRI的计算相关推荐
- 【R语言】净重新分类指数NRI及其R语言实现
0. 引言 净重新分类指数 (net reclassification imporvement, NRI)是运用较多的比较新旧模型预测效果的方法 NRI关注的是在诊断截点处,通过考察使用新模型 后个体 ...
- R语言临床预测模型的评价指标与验证指标实战:净重新分类指数NRI(Net Reclassification Index, NRI)
R语言临床预测模型的评价指标与验证指标实战:净重新分类指数NRI(Net Classification Index, NRI) #净重新分类指数NRI
- R语言临床预测模型的评价指标与验证指标实战:自定义的净重新分类指数NRI(Net Reclassification Index, NRI)函数
R语言临床预测模型的评价指标与验证指标实战:自定义的净重新分类指数NRI(Net Reclassification Index, NRI)函数 目录
- R新旧模型、计算净重新分类指数(NRI)和整体鉴别指数(IDI)详解及实战
R计算净重新分类指数(NRI,net reclssification improvement)和整体鉴别指数(IDI,integrated discrimination improvement) NR ...
- 凝心聚力,共赢绿色计算新时代 ——2020绿色计算产业峰会在京召开
11月26日,由绿色计算产业联盟(Green Computing Consortium,简称GCC)主办的"2020绿色计算产业峰会"在北京召开.峰会按照"1+3+1&q ...
- Java SE 8新功能介绍:使用新的DateTime API计算时间跨度
使用Java SE 8新的DateTime API JSR 310-可以实现更清晰,可读且功能强大的编码. Java SE 8,JSR 310 在上一篇文章" 使用Streams API处理 ...
- 9、本金10000元存入银行,年利率是千分之三。每过一年,将本金和利息相加作为新的本金。计算五年后,获得的本金是多少?
import java.util.Scanner; public class Zuoye1 { /** * * 本金10000元存入银行,年利率是千分之三. * 每过一年,将本金和利息相加作为新的本金 ...
- 10、(加强版)本金10000元存入银行,年利率是千分之三。每过一年,将本金和利息相加作为新的本金。计算五年后,获得的本金是多少?
import java.util.Scanner; public class Zuoye1 { /** * * 升级版:本金10000元存入银行,年利率是千分之三. * 每过一年,将本金和利息相加作为 ...
- #今日论文推荐# 用GNN做CV三大任务的新骨干,同计算成本性能不输CNN、ViT与MLP|中科院华为诺亚开源
#今日论文推荐# 用GNN做CV三大任务的新骨干,同计算成本性能不输CNN.ViT与MLP|中科院&华为诺亚开源 用图神经网络(GNN)做CV的研究有不少,但通常是围绕点云数据做文章,少有直接 ...
最新文章
- 程序员必须知道的10大基础实用算法及其讲解
- 企业局域网离不开交换机/路由器/防火墙—Vecloud
- 2000/XP自动网络GHOST+全自动改IP
- 漫谈Google的Native Client(NaCl)技术
- Scala _10Actor Model
- 简单搜索 poj1321
- 多线程的那点儿事(之多线程调试)
- ENS与Cloudflare合作推出改进的ETH.LINK服务
- Itil v3 process model
- 华为服务器修改密码命令,华为IAD命令行配置命令
- 一文详解高精地图构建与SLAM感知优化建图策略
- 如何深入学习c语言,如何深入学习C语言?
- 修改了Excel默认打开方式后仍然使用WPS打开的解决办法
- win10怎么更改照片分辨率和大小?图片dpi修改方法
- 专访商汤科技联合创始人林达华丨一名AI人才,需要多少栽培?
- 计算机科学 vs 计算机技术
- 手机那点事!已有高人把常见的不常见的坑都给找出来了,我就随便转一下了
- C语言之大端模式与小端模式
- mysql事务的四个特点和实现原理
- mysqld_multi 没法stop