净重新分类指数NRI的计算

这篇具有很好参考价值的文章主要介绍了净重新分类指数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 == 210)

# 两个只由预测变量组成的矩阵
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%时,认为是重新分类。

上面的代码运行后结果是这样的:

UP and DOWN calculation:
  #of total, case, and control subjects at t0:  232 88 144

  Reclassification Table for all subjects:
        New
Standard < 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:
        New
Standard < 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:
        New
Standard < 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:
                  Estimate
NRI            0.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

Now in bootstrap..

Point & Interval estimates:
                  Estimate   Std.Error        Lower       Upper
NRI            0.001893939 0.027816095 -0.053995513 0.055354449
NRI+           0.022727273 0.021564394 -0.019801980 0.065789474
NRI-          -0.020833333 0.017312438 -0.058823529 0.007518797
Pr(Up|Case)    0.034090909 0.019007629  0.000000000 0.072164948
Pr(Down|Case)  0.011363636 0.010924271  0.000000000 0.039603960
Pr(Down|Ctrl)  0.013888889 0.009334685  0.000000000 0.035211268
Pr(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后得到的估计值,并且有标准误和可信区间。

最后还会得到一张图: 净重新分类指数NRI的计算

这张图中的虚线对应的坐标,就是我们在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)
                 )
##  _________________________________________
##  
##      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==210# 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.30.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:
        New
Standard < 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:
        New
Standard < 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:
        New
Standard < 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:
                Estimate
NRI           0.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

Now in bootstrap..

Point & Interval estimates:
                Estimate        Lower      Upper
NRI           0.05377635 -0.082230381 0.16058172
NRI+          0.03748660 -0.084245197 0.13231776
NRI-          0.01628974 -0.030861213 0.06753616
Pr(Up|Case)   0.07708938  0.000000000 0.19102291
Pr(Down|Case) 0.03960278  0.000000000 0.15236016
Pr(Down|Ctrl) 0.04256352  0.004671535 0.09863170
Pr(Up|Ctrl)   0.02627378  0.006400463 0.05998424
净重新分类指数NRI的计算
Snipaste_2022-05-20_21-49-38

结果的解读和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==210# 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 多平台发布文章来源地址https://www.toymoban.com/news/detail-438818.html

到了这里,关于净重新分类指数NRI的计算的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 本文带你了解透彻云计算(前世,今生,未来)

    作者简介:一名在校云计算网络运维学生、每天分享云计算运维的学习经验、和学习笔记。   座右铭:低头赶路,敬事如仪 个人主页:网络豆的主页​​​​​​ 对于云计算,我们将会通过云计算的前世,今生,未来,特点,原理等几个方面进行讲解。 云计算的一个 核心

    2024年01月21日
    浏览(41)
  • “中国法研杯”司法人工智能挑战赛:基于UTC的多标签/层次分类小样本文本应用,Macro F1提升13%+

    相关文章推荐: 本项目主要完成基于UTC的多标签应用,更多部署细节请参考推荐文章。本项目提供了小样本场景下文本多标签分类的解决方案,在 UTC的基础上利用提示学习取得比微调更好的分类效果,充分利用标注信息。 项目背景: 近年来,大量包含了案件事实及其适用法

    2024年02月05日
    浏览(54)
  • 矩阵指数函数的计算及例子

            最近在学习矩阵指数函数相关的知识,遇到了一些比较有意思的计算方法,写出来供大家一同学习。         矩阵指数函数的计算是线性系统学习过程中非常重要的部分,下面我将给出四种常见的计算矩阵指数函的方法: (1)定义法: 对于给定的矩阵,计算的

    2024年02月03日
    浏览(42)
  • 【MATLAB第70期】基于MATLAB的LightGbm(LGBM)梯度增强决策树多输入单输出回归预测及分类预测模型(全网首发)

    (LGBM)是一种基于梯度增强决策树(GBDT)算法。 本次研究三个内容,分别是回归预测,二分类预测和多分类预测 参考链接: lightgbm原理参考链接: 训练过程评价指标metric函数参考链接: lightgbm参数介绍参考链接: lightgbm调参参考链接: 1.数据设置 数据(103个样本,7输入1输出)

    2024年02月10日
    浏览(41)
  • 现代控制理论——矩阵指数函数的计算方法

    矩阵指数函数的计算方法有很多,至少有十几种,本次只介绍4钟; 方法1,根据矩阵指数函数的定义直接计算 将系统矩阵带入定义直接计算,注意我们只取有限项,因此得到的解并不是精确解 方法2将系统矩阵化为对焦标准型或者约当标准型 方法3;拉氏变换 方法四:应用凯莱

    2024年02月06日
    浏览(43)
  • 【MATLAB第70期】基于MATLAB的LightGbm(LGBM)梯度增强决策树多输入单输出回归预测及多分类预测模型(全网首发)

    (LGBM)是一种基于梯度增强决策树(GBDT)算法。 本次研究三个内容,分别是回归预测,二分类预测和多分类预测 参考链接: lightgbm原理参考链接: 训练过程评价指标metric函数参考链接: lightgbm参数介绍参考链接: lightgbm调参参考链接: 1.数据设置 数据(103个样本,7输入1输出)

    2024年01月20日
    浏览(37)
  • L1-012 计算指数(Python实现) 测试点全过

    前言: {color{Blue}前言:} 前言: 本系列题使用的是“PTA中的团体程序设计天梯赛——练习集”的题库,难度有L1、L2、L3三个等级,分别对应团体程序设计天梯赛的三个难度,如有需要可以直接查看对应专栏。发布个人的刷题笔记的同时,也是希望可以帮助到有需要的人,我

    2024年02月10日
    浏览(38)
  • 利用GEE计算遥感生态指数(WBEI)——Landsat 8为例

    基于GEE平台,实现顾及水效益的生态环境质量评价方法。 《Water Benefit-based Ecological Index for Urban Ecological Environment Quality Assessments》 论文链接:https://ieeexplore.ieee.org/document/9492814 运行结果: 代码如下(示例): 代码如下(示例): 代码如下(示例): 代码如下(示例): 代

    2024年02月12日
    浏览(37)
  • Python空间分析| 01 利用Python计算全局莫兰指数(Global Moran‘s I)

    空间自相关(spatial autocorrelation)是指一些变量在同一个分布区内的观测数据之间潜在的相互依赖性。Tobler(1970)曾指出“地理学第一定律:任何东西与别的东西之间都是相关的,但近处的东西比远处的东西相关性更强” 全局莫兰指数(Global Moran’s I)是最常用的空间自相关指

    2024年01月21日
    浏览(37)
  • 金融时间序列分析:Python基于garch模型预测上证指数波动率、计算var和var穿透率、双尾检验

    目录 一、收益率波动效应的分析 1.1  收益率序列平稳性检验 1.2 建立AR(p)模型 1.3 Ljung-Box混成检验残差序列的相关性,判断是否有ARCH效应 1.4 建立ARCH模型  二、GARCH模型与波动率预测 2.1 建立GARCH模型 2.2 波动率预测 三、正态分布的假设下通过波动率计算VaR  四、厚尾分布的假

    2024年02月04日
    浏览(53)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包