- 在线时间
- 1630 小时
- 最后登录
- 2024-1-29
- 注册时间
- 2017-5-16
- 听众数
- 82
- 收听数
- 1
- 能力
- 120 分
- 体力
- 564676 点
- 威望
- 12 点
- 阅读权限
- 255
- 积分
- 174626
- 相册
- 1
- 日志
- 0
- 记录
- 0
- 帖子
- 5313
- 主题
- 5273
- 精华
- 3
- 分享
- 0
- 好友
- 163
TA的每日心情 | 开心 2021-8-11 17:59 |
|---|
签到天数: 17 天 [LV.4]偶尔看看III 网络挑战赛参赛者 网络挑战赛参赛者 - 自我介绍
- 本人女,毕业于内蒙古科技大学,担任文职专业,毕业专业英语。
 群组: 2018美赛大象算法课程 群组: 2018美赛护航培训课程 群组: 2019年 数学中国站长建 群组: 2019年数据分析师课程 群组: 2018年大象老师国赛优 |
; X0 i5 A' ^/ @( |# c: E
净重新分类指数NRI的计算
& i1 J' m! c0 K: G$ p1 l9 `, ]“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。$ O) s6 F% m. b4 Z/ A& B
NRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!
7 n6 R$ w8 j, ^2 Y9 _" F8 J
: U) r! `$ I/ h# m1 V6 A在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens包,PredictABEL可以计算logistic模型的净重分类指数,survNRI可以计算cox模型的净重分类指数。/ J9 T" p) h# l: h
8 [, @8 d; |, A# slogistic的NRI
- J7 |9 l% ?1 knricens包; z: i+ K/ f' [
PredictABEL包
* R1 ^) |% B; d! n' O8 u生存分析的NRI7 t& o* z5 v1 z/ V
nricens包+ n! ~# W( l W+ c0 u
survNRI包
/ a9 h! b6 |+ {8 A. |$ Mlogistic的NRI9 y6 I8 g. n' E; A7 `
nricens包
( u4 |+ P" X9 r9 @* Y0 i* e#install.packages("nricens") # 安装R包/ l6 o. h; b! X _# }* k# v$ ^) Y2 s
library(nricens)& f M4 t2 ]2 e1 w& h: a
1
1 M! S! {8 [+ c## Loading required package: survival; R$ N& ^! N( Z1 V& S7 V/ i
1
6 U, K/ G9 S6 Y& U# M使用survival包中的pbc数据集用于演示,这是一份关于原发性硬化性胆管炎的数据,其实是一份用于生存分析的数据,是有时间变量的,但是这里我们用于演示logistic回归,只要不使用time这一列就可以了。. ?8 t: h: J: P$ g3 B
% u) R+ i* I/ c5 f$ Rlibrary(survival)8 F, G6 ^; E+ t& t" P) f9 ^, w
2 S. ?" ]# n, D( ^
# 只使用部分数据! a' p' }, u* D3 G
dat = pbc[1:312,] 0 [$ y' X2 b( E6 M4 E/ q5 ?
dat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]
$ F3 t. h& [' A$ d5 c3 [2 w& U2 g
str(dat) # 数据长这样2 }% V7 |3 V) g" G
12 n4 C! K, i8 p8 T0 b
## 'data.frame': 232 obs. of 20 variables:
1 [4 V/ k3 o8 T' G$ M( H## $ id : int 1 2 3 4 6 8 9 10 11 12 ...
9 u" t# J5 a/ {+ @## $ time : int 400 4500 1012 1925 2503 2466 2400 51 3762 304 ...
3 m& |/ p+ b9 A" |/ D## $ status : int 2 0 2 2 2 2 2 2 2 2 ...
5 n6 J! h) ?( f! }3 w9 F## $ trt : int 1 1 1 1 2 2 1 2 2 2 ...0 D2 h' s( Y3 r% `3 Z! I2 a
## $ age : num 58.8 56.4 70.1 54.7 66.3 ...2 J4 t9 ?8 |$ Z! B# @, S& ?
## $ sex : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
) o0 h& j+ D& K. V2 J## $ ascites : int 1 0 0 0 0 0 0 1 0 0 ...
0 G5 v2 I0 }% u## $ hepato : int 1 1 0 1 1 0 0 0 1 0 ...
/ k; a9 U k0 x, r5 v9 A## $ spiders : int 1 1 0 1 0 0 1 1 1 1 ...
1 R4 b7 K8 o" K6 c+ X0 k( T## $ edema : num 1 0 0.5 0.5 0 0 0 1 0 0 ...) p# V; P2 ?( }; z
## $ bili : num 14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ... W- R* W) ]; O9 M2 E
## $ chol : int 261 302 176 244 248 280 562 200 259 236 ...3 m2 i, q2 O( ~* S' x+ y: g8 s
## $ albumin : num 2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...
( h& q B5 V. @ `# k, A V## $ copper : int 156 54 210 64 50 52 79 140 46 94 ...1 U- H/ }8 _. Q3 |& d, y) } T9 Q
## $ alk.phos: num 1718 7395 516 6122 944 ... }. F8 o3 a1 u' F
## $ ast : num 137.9 113.5 96.1 60.6 93 ...! p% X/ R# A, i N9 g
## $ trig : int 172 88 55 92 63 189 88 143 79 95 ...3 h7 X5 n& F5 o- k! A! p4 Y
## $ platelet: int 190 221 151 183 NA 373 251 302 258 71 ...
5 w" T: @" m% {+ L## $ protime : num 12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...% m* A& s9 n/ d6 Q1 l1 ?) [( j1 `5 a
## $ stage : int 4 3 4 4 3 3 2 4 4 4 ... R" r) H' v7 {4 Q: o6 i0 j& ?
9 ^0 o1 v( G \6 c7 d; q' ` [ d1
4 s: Z' L9 ~% {, v9 ^$ O: x: m. Mdim(dat) # 232 20) t, _' C5 ~9 @3 ?0 L
1* m% ^( R x: x
## [1] 232 20
' q/ y" [2 W3 I0 Q" I1
% F' g) T, `9 ?然后就是准备计算NRI所需要的各个参数。+ z6 x; J' \% s$ t( |/ Q
) z& I/ ?' D( P# {* ~
# 定义结局事件,0是存活,1是死亡
/ L" V- e1 }7 L9 revent = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)8 _4 q9 R1 q6 m7 I, S
4 y; A" e" p- ]& }# 两个只由预测变量组成的矩阵4 {& I. m5 e, s- h5 w9 P
z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
7 s' F# o' J: Mz.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))
! R. P& J- h" B5 r4 i. d: M) ?* J2 g0 b: C0 i
# 建立2个模型: D, w3 ~# V4 u& ]. B
mstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)
6 L$ j* O: H Smnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)2 r! Y, u4 x8 {9 c9 ? R0 o; E
" y @8 p: c+ f2 j3 e# 取出模型预测概率
/ o4 R7 p* C% p* H2 F6 Y3 S. Yp.std = mstd$fitted.values
3 |! c6 [, v3 C( Fp.new = mnew$fitted.values
% q0 z' P6 _9 ^; H" r. M! t. g# |# Z" X# }$ @; i' ]- H
1
8 d1 v+ c$ j% ]: Q2 W3 k+ s4 Q/ \: P4 W然后就是计算NRI,对于二分类变量,使用nribin()函数,这个函数提供了3种参数使用组合,任选一种都可以计算出来(结果一样),以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。
7 ?5 i+ _" v% H5 z7 b0 e
+ J2 K- [$ @* A( c9 n' |# 这3种方法算出来都是一样的结果" A' @4 f+ r! j4 a) l
$ n, ]( v- ?2 `8 A7 ^2 A# L+ d: k N- ]# 两个模型
. h6 h2 l; G$ U4 ]nribin(mdl.std = mstd, mdl.new = mnew, " D8 L$ N, v V; _. ?* B, M' y( h
cut = c(0.3,0.7),
- b+ s; b) s0 ?# x! }! ?! _- L3 ~ niter = 500,
- V" z8 c/ |' V updown = 'category')' R T% e v0 J- k
( K8 q7 D' h8 \( q) _# 结果变量 + 两个只有预测变量的矩阵
4 ?( Q/ I" z: u9 `% ~; Ynribin(event = event, z.std = z.std, z.new = z.new, . q8 v. ?/ Z D3 B! X! _
cut = c(0.3,0.7), 2 D2 O1 h9 I- h& b p Y
niter = 500, . S! o$ B, i H; x: `
updown = 'category')
6 k+ D9 {7 t8 i/ k: w) H. {, E/ d' z: o# d7 v- E/ \8 n
## 结果变量 + 两个模型得到的预测概率
4 f; N* `( `' h; o& `9 Znribin(event = event, p.std = p.std, p.new = p.new, * S0 C6 P0 e9 s, y% t
cut = c(0.3,0.7),
0 x; |# I; _" }! v$ _# ?9 V niter = 500, + _8 V4 }8 @, [ Q" C6 N
updown = 'category') o3 L# K6 X# ~
% O- ~3 s# E+ j/ w' m
1
9 G, _# C4 V1 K M其中,cut是判断风险高低的阈值,我们使用了0.3,0.7,代表0-0.3是低风险,0.3-0.7是中风险,0.7-1是高风险,这个阈值是自己设置的,大家根据经验或者文献设置即可。
. b' ^; j6 l+ ]# L2 {) y8 _+ o# d# A! T; a
niter是使用bootstrap法进行重抽样的次数,默认是1000,大家可以自己设置。. f% V) {- P: C! j$ b; F6 D
9 ]! I9 T0 n7 a3 e# a Tupdown参数,当设置为category时,表示低、中、高风险这种方式;当设置为diff时,此时cut的取值只能设置1个,比如设置0.2,即表示当新模型预测的风险和旧模型相差20%时,认为是重新分类。. m) ]2 X: ?- d# v3 u( F3 {
6 S3 L0 O u- H4 v3 Q" ?
上面的代码运行后结果是这样的:
N( t) R$ L1 B# D& W% I: B& r" N/ \5 x
UP and DOWN calculation:
9 y" U5 R0 A0 X. m" j( i% W; V3 w #of total, case, and control subjects at t0: 232 88 144: D; T$ C) Q; Q$ x* G' x! p
# h0 D+ f8 p N5 ^
Reclassification Table for all subjects:
; V+ y" p% |$ Z5 D. q6 V$ K New
" N6 {+ i2 e6 Z: U7 zStandard < 0.3 < 0.7 >= 0.7# y- O! K* [( \ u$ V3 q
< 0.3 135 4 07 U# \6 K1 o; F* ?' g P
< 0.7 1 31 4
3 Z0 I2 @9 k3 m+ G/ f >= 0.7 0 2 55
$ c+ I3 F9 M7 k& F; j4 e% a
+ g5 W9 U$ X# z1 F6 |) Q Reclassification Table for case:+ {& o2 c- ~6 `% u
New4 S% M$ H! {0 F' z' A K: M
Standard < 0.3 < 0.7 >= 0.7
; g' X+ W7 p0 g7 P# Y6 B" ] < 0.3 14 0 0
2 [3 `( ^3 A5 A9 E5 t5 P3 @ S < 0.7 0 18 3
! e! K- O8 W& ~5 b& @ >= 0.7 0 1 52$ \' g1 w- \- {/ `7 h0 v
! D5 e2 s$ F8 i% l" q6 q3 ~% @ Reclassification Table for control:" X3 w7 y3 L& X& o% a$ P; ?1 c
New6 ?/ P2 U$ F+ p8 y: s1 h
Standard < 0.3 < 0.7 >= 0.7
( B# H2 p% F. j < 0.3 121 4 0; h' V6 x: z" L0 Y: K
< 0.7 1 13 16 n# K2 s% u/ R' {8 |& C0 R; Y
>= 0.7 0 1 3/ T% F: M2 K" P# d4 Q/ p ?
" V* ]5 l9 w2 M. O9 Y' J
NRI estimation:# L: {. D4 f+ i+ X8 A0 |9 [0 q
Point estimates:
+ I4 k0 J4 t) H; ~: Z Estimate8 v' o8 p5 Z: p
NRI 0.001893939, t) h( Z1 e( x. Q9 R
NRI+ 0.022727273! k* l( L& O6 T+ j/ z* k" D
NRI- -0.020833333
. @1 ?+ g, ?, L# uPr(Up|Case) 0.034090909: c4 X6 z; r! e# Q) K" Z4 w6 R
Pr(Down|Case) 0.0113636361 y+ V% H) O. L# O8 \- ]. ]" A
Pr(Down|Ctrl) 0.013888889
3 l. W. M( E0 u$ A( H4 XPr(Up|Ctrl) 0.034722222, M, S E$ s3 h5 C6 M L
: ?8 `. F" C4 L f( n ]Now in bootstrap..( A9 ^% n) k7 ?: h: U
7 N2 D2 z4 J. Z, |4 X: gPoint & Interval estimates:3 M8 E, @. s# x3 s X3 Z0 Y- Y3 a
Estimate Std.Error Lower Upper8 u6 j9 _9 _4 J+ Z1 f: w) @
NRI 0.001893939 0.027816095 -0.053995513 0.055354449' ?3 F4 Z( {1 a( H
NRI+ 0.022727273 0.021564394 -0.019801980 0.065789474
7 R' R7 e( I* a/ q3 Y6 `NRI- -0.020833333 0.017312438 -0.058823529 0.007518797
! u/ F# \+ Q& H6 k4 i! |; H6 fPr(Up|Case) 0.034090909 0.019007629 0.000000000 0.072164948* J3 O8 n) N1 a' o, G
Pr(Down|Case) 0.011363636 0.010924271 0.000000000 0.039603960
+ G6 d6 O3 D7 \& N$ q r' E/ \Pr(Down|Ctrl) 0.013888889 0.009334685 0.000000000 0.035211268
* _# V+ u6 R& r. V' V" MPr(Up|Ctrl) 0.034722222 0.014716046 0.006993007 0.066176471
2 J0 O' k% K% A7 ?- F4 l' H- {7 ^9 Z; F& `
14 J! m L3 M& ? Y" a7 M) J) U# t, g
首先是3个混淆矩阵,第一个是全体的,第2个是case(结局为1)组的,第3个是control(结局为2)组的,有了这3个矩阵,我们可以自己计算净重分类指数。
5 A0 i% Q: M, x6 E5 R( D. c2 |
4 F/ T' y; A4 E看case组:
3 l- A8 T5 y# m4 {
- d2 o" e# u( [/ q; ~' @5 Z7 d6 I净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.0227272735 r2 X3 S! J5 q; L3 B9 B+ Q8 R0 X5 [
8 l' r4 b. C% v$ M再看control组:7 Y* c4 A. r2 L! w. g# y& \4 v
0 p& c6 {" _& o# ?$ C/ i净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333. D( @* {* A& n; v" ]# B9 s9 Q
% L4 G+ y; [, ~+ r7 } o相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.000315657$ h; |% k, z Q' u" b
6 E( J( B9 y6 `% o
再往下是不做bootstrap时得到的估计值,其中NRI就是绝对净重分类指数,NRI+是case组的净重分类指数,NRI-是control组的净重分类指数(和我们计算的一样哦),最后是做了500次bootstrap后得到的估计值,并且有标准误和可信区间。
9 C! [4 a2 L) m! I0 X& c0 Z% K- i# Y; p: `
最后还会得到一张图:5 L) V0 r/ W" F
5 A$ H- ` \6 W K* T/ ?' \& l
这张图中的虚线对应的坐标,就是我们在cut中设置的阈值,这张图对应的是上面结果中的第一个混淆矩阵,反应的是总体的情况,case是结果为1的组,也就是发生结局的组,control是结果为0的组,也就是未发生结局的组。
8 C% }; p% e! z1 X
3 T s' M; E0 i8 o7 H# nP值没有直接给出,但是可以自己计算。
8 \# Y$ \& P+ K) k8 [ h
+ ~# j; x6 L; L- [# 计算P值
7 Z8 N, {# x6 M; \ z1 s1 ]z <- abs(0.001893939/0.027816095)
; y# ?& b" p) ~7 y: i) ap <- (1 - pnorm(z))*28 _$ s7 K# K- Q. H1 O" L
p& z9 r# g- s% P
1& @( P- n4 ~4 C+ F( L- D8 O R
## [1] 0.94571570 e, O M& c+ ]: N. M* f
16 ?/ v: t" y! F4 K
PredictABEL包
& u I3 X$ [& H% o2 p6 D! ^' y; c$ {#install.packages("PredictABEL") #安装R包
/ O/ J0 b5 i2 k4 v2 u% W6 hlibrary(PredictABEL) " H& o3 }5 G1 ^0 u3 C2 @
& r9 P6 y* _+ q- x5 m! _+ `3 V+ L8 E
# 取出模型预测概率,这个包只能用预测概率计算
. d6 p8 y; K: L+ Y9 G6 xp.std = mstd$fitted.values
" _- P& b8 \2 i+ r" gp.new = mnew$fitted.values , Y7 j% D9 h* d1 [
1* C( s, c) i3 o% ]
然后就是计算NRI:
4 Y4 U3 s. K$ Y, n2 {, s$ D2 [# ?! Z; q7 h6 Y0 U/ \# p* F3 p$ N
dat$event <- event
% _7 S. M' w) y9 j( M' v- c" e% q& O7 f
reclassification(data = dat,
/ u+ ?0 q) n" y cOutcome = 21, # 结果变量在哪一列
^, Q X5 x [7 o- H* B predrisk1 = p.std,
+ n, Z# X& |6 V& W: `, E& ^9 R" R predrisk2 = p.new,
2 K, ~* Z. @/ N& K4 ^- k. t cutoff = c(0,0.3,0.7,1)
+ V: m% `; \5 Y0 Y+ Q )9 m3 U- a9 o2 j# p" Q" C
1! M5 Z" a# N0 e {
## _________________________________________: F: D% P/ z! M3 C
##
1 l# C9 ]- v* Y## Reclassification table 9 p- B7 }) t, k( a- d; d/ k
## _________________________________________% y8 i6 Z C& Z
## & ^9 y* n" Y$ d6 Y/ c; a
## Outcome: absent
! v7 k4 @1 m5 B: O6 V, W## 8 R; ^: z h* @1 x6 Q. W
## Updated Model
1 B# y2 Q' W/ u0 R o0 g## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
4 ~ p4 n L; A+ u## [0,0.3) 121 4 0 3
0 T- A. U' I+ v4 M/ W8 Z, z8 I## [0.3,0.7) 1 13 1 13+ Y) S+ {. h5 V8 U$ G1 Z0 b1 W
## [0.7,1] 0 1 3 25( @; A. U4 U! K
##
: u. U8 z( q& w( o## 4 n6 C7 Z1 M6 e" q: e' V- V$ r) \
## Outcome: present
6 P% F/ c. P' E) s, H## $ E2 k9 [& H, o4 h$ b v# {& N
## Updated Model- c) P# @8 }) a7 m# V6 R/ Z
## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
, Z. M' C6 F6 |: J" h; y) g## [0,0.3) 14 0 0 0# p- V8 n# x8 A% m ?( b# F3 {
## [0.3,0.7) 0 18 3 14
0 O) a) h- @6 i: j; v( ?## [0.7,1] 0 1 52 2: o- n- ?' E v* K5 P0 t
##
, P4 l6 S3 p; n$ V, Y( F7 Q##
1 |' U- S3 \# G7 i1 t, }; M## Combined Data
) Q% H' Z. u( |2 i: X2 |## . L$ M" f, m) E4 M9 y0 m
## Updated Model6 y. G7 X) G3 o/ ]+ s% O
## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified5 x$ B$ {: Z( k4 A
## [0,0.3) 135 4 0 3
$ ^7 [4 H$ l8 S2 G$ X! P: N) r: @## [0.3,0.7) 1 31 4 14
* Q K( J, _, }: L0 m+ Q9 C## [0.7,1] 0 2 55 46 Z5 e+ o! t# O3 Z" I H2 q
## _________________________________________
2 M7 H: R5 r7 [6 r. T## 3 u4 c: [" e' H( @0 j, E
## NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806 . ^, c4 t% t8 N3 z9 _# |
## NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048 & ^6 i* A) J5 z% T5 z1 c/ g, `" Z
## IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.283966 m \* n9 r8 h5 l! `% v
9 B2 e1 T# L" H! v
17 v" k% p1 S2 x( x
结果得到的是相加净重分类指数,还给出了IDI和P值。两个包算是各有优劣吧,大家可以自由选择。
/ d5 O+ I% I3 O4 \/ M" }
4 S5 D8 l2 Y- x5 b8 Z9 s _生存分析的NRI
3 P# j$ W [7 n, J还是使用survival包中的pbc数据集用于演示,这次要构建cox回归模型,因此我们要使用time这一列了。
2 S) y! C& E5 p" n
( j4 H0 Z3 e, }8 ^nricens包
* Z9 i4 K# r m! z9 c$ Glibrary(nricens)
# s# C' G9 ?5 R( olibrary(survival)
4 p# X+ _0 n( g) C, g: ~6 z5 _0 T& [# {9 d$ g; [6 w" o
dat <- pbc[1:312,]
' ?0 J* i% z% e: J# Bdat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡- u) J2 s R6 E2 b! b( `
1
2 w- [$ I3 X+ n6 D( g7 h1 |' w然后准备所需参数:
! _1 T% Y* }+ G# K
3 M$ J1 r5 u2 S; }' t' O9 m H5 f# 两个只由预测变量组成的矩阵
3 T, C" w m6 I' S6 e" J3 Jz.std = as.matrix(subset(dat, select = c(age, bili, albumin))). k- n2 Y5 ?1 _; V8 q2 g% G" T' m
z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))
a' F. X' o% u# E' ^! O# _
7 g3 h/ d; ?* n7 P# 建立2个cox模型
5 |6 K' o( Q* r2 m2 ?; b% amstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE)# n8 \; g* g# x$ P& \
mnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE)
( w1 t" a1 c4 n1 e3 ]% j2 {" m& A0 A8 C1 o* U
# 计算在2000天的模型预测概率,这一步不要也行,看你使用哪些参数
9 I/ \0 q* z8 w9 k( j& R9 K* O3 Mp.std <- get.risk.coxph(mstd, t0=2000)9 X' X7 w" b; m( X
p.new <- get.risk.coxph(mnew, t0=2000)
1 @/ P( N u) `7 {, Z1
; @9 Z6 d- h p$ u计算NRI:
, k. A2 n* V( s3 _" C! F9 |) k; I: W
nricens(mdl.std= mstd, mdl.new = mnew, 3 S. u, c* l# ^. S% N! W) z1 `
t0 = 2000,
& ~7 G5 P7 D) J7 z3 D cut = c(0.3, 0.7),1 N6 d6 |( @; m7 x x W; z
niter = 1000,
) H! M1 `' O; n2 }/ H updown = 'category'), g, P: l2 Z5 D2 c# ]' ~) P f' x& h
! W$ m* d+ \7 H1 [2 HUP and DOWN calculation:& }' i! I0 {9 E A2 j
#of total, case, and control subjects at t0: 312 88 144% d: r' z) h0 _- \, V
; t. H& P" q2 C4 H* s4 j Reclassification Table for all subjects:
& |2 Y+ f! @3 X1 C New
2 Q) I+ y# |, yStandard < 0.3 < 0.7 >= 0.7 j- K) j* w. L1 L5 K$ F2 i0 L1 p
< 0.3 202 7 0. m' c/ }& N8 N; Q
< 0.7 13 53 6
3 c' m& j, `$ H, O; Q >= 0.7 0 0 317 H# u0 g, y0 U) H% o
% Y8 Y1 L9 T8 N7 X' F
Reclassification Table for case:+ }' @+ Q: Y4 e+ p; t
New. E/ a. P/ V1 {0 m4 l% f: p- i
Standard < 0.3 < 0.7 >= 0.7
" F! p1 G2 G1 U8 d z < 0.3 19 3 0
( P* t7 c& C. ^ < 0.7 3 32 4
^; d) L- u$ [6 u >= 0.7 0 0 272 L* C; n6 U, n& d
" I$ n+ Z2 x2 m% W; ~
Reclassification Table for control:
# u3 D& C! @; q. l% z( P New: F I! e9 A0 Z9 h
Standard < 0.3 < 0.7 >= 0.76 h" P+ R( u' N) B( z
< 0.3 126 3 0
5 |, u' ^ ^* s) [( V8 X < 0.7 5 7 2: b7 Y4 G4 ?0 z6 z" w5 P
>= 0.7 0 0 1
+ K; O0 J7 N5 B3 V* e, T' k5 n# s2 J Y( c7 ?- `0 U
NRI estimation by KM estimator:
. X. v8 A% B3 A N' k; N0 z) B! t% _/ k
Point estimates:
! M/ L" o9 s6 N" }/ _9 I% M6 Y1 l1 N Estimate6 Y" d2 q9 w6 a4 ]9 A
NRI 0.053776358 n7 f$ H7 [6 Y7 {6 d; |
NRI+ 0.03748660
, T, f' w& G" h1 ONRI- 0.016289742 x' h3 z- z1 [ z- r# G
Pr(Up|Case) 0.07708938
# k2 x+ f6 t$ J K/ `Pr(Down|Case) 0.03960278; U/ n/ j; y' a) c7 M- u* w
Pr(Down|Ctrl) 0.04256352/ ?4 Z, q/ k1 g6 k! P8 n4 V
Pr(Up|Ctrl) 0.02627378
6 f4 q8 K+ h y% u/ \: O1 A! B2 R+ E$ X( N: R8 o6 G- t
Now in bootstrap..( M- G# b5 R7 V ?% m# ]3 j- H
# D8 v. A# w3 J/ EPoint & Interval estimates:
O" P* q! N+ _$ ~! b Estimate Lower Upper# n; y3 b" C0 C
NRI 0.05377635 -0.082230381 0.16058172
i" V; W5 e% u3 ?NRI+ 0.03748660 -0.084245197 0.13231776* I, g) A. L+ W& ]( Q- s& y
NRI- 0.01628974 -0.030861213 0.06753616. L7 p! ?+ V4 Q- m1 c% A' M
Pr(Up|Case) 0.07708938 0.000000000 0.19102291
, V# k; X9 M4 oPr(Down|Case) 0.03960278 0.000000000 0.15236016
+ C. O5 x5 W" \- S/ n& n4 sPr(Down|Ctrl) 0.04256352 0.004671535 0.09863170% s& I3 p0 n4 Z
Pr(Up|Ctrl) 0.02627378 0.006400463 0.05998424
" z9 d4 W D. i; W( U/ _3 g" U
0 M0 E# t9 i( f U1
' Z9 |4 F1 M: G" Z
4 j$ M4 A5 N l( Z8 p% }" Q: M" r* zSnipaste_2022-05-20_21-49-38
5 d: G' ]8 v, K结果的解读和logistic的一模一样。) l% C8 H/ W0 d$ l
3 J* N% V6 D) _4 s; b+ k/ usurvNRI包
4 ~+ x; w4 g! u# 安装R包1 f) {, I0 I5 y# t4 q, o1 a
devtools::install_github("mdbrown/survNRI")8 [/ Z3 m4 n0 Y' A6 P
11 W- C7 f: F2 ^0 g2 v
加载R包并使用,还是用上面的pbc数据集。
: f& F7 o& j) \( r( j w& F
9 P9 p1 v* M- vlibrary(survNRI)
! ]) S% N* {! ~2 W3 q' {" t* C/ W# |4 }1
/ y* r! E& R8 I% E5 y## Loading required package: MASS- O! X( f1 M& t, Y, P
1
F1 H- N6 @7 s6 S0 n3 Alibrary(survival)
2 n% J$ ` N9 y7 R# c# a' k. H
* I" e0 I" E1 r" Z! B9 w# 使用部分数据
7 x1 F6 v; u* p: o1 p, k+ hdat <- pbc[1:312,]
& ?7 }" y, @3 k- ^2 _. k: Cdat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡$ V4 g2 \; X! ~% ]# q
4 e/ O7 r; H Y- k7 o$ Y3 ?
res <- survNRI(time = "time", event = "status", # ?; l g; n! H0 Z
model1 = c("age", "bili", "albumin"), # 模型1的自变量
3 H/ [; w+ O/ r2 N8 g model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量
5 j" F+ H7 s4 O9 Z1 F/ K data = dat,
! B- W- J4 X8 T9 h0 z# ? predict.time = 2000, # 预测的时间点
, v9 b4 v) ]/ x S method = "all", 5 Z, F9 A1 v. U: C" {" v/ {' y
bootMethod = "normal", 6 n0 U7 m4 b, e, a* l$ x
bootstraps = 500, ! z& R# L& [* O9 ?, o) ~- x! i/ \
alpha = .05)
* j, o5 `! z* ~& D
& a! O- v9 F c5 |14 |4 L7 _ { \! W) A6 A
查看结果,$estimates给出了不同组的NRI以及总的NRI,包括了使用不同方法(KM/IPW/SmoothIPW/SEM/Combined)得到的结果;$CI给出了可信区间。% n( T! H, q9 h. n6 M
' k6 {- X. g# ?) v% A" }: V( Nres0 `9 C2 A. S" v
1
5 W' h6 h9 X7 w t) A## $estimates& k; u9 k+ t) X }( p+ m
## NRI.event NRI.nonevent NRI
' ~2 G- X; `. P## KM 0.20445422 0.3187408 0.5231951
1 C5 O8 n5 p5 J0 x$ i3 ~' @3 B## IPW 0.22424434 0.3273544 0.55159872 V6 O- K8 Q! w% O# \; I' o
## SmoothIPW 0.19645006 0.3144263 0.5108763
" [, B1 x h. @- c; v$ k5 B## SEM 0.07478611 0.2632127 0.33799881 o/ @' o- e2 [; W7 I, } |" I3 E
## Combined 0.19633867 0.3143794 0.5107181% x1 z+ S: j/ j) ? G/ w9 s
## * A( O+ |% E& r8 C: X5 Y3 M* w0 F
## $CI
7 n5 j' O0 L3 p3 G## $CI$NRI.event+ u6 l3 a( ~5 {6 i; a
## KM IPW SmoothIPW SEM Combined
3 b: w: a. O- e4 Z& k" f## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723( D5 r* P/ |1 b* X
## upperbound 0.44806768 0.47033936 0.44014214 0.2658309 0.4400496. F" s; B7 K6 d% r: X. l
##
* g: K$ ~/ w: \& @' q" E0 ^## $CI$NRI.nonevent
& J K; @ |# w6 r## KM IPW SmoothIPW SEM Combined
' i3 \7 H/ C5 @6 b% l6 D- A3 K## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.1286426
0 R) y8 h P7 _. S4 S4 Z## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.69645492 p/ n$ I% `$ g" D5 p
## : O6 B# J6 B, ]. H( V z( d; }8 O
## $CI$NRI
/ a/ Q! z; a# U' U## KM IPW SmoothIPW SEM Combined% k$ e; y! H2 @2 r1 t, r
## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409
3 v) j' x( V% ^! ^## upperbound 0.89306122 0.92464359 0.87970125 0.64253510 0.87953153
5 p, }8 M- J$ {5 }7 W##
0 l) F' P4 F1 E( J2 O% Q6 O##
& Z* b$ l( L: j' r) {' D## $bootMethod
_0 r( r/ `0 i4 L- G' _## [1] "normal"
/ m0 t# A. E& V7 ~## ) F. q7 r, e* h& C; T
## $predict.time' j: v8 O' v4 u9 k& a
## [1] 2000( L% _! q$ y( r9 e% U5 T2 k
##
2 ]2 B4 g( V2 l( ^9 b. b6 b# i" R## $alpha1 p0 t. n( o- d/ `# k
## [1] 0.05- [- M. X( W; y/ Y2 [1 K, \ ]% v5 `
##
6 i* v: b; Q& Q& |1 c, L## attr(,"class")
; w+ J H6 |/ E## [1] "survNRI"
7 h( e2 f' Y4 _2 _/ N' }- V0 k7 D9 W3 H' V, ?
18 ?6 C8 S$ a5 q7 R7 R$ u/ k) K6 s
OK,这就是NRI的计算,除此之外,随机森林、决策树、lasso回归、SVM等,这些模型,都是可以计算的NRI的,后面会继续介绍。大家如果有问题欢迎在评论区留言。2 w% i) Q# s# t8 i; N/ Y. _
; d4 j8 ?+ ?- u2 L2 [. I本文首发于公众号:医学和生信笔记
- O+ f% v9 N" w/ I0 X" U, |+ G3 z- B' J2 g" J- j% W3 g/ a
“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
1 M7 i; G& A# F1 J) W- s本文由 mdnice 多平台发布
; M, p6 v& v9 ?; b————————————————
2 E, S% d; H N( R+ d) S! C版权声明:本文为CSDN博主「阿越就是我」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
: J5 B1 _" e% p/ }( I原文链接:https://blog.csdn.net/Ayue0616/article/details/126768006
: n4 C# ?. a! g) q% e
0 u7 `' {, r6 D1 q# E8 k* F2 M( H
$ q0 u% K) ~* M2 Z, D$ H( W |
zan
|