数学建模社区-数学中国
标题:
净重新分类指数NRI的计算
[打印本页]
作者:
杨利霞
时间:
2022-9-12 18:43
标题:
净重新分类指数NRI的计算
+ g4 e* a* a/ _8 ~* \: B7 w
净重新分类指数NRI的计算
" \% W" r3 U( N
“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
1 ]5 A# w' y1 u: D) i4 }1 }4 g, z" D
NRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!
# h5 q: M8 K: G# k! c
- s5 X% H! |! P( c/ \0 ~
在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens包,PredictABEL可以计算logistic模型的净重分类指数,survNRI可以计算cox模型的净重分类指数。
$ e) s: s8 w: ?( H* \ h
2 x& S% a9 B" D0 V8 O& q" m
logistic的NRI
+ c1 J$ l# |9 S2 x
nricens包
& f$ I7 |3 O& J# C& k' k( ~
PredictABEL包
# |# A$ r; B5 y0 t/ d$ l: k
生存分析的NRI
8 Q, G' j) Z' A+ |$ a
nricens包
) W. X' r1 P0 V2 w; w( T: ^
survNRI包
( L+ ?+ k: L2 ?+ b- L/ Q
logistic的NRI
( B2 C& p! j* l( j" c" C
nricens包
* f, T# ]$ [+ d' ~/ V6 ~
#install.packages("nricens") # 安装R包
2 T6 h |8 I4 Z( d/ p4 c( c
library(nricens)
' Z) A9 l3 l' `. E
1
% B5 y& t; a- d& X& s+ h
## Loading required package: survival
3 F2 v1 g z( O, Z0 C7 J+ {$ P
1
( D* m7 a m. e$ T, W
使用survival包中的pbc数据集用于演示,这是一份关于原发性硬化性胆管炎的数据,其实是一份用于生存分析的数据,是有时间变量的,但是这里我们用于演示logistic回归,只要不使用time这一列就可以了。
4 i9 x k1 R9 b! u8 b; b
* C2 o# U: E. i( {
library(survival)
) `0 t/ q- @8 J4 Y m
2 H: G7 e. V4 s8 F" L
# 只使用部分数据
- {4 c+ d8 s% `% _/ c: d2 H
dat = pbc[1:312,]
# R5 ?: Y; [' q1 s8 J# l
dat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]
0 Q; b4 \; ?, ^% j
$ j* x* F+ y* k, l9 [
str(dat) # 数据长这样
" r/ d7 {% Q* ?. Q; {% J
1
, W7 z# v" f4 o) F# g8 h
## 'data.frame': 232 obs. of 20 variables:
3 a# ^+ H" h, f2 D6 T, _
## $ id : int 1 2 3 4 6 8 9 10 11 12 ...
$ n- d+ j3 k; T% K( @' W
## $ time : int 400 4500 1012 1925 2503 2466 2400 51 3762 304 ...
3 @" u8 Y! N, [) C* f1 ?% i
## $ status : int 2 0 2 2 2 2 2 2 2 2 ...
# R9 l( O& o! C/ u9 T& [
## $ trt : int 1 1 1 1 2 2 1 2 2 2 ...
: a+ T" f# _/ ]2 p$ }. R
## $ age : num 58.8 56.4 70.1 54.7 66.3 ...
* U1 G, @' w6 l3 R: X! g
## $ sex : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
4 t. x `# v) q/ e' T; w! W
## $ ascites : int 1 0 0 0 0 0 0 1 0 0 ...
3 r, V Q1 d% v: W
## $ hepato : int 1 1 0 1 1 0 0 0 1 0 ...
. ^2 k ^) N- p, u+ P- C
## $ spiders : int 1 1 0 1 0 0 1 1 1 1 ...
6 {/ v2 M& w! p" k7 Y; a2 n
## $ edema : num 1 0 0.5 0.5 0 0 0 1 0 0 ...
7 {4 h2 T: T6 c5 r2 J4 e4 [
## $ bili : num 14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...
4 h' a8 E4 |6 }$ m8 s
## $ chol : int 261 302 176 244 248 280 562 200 259 236 ...
+ ^! A) w6 R" v
## $ albumin : num 2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...
1 s! o$ z9 j. N1 ^/ a' x$ I: {
## $ copper : int 156 54 210 64 50 52 79 140 46 94 ...
) z+ o' |0 w/ p
## $ alk.phos: num 1718 7395 516 6122 944 ...
0 l' R x' s+ N& ~5 T
## $ ast : num 137.9 113.5 96.1 60.6 93 ...
9 A; M) \1 o8 ^8 P3 {. c8 _
## $ trig : int 172 88 55 92 63 189 88 143 79 95 ...
& g5 E0 n) X5 E i
## $ platelet: int 190 221 151 183 NA 373 251 302 258 71 ...
3 d7 G6 y- \! M# M% a$ M
## $ protime : num 12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...
& ^# o+ W8 _# \" ~* h# |
## $ stage : int 4 3 4 4 3 3 2 4 4 4 ...
; Q# o! m5 y/ ^/ j
) r$ E3 |3 Y* A1 L0 z: u2 F
1
2 J8 f* J( s3 |# d: x" w0 `9 H
dim(dat) # 232 20
) y% h/ `& U+ _) A+ {! k
1
9 T) n( {0 k1 _& x8 U
## [1] 232 20
8 W2 K( L/ i) N# w
1
8 ~, @. f j7 i& d S5 X, F, b3 _
然后就是准备计算NRI所需要的各个参数。
9 ^( j, `2 _4 ]: j C
( P; Y# y1 [+ J$ j' z
# 定义结局事件,0是存活,1是死亡
. [# t( ]0 {8 ? d
event = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)
8 ~, s) ]* v+ w# B% a6 p3 ?
/ ~3 ^1 E1 v, p& {0 ^1 b1 i7 D
# 两个只由预测变量组成的矩阵
$ h) \# i) F* s/ }
z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
$ i2 f1 b" `2 w0 R2 O
z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))
" q) J4 |9 x7 W. V
, y( ~' n- P5 R9 h
# 建立2个模型
- r9 A2 _" P0 } h1 v
mstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)
9 ?9 s ], \8 `% O3 @% R
mnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)
% C% \" h. b& h1 F+ W' }# Y
! o, Z v, ~8 G7 Y1 s% w
# 取出模型预测概率
' p p' { L+ Q2 I7 t9 C+ Y
p.std = mstd$fitted.values
8 D. [& U! A6 U7 k
p.new = mnew$fitted.values
2 z: j( n% ~; b* q" Q$ ~
% q x% d9 _4 _2 c& t) v
1
7 a; I7 L0 j" ~9 D2 d: H/ a
然后就是计算NRI,对于二分类变量,使用nribin()函数,这个函数提供了3种参数使用组合,任选一种都可以计算出来(结果一样),以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。
" P9 r: e e' b, V9 x( P
' l* p1 T9 |7 o5 {( p3 Z9 W1 `5 c3 I
# 这3种方法算出来都是一样的结果
' _0 x* j, v& {( N9 @3 `* F. e
! J. V( V- e" R6 [% s0 n3 h9 x4 h
# 两个模型
4 Z/ \8 z" l1 M9 v3 F* k( p8 p& E1 y
nribin(mdl.std = mstd, mdl.new = mnew,
- o" Z2 y9 y( A" \- v6 e
cut = c(0.3,0.7),
! t, M# S: |7 T' Q, r
niter = 500,
; o. w/ U6 ]' K3 T
updown = 'category')
3 z; t E9 m- S# `( j1 Q
& r. c* ~2 d; F r% F( S# G
# 结果变量 + 两个只有预测变量的矩阵
6 F6 j+ p- |1 b0 t. w# b; H* c
nribin(event = event, z.std = z.std, z.new = z.new,
& q% z( B$ Z# Y) D" e4 f6 Z
cut = c(0.3,0.7),
! h7 D5 _" f! p+ U( b0 D3 X
niter = 500,
- a# G) ]7 D! z: i6 W! x
updown = 'category')
& s8 Q6 U1 f% _& q9 n1 G" M5 k3 }. C
' T: t/ \9 _2 Z7 V* v4 Y0 q1 t5 J
## 结果变量 + 两个模型得到的预测概率
3 L! J9 ]9 a; C0 {# J
nribin(event = event, p.std = p.std, p.new = p.new,
! C9 F( X6 v2 p/ _ N) ]
cut = c(0.3,0.7),
8 W. @6 P6 `: O9 L( h
niter = 500,
$ j9 F/ Y9 X% ` G, J; `
updown = 'category')
& J' C7 G# L! f' A, ?! W/ m5 _1 b
& n ~, Y% g7 z( G+ P" ^
1
2 ~2 p( ]0 {+ \8 s
其中,cut是判断风险高低的阈值,我们使用了0.3,0.7,代表0-0.3是低风险,0.3-0.7是中风险,0.7-1是高风险,这个阈值是自己设置的,大家根据经验或者文献设置即可。
* K. x. Q. C1 {3 p% t. M8 ~# t* b
; |& f4 ~4 P" ` ]; F& {" u
niter是使用bootstrap法进行重抽样的次数,默认是1000,大家可以自己设置。
: I3 x% }. D+ U. p" W; l( s
& M) d! Z8 _( U9 e+ G9 t4 O
updown参数,当设置为category时,表示低、中、高风险这种方式;当设置为diff时,此时cut的取值只能设置1个,比如设置0.2,即表示当新模型预测的风险和旧模型相差20%时,认为是重新分类。
: k" o2 A) k( i# _
5 w& G, C' X* E& s8 L- ^5 N/ b" Z0 N
上面的代码运行后结果是这样的:
" r' [% l% V' D* S
, W1 b9 q2 ^# Z0 H! t8 k* V7 w. e
UP and DOWN calculation:
1 H% F! @# s9 I- @8 }1 i4 a
#of total, case, and control subjects at t0: 232 88 144
' S0 e; J: H7 \1 v1 `1 R4 M( [: N& S9 E
0 t( m# c9 H7 V" ^
Reclassification Table for all subjects:
5 z7 }; h' z5 @' B) t2 D4 R* L
New
" F& q" `3 Z3 I0 f ^; u, R
Standard < 0.3 < 0.7 >= 0.7
# A+ H5 W, r2 a M5 C, W
< 0.3 135 4 0
- R/ `, G s' ]2 l
< 0.7 1 31 4
* H& k# z! A7 s' i3 p. p. S
>= 0.7 0 2 55
0 L e8 W5 G* g5 n7 g3 l; p' x
" T( [: R2 d8 ?# ]; f7 f' E
Reclassification Table for case:
+ c, h5 E4 {0 k$ V6 y" r% j& }
New
9 O3 J4 x' w7 a5 N
Standard < 0.3 < 0.7 >= 0.7
* Z3 R; t! J; s
< 0.3 14 0 0
, g( B# i( a p# V3 w ?
< 0.7 0 18 3
- h) Z. F( ^- `& j u; Q, Z9 C9 x
>= 0.7 0 1 52
( n7 @7 F$ D" w, |5 n
1 X" ]1 C+ M: x5 w8 w: r' k5 p
Reclassification Table for control:
0 x/ t' M3 B: R0 ^7 |3 z
New
8 }# D* Y/ q/ C; V
Standard < 0.3 < 0.7 >= 0.7
! W6 U3 \ O# N5 J) Y4 r
< 0.3 121 4 0
/ T" g9 n( v, ~# V0 o9 ?- K. O
< 0.7 1 13 1
9 u# K8 T# L( a1 {" j. T% F
>= 0.7 0 1 3
^) w2 E' E, K6 Q. }* i+ G
; l' A7 \: _4 ~+ |- @
NRI estimation:
+ ]1 a4 k! e2 S/ k4 R
Point estimates:
) X/ ^, _( C, p) ?
Estimate
' q# J% S4 c% ]* x. q) @
NRI 0.001893939
- e6 G8 @4 ]- \
NRI+ 0.022727273
$ V$ x d& i3 c. V6 Y0 N$ p
NRI- -0.020833333
8 a8 c$ m/ {( v+ R! d: `0 T: Q
Pr(Up|Case) 0.034090909
Z$ r* u3 l; y2 }
Pr(Down|Case) 0.011363636
9 v4 x& z8 ?( z
Pr(Down|Ctrl) 0.013888889
! G6 {2 q' L6 x( [) i
Pr(Up|Ctrl) 0.034722222
9 } O2 Q% x7 l: j6 O7 X4 {6 m. G
$ C3 ~/ R0 ` f+ i" k( V
Now in bootstrap..
' d% J) c g( @8 }) D
" [& k5 F7 b9 j. s& E& S
Point & Interval estimates:
2 g* I+ W$ R q) d/ \3 Y! G
Estimate Std.Error Lower Upper
; i& h' f2 T. n4 [
NRI 0.001893939 0.027816095 -0.053995513 0.055354449
- G# X0 Q1 g5 r" `' i7 @' ^" {
NRI+ 0.022727273 0.021564394 -0.019801980 0.065789474
0 C ?3 t- B7 \) H% |
NRI- -0.020833333 0.017312438 -0.058823529 0.007518797
1 X/ Z, L7 V7 R5 _
Pr(Up|Case) 0.034090909 0.019007629 0.000000000 0.072164948
* [2 r6 E4 a8 e, l" b; ]
Pr(Down|Case) 0.011363636 0.010924271 0.000000000 0.039603960
1 y0 Y! s k( }9 x( z) Y
Pr(Down|Ctrl) 0.013888889 0.009334685 0.000000000 0.035211268
! y# B" ^; x* E
Pr(Up|Ctrl) 0.034722222 0.014716046 0.006993007 0.066176471
/ U5 r4 m9 N' v% n$ Q4 J r* T
: m' a Q3 u6 T0 ~4 _- y& m% }
1
C2 y6 \) z1 D) j# e
首先是3个混淆矩阵,第一个是全体的,第2个是case(结局为1)组的,第3个是control(结局为2)组的,有了这3个矩阵,我们可以自己计算净重分类指数。
6 d- n) }+ ?/ \9 p% c
6 H1 V* u2 R$ b3 R2 R; ?- g
看case组:
6 d5 k% X" Q4 Z+ X
- Z& c# C( i, e$ L% U: p/ v5 b9 j
净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273
- z4 h! J+ M, E, J0 b
8 l" P- Z6 j/ |/ U
再看control组:
2 n( O; x! J5 g7 \2 T. k- `
5 j' q' \% t: m. M) ~ l
净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333
, x) A8 }- }4 i+ o" D# k, w2 M6 ^. c+ _
" N6 M: l1 V3 j9 n0 F) ~) P- ~
相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.000315657
% J% v p J; K" P
! k1 h+ N$ t% v
再往下是不做bootstrap时得到的估计值,其中NRI就是绝对净重分类指数,NRI+是case组的净重分类指数,NRI-是control组的净重分类指数(和我们计算的一样哦),最后是做了500次bootstrap后得到的估计值,并且有标准误和可信区间。
6 x. C' Y$ t, y- C- K3 w
( K$ K6 \( k) O( f/ k: l/ U& K
最后还会得到一张图:
. z# l( I9 d3 v3 C: g
$ ]& R- E$ ]+ }4 e7 o* p) e! m6 U
这张图中的虚线对应的坐标,就是我们在cut中设置的阈值,这张图对应的是上面结果中的第一个混淆矩阵,反应的是总体的情况,case是结果为1的组,也就是发生结局的组,control是结果为0的组,也就是未发生结局的组。
/ G& W2 C' @2 q8 M) u3 P# h6 p8 ?# n
6 a6 D. \2 d4 g+ T
P值没有直接给出,但是可以自己计算。
* g4 B! o$ i/ ]: J& t
) V: u2 t9 K; K( Q
# 计算P值
7 q7 F3 ^* C( J5 N2 [1 n- e
z <- abs(0.001893939/0.027816095)
/ H7 N) O8 T0 l9 ^) |' w! ?( ?( I
p <- (1 - pnorm(z))*2
! ?! W# O6 F, a K, W8 X* O
p
, W/ V2 J W9 D! K/ W& B
1
$ b, v: U7 J7 r$ \( ?2 {
## [1] 0.9457157
; w( a) D. y# X2 |# F+ N
1
9 \, L( }+ K' F( Y
PredictABEL包
, m* B6 O' y/ u* E3 v
#install.packages("PredictABEL") #安装R包
5 ^9 E0 ]3 g$ o. }
library(PredictABEL)
; r1 V# b3 h1 N/ |" S2 z$ H9 G- e
" Z) T+ U9 P# f
# 取出模型预测概率,这个包只能用预测概率计算
4 {; w; a9 v' \% X2 e" |( Z
p.std = mstd$fitted.values
; r% j/ K: }$ j6 ~4 O8 q5 C8 g3 `! a& O
p.new = mnew$fitted.values
6 z3 E) g! `* Y5 \5 h9 } c% S
1
# d" _2 g: d. f! t
然后就是计算NRI:
) u4 ?' z8 X1 y6 @6 S# |5 p) `
% G. g% f( e6 D! G
dat$event <- event
A/ \( i/ R u. p$ F C" K$ g3 y
' w+ q% l6 Y* ?$ @) `
reclassification(data = dat,
% [3 p" j5 a5 M, M3 M
cOutcome = 21, # 结果变量在哪一列
4 L1 Q+ h* Y8 o, V
predrisk1 = p.std,
- e4 p) z2 S, Z
predrisk2 = p.new,
3 _# q8 u4 w5 v6 t6 C7 G' @
cutoff = c(0,0.3,0.7,1)
) J% x% k2 z5 `; K7 T1 r
)
* {+ c4 c1 A; a9 z2 D$ ^% F, ?9 k
1
) M; h) A) p1 ]0 k! C7 C0 @( p G9 a
## _________________________________________
3 g- J7 ^8 }! x+ k5 b* ]
##
1 w1 |9 ~3 c/ r4 _6 H
## Reclassification table
' z0 D2 _& _8 m9 R4 ^. e$ M. j
## _________________________________________
) A9 {( ]4 f, m+ o$ d( m; {/ l
##
4 d8 A! r( K- [
## Outcome: absent
Z. v: E$ a4 _8 v
##
; }& V% w( E/ J
## Updated Model
0 Q$ C7 o8 F9 E' r: n
## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
" @' w2 ]( B% @" D7 h- q% n9 z$ W
## [0,0.3) 121 4 0 3
( |# O; l. w5 t
## [0.3,0.7) 1 13 1 13
2 c/ J2 x; a. q* A, `' O
## [0.7,1] 0 1 3 25
% Y& m3 ^: c2 c* U o/ t; e
##
2 s4 z# i; j6 U' e5 m! M4 x Q) ~
##
k* Y9 Q8 O2 e: D* q
## Outcome: present
5 G3 L: d3 M' g9 k( J' I
##
. T3 |5 b, A; p* E# r5 k4 Q. f
## Updated Model
1 P# [5 s% K% j" q) w; A( c& m* \
## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
" I/ ~0 n4 _- k& q+ |% a
## [0,0.3) 14 0 0 0
( l& D5 a7 K: ]
## [0.3,0.7) 0 18 3 14
% S- D& Z$ g9 H2 {6 P1 z a( J
## [0.7,1] 0 1 52 2
1 N. m4 v' h |# ^$ L2 x |
##
# a' B; a' w( @* h( B- o
##
( }: V; q6 ]6 t8 }! \
## Combined Data
1 V) S+ U+ ^4 A x- E' z! r
##
' n! L" ]! d3 \, \
## Updated Model
0 Q/ b" t0 Q) ^4 D
## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
: `( r# m% O) ?. L9 e+ S
## [0,0.3) 135 4 0 3
# U. Z1 N. ~; C1 ~7 ?
## [0.3,0.7) 1 31 4 14
: c, G. ~0 e8 b0 l
## [0.7,1] 0 2 55 4
. a. z2 |. G. @/ T" v. c
## _________________________________________
/ k# f. k, s, e& g
##
) z" B# ?! a! `! n
## NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806
& d( e& Q) d/ F B$ P2 K8 w
## NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048
2 J0 `* V: O, Q9 a+ X1 n6 [% Y
## IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.28396
( x- G# M5 P+ _' g) d) g$ |
. r6 d9 W- r" m& b
1
: v: Q0 d" l6 P# h, v
结果得到的是相加净重分类指数,还给出了IDI和P值。两个包算是各有优劣吧,大家可以自由选择。
1 d& O' |, `, Q1 F" }1 c8 b
# z% k# P9 u( N! [ d
生存分析的NRI
( W& j9 e% H8 [! w( [
还是使用survival包中的pbc数据集用于演示,这次要构建cox回归模型,因此我们要使用time这一列了。
# S0 g- ~- X7 ~: L3 ~# x
/ b/ b7 l+ h) v9 z! O
nricens包
: O1 i& D8 i: |5 m
library(nricens)
0 o0 g- f, \3 {' o' G# T
library(survival)
# q! }; V- ~* u* g9 a/ P+ Q, A. z
2 T0 f1 J% ^0 w+ g" y
dat <- pbc[1:312,]
& Q, r5 u* M! o% F- s
dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡
g# S- Q u! ]
1
6 D! H6 I8 z& q4 h
然后准备所需参数:
! f* D4 Q1 q& n
2 f/ h% q$ u3 ^
# 两个只由预测变量组成的矩阵
; P! H4 i7 H" N# g2 M
z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
9 U& ] W. H6 A) J. v
z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))
$ X( N/ X+ u3 {8 v9 X
0 Z7 j4 q/ I0 x$ h- g7 t
# 建立2个cox模型
4 [4 |: M8 j! d; `4 l& N
mstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE)
( b: w' P; A/ z' B) Y6 X+ T
mnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE)
( C2 ]* b k3 I# Q( b
' {+ F \/ L! w
# 计算在2000天的模型预测概率,这一步不要也行,看你使用哪些参数
: |: \' y5 G/ T7 J7 |/ J/ }. A
p.std <- get.risk.coxph(mstd, t0=2000)
5 N; V6 L O6 a: k8 j) `5 y0 M
p.new <- get.risk.coxph(mnew, t0=2000)
# x5 ?8 _5 u1 A% [9 Q9 |$ k
1
+ m8 U- _- q ?7 y! ^" l6 c- u3 `0 F
计算NRI:
! z8 p- E, T) v+ t( E- u
3 D) C: X$ X7 i% {
nricens(mdl.std= mstd, mdl.new = mnew,
: I4 n$ M! v. Z
t0 = 2000,
9 Y8 Y! u+ a# v. P) K
cut = c(0.3, 0.7),
2 D( g. C. q1 p. ]; |* y
niter = 1000,
, i+ f1 x5 W" R8 ~+ Y" r
updown = 'category')
0 d* H0 Z. H0 E! v8 J8 h% f: n
0 S+ l! Z5 T9 ?8 ]/ P# H) o+ l7 ~
UP and DOWN calculation:
8 u' ]' J: Y, U7 }
#of total, case, and control subjects at t0: 312 88 144
) ~; _4 A" ] n* w0 v2 V6 Q+ y
& I9 W G% Z0 `2 G
Reclassification Table for all subjects:
) t# d" H! U; ]( U, @$ C5 Y' t$ m
New
" Y' f; a9 R& E7 O3 f8 X, |
Standard < 0.3 < 0.7 >= 0.7
; D: j2 `5 } s$ j k
< 0.3 202 7 0
& m+ Y& v r, Y: v8 i/ a
< 0.7 13 53 6
# i, x8 j3 m+ ?$ l: {
>= 0.7 0 0 31
+ z3 N Y7 C) P! ?
$ S N. U0 x: i7 W/ F3 y6 u$ p
Reclassification Table for case:
- R8 O3 v% o# F1 ], _# P
New
s9 c) v5 V. S7 [
Standard < 0.3 < 0.7 >= 0.7
% g5 _! h4 B* E! Z6 N% x( K
< 0.3 19 3 0
4 g/ o5 v/ u: T+ N; w2 V# ?
< 0.7 3 32 4
]3 T R1 q7 g* ?/ g
>= 0.7 0 0 27
3 N8 i' m0 r6 O1 V0 |( U
7 O9 V+ _0 j$ S' y/ Q
Reclassification Table for control:
) ^# O/ E& k) ~
New
- t6 _5 F( ], s. o7 R; A3 |6 q
Standard < 0.3 < 0.7 >= 0.7
6 b: r. z5 M$ e' g5 z# C( U4 F" d+ \
< 0.3 126 3 0
4 c8 ~0 t/ A D
< 0.7 5 7 2
* ~6 x9 @$ _9 _; K# ^
>= 0.7 0 0 1
' ?: B9 U$ {7 `+ ]8 w# R% a. J. ]
. n7 ^' o: J) p$ o7 q
NRI estimation by KM estimator:
4 Y, m2 c2 `. f8 H
: j* M( ~* a8 W, a, @3 W
Point estimates:
$ Y- M5 v7 U7 u& ^* f4 Q( N( D* I
Estimate
2 z+ O) I) L* q. d$ o1 y, y
NRI 0.05377635
: U9 H' l8 Z, {2 J
NRI+ 0.03748660
) e6 f$ S4 B$ Y( }5 B
NRI- 0.01628974
4 G8 N9 W0 q$ H: {% {- ]2 Y0 m K
Pr(Up|Case) 0.07708938
* S0 |8 _$ k0 e) e& K* _" w+ A
Pr(Down|Case) 0.03960278
) c+ z% z0 E* }* z5 ~# x
Pr(Down|Ctrl) 0.04256352
' e$ z: Z0 J2 O; A
Pr(Up|Ctrl) 0.02627378
: m4 }4 G7 P7 o; E' ~2 E% X: H
2 E* Z( c" @4 o0 I( P( s5 t
Now in bootstrap..
+ y5 T, e% b6 a, H( W: A3 Q- v9 Q
& j& Y( B: `' K- u- e% N, r
Point & Interval estimates:
0 Y; a& B% Z i/ E% D1 \
Estimate Lower Upper
8 g! g7 i0 H5 w/ `3 b
NRI 0.05377635 -0.082230381 0.16058172
9 {, \" d) O f0 [" X( S
NRI+ 0.03748660 -0.084245197 0.13231776
- B* K6 ^, V r/ a( e
NRI- 0.01628974 -0.030861213 0.06753616
5 |4 F/ d9 t) L: P6 i/ t6 P
Pr(Up|Case) 0.07708938 0.000000000 0.19102291
& T* t* ]' x3 l+ x
Pr(Down|Case) 0.03960278 0.000000000 0.15236016
4 r& B! x" B+ F) C$ L H9 D
Pr(Down|Ctrl) 0.04256352 0.004671535 0.09863170
4 N. A U0 P/ Z, V* V' N, J
Pr(Up|Ctrl) 0.02627378 0.006400463 0.05998424
, F% ~' Z3 q, |# {- o" B: V) m& }, `
* r5 c! L E2 ^& i
1
5 E4 _( Q6 S9 P
/ u7 e5 \3 a ?( z L7 l" b, h q
Snipaste_2022-05-20_21-49-38
, W) V. F0 G5 {) n
结果的解读和logistic的一模一样。
8 W) _ d; y) f* C( R' F
' _7 ?: H, d: Y
survNRI包
" A# s2 g7 p8 c
# 安装R包
9 _5 q5 ^3 I( N: b4 r' A
devtools::install_github("mdbrown/survNRI")
- n$ N2 F$ z7 M0 [2 s6 ?. p
1
3 O6 i+ n0 E% U. n D
加载R包并使用,还是用上面的pbc数据集。
7 [# p2 S8 a8 t+ b- @7 ^# I c
& U$ v7 ?* L' g5 Z7 G9 v8 `% C
library(survNRI)
6 k4 R5 \* T6 G) V" j% L& g
1
X: t W; X2 t& C" J q3 p/ N6 R
## Loading required package: MASS
9 d- H0 P4 ^8 _/ M/ m
1
7 _! h. ^) s2 s( t+ `
library(survival)
0 j$ M- y9 }9 P6 ? H0 R
% Q1 n4 V0 j/ _& s. {) \
# 使用部分数据
) V' N/ {0 X; s: Y' P% `
dat <- pbc[1:312,]
" O& I5 [( Q: U; m, H3 z/ d- b1 f
dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡
: K2 j, H' Z! I: G0 V8 `
0 o3 j+ W# }$ `
res <- survNRI(time = "time", event = "status",
0 q1 o% I0 _, y9 p& D
model1 = c("age", "bili", "albumin"), # 模型1的自变量
% a3 T9 [0 M' e6 z7 |- G
model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量
, g+ s! K3 \5 e( p) w
data = dat,
0 o- U% g7 o" e! B7 I# I. M
predict.time = 2000, # 预测的时间点
! D5 x2 B- U3 W1 e; \2 b) d
method = "all",
- ]2 k U( R ?3 z! a
bootMethod = "normal",
. j. f) W- P+ ~' C9 a3 t
bootstraps = 500,
# G+ f% s& z) r0 ^( j) p
alpha = .05)
4 J9 a4 G, \* R9 T2 i6 ]
' A6 m! x$ o0 p6 Y Y
1
9 }/ f) O6 N# K6 z
查看结果,$estimates给出了不同组的NRI以及总的NRI,包括了使用不同方法(KM/IPW/SmoothIPW/SEM/Combined)得到的结果;$CI给出了可信区间。
8 T6 a6 C$ B7 Y4 o6 m6 o! q( |0 X5 E! Z
9 s! g0 g# A3 H- o+ }% o6 U2 u% o8 q
res
3 T ?# U) t. ]7 I2 Q( ?. W
1
) H5 t0 e) [# H* F: p: Z/ D
## $estimates
0 R, L; ]% s V. [
## NRI.event NRI.nonevent NRI
5 |( Q, k" h4 D: V$ `' v
## KM 0.20445422 0.3187408 0.5231951
5 z/ r) A( d, V2 J
## IPW 0.22424434 0.3273544 0.5515987
. v. z6 L" `' V4 V0 Z) \' B
## SmoothIPW 0.19645006 0.3144263 0.5108763
/ c0 `4 h7 A4 k8 j7 }
## SEM 0.07478611 0.2632127 0.3379988
" N2 H6 G: i7 z' b! e
## Combined 0.19633867 0.3143794 0.5107181
# H- o) ^; V; _( a9 j& Z2 W. t
##
& V) E5 E/ B4 g
## $CI
* Z' U5 h7 V- x
## $CI$NRI.event
$ T, P( W0 w- `1 E9 T
## KM IPW SmoothIPW SEM Combined
1 b% e1 { F' k; \; h
## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723
0 y* f8 m; ?( h& e. t5 s
## upperbound 0.44806768 0.47033936 0.44014214 0.2658309 0.4400496
6 F0 N# h5 ^4 b) g8 h6 h: W$ _2 t
##
3 K* o. m# H( c
## $CI$NRI.nonevent
. _5 B3 u f" O. r0 H: A8 o
## KM IPW SmoothIPW SEM Combined
' D( s. ]6 {# o) k, W( s$ H
## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.1286426
3 G) r2 @7 V( N! Z1 N% g
## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.6964549
* [" O @ J2 t) h9 V. p1 J
##
" e- O1 K* [; H
## $CI$NRI
3 M* g0 k6 A' G6 y7 V \
## KM IPW SmoothIPW SEM Combined
' _% H) I& N3 @( X/ w3 f
## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409
% J- T4 e2 _$ U& A; X, @
## upperbound 0.89306122 0.92464359 0.87970125 0.64253510 0.87953153
/ z, }& V( S6 S# r/ d2 E
##
H# i, A3 e/ d: h0 Z
##
1 T/ {- o# j- m8 g
## $bootMethod
( x) }1 X, d% l; x) T
## [1] "normal"
3 r3 V# s( S! r2 ~) e9 b
##
( m0 K) S2 I$ p% @ a
## $predict.time
1 X5 ^1 c) L7 [# T/ I& N/ |
## [1] 2000
+ F: ^- {$ E! W, h9 M4 O. }8 s, O
##
3 u% K: _1 r. q' _) U' o" P
## $alpha
. y- r' Y M4 V
## [1] 0.05
% ?" ?1 K+ a5 v* a; K
##
5 j7 c. O5 ^ |' F
## attr(,"class")
" B: H5 s8 K5 _
## [1] "survNRI"
* Z% }9 `, m6 m1 P" }, U
, H& u# J# B d( B" A! [8 {
1
1 ~- D4 |; W8 X$ @$ n$ n2 d0 \1 R8 H
OK,这就是NRI的计算,除此之外,随机森林、决策树、lasso回归、SVM等,这些模型,都是可以计算的NRI的,后面会继续介绍。大家如果有问题欢迎在评论区留言。
. ^3 T. k. w( y# ?% S- b, F
# D9 ~, X; e, [/ V4 p
本文首发于公众号:医学和生信笔记
7 ~! O2 f) }+ m# C
J; j3 W! R, O5 z0 K0 r7 b
“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
' O% g6 _0 B3 _# e
本文由 mdnice 多平台发布
" B7 O* w7 i. i* z) \! Z
————————————————
8 d% {( k& e7 n; {/ N
版权声明:本文为CSDN博主「阿越就是我」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
! p' R* c, o+ {+ y3 k
原文链接:https://blog.csdn.net/Ayue0616/article/details/126768006
/ ]$ }8 ?9 z2 O
- D' t# q8 E: z/ ` l7 h
A0 k5 _! d( M( R/ m5 B& z
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5