NRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!
在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens
包,PredictABEL
可以计算logistic模型的净重分类指数,survNRI
可以计算cox模型的净重分类指数。
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.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后得到的估计值,并且有标准误和可信区间。
最后还会得到一张图:
这张图中的虚线对应的坐标,就是我们在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==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: 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
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==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的,后面会继续介绍。大家如果有问题欢迎在评论区留言。