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