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