数学建模社区-数学中国

标题: 净重新分类指数NRI的计算 [打印本页]
作者: 杨利霞    时间: 2022-9-12 18:43
标题: 净重新分类指数NRI的计算
$ Z8 i  S% a3 F+ k
净重新分类指数NRI的计算
/ T) L: g  N- d+ i“ 医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。! p0 R+ a4 l3 h% o& M
NRI，net reclassification index，净重新分类指数，是用来比较模型准确度的，这个概念有点难理解，但是非常重要，在临床研究中非常常见，是评价模型的一大利器！& o; j% r' D0 {; O8 q, S

, ?+ \# x4 ~2 n; K$ _$ N在R语言中有很多包可以计算NRI，但是能同时计算logistic回归和cox回归的只有nricens包，PredictABEL可以计算logistic模型的净重分类指数，survNRI可以计算cox模型的净重分类指数。
9 ^9 Z9 N9 F8 w0 k) q) V
# c% R0 O* _0 e: P2 ilogistic的NRI
" g6 w$ L; {4 G3 qnricens包
# J  C) l) C" i4 G( g, X( m. UPredictABEL包' L' i+ c0 \$ [6 y% h
生存分析的NRI' ?3 j) Y$ C2 w5 ^9 C8 [- V; i) l
nricens包9 Y, j* j5 |6 U! M( y6 g
survNRI包
8 q- [$ {5 G" ]) y1 z# N- ]logistic的NRI
  }5 C4 t' b3 V5 a4 Snricens包
$ ]' _9 q% H5 l  R- y! @' B* k#install.packages("nricens") # 安装R包* p6 h+ s# e2 M' U1 h8 x. Q, C
library(nricens)
  h8 Q! a) B% t6 j: c1
5 F6 M$ M4 Y& B& g## Loading required package: survival( U! K5 e/ g# v
15 v3 P9 I. f* L, f% a3 ]7 p
使用survival包中的pbc数据集用于演示，这是一份关于原发性硬化性胆管炎的数据，其实是一份用于生存分析的数据，是有时间变量的，但是这里我们用于演示logistic回归，只要不使用time这一列就可以了。& i8 {0 g1 l/ P$ F& S3 x- D

3 C- q2 X- k2 Alibrary(survival)
4 c$ u3 ^: |( h" V" C- C) G/ g$ T! A  ~' Y
# 只使用部分数据6 N$ _6 ?. n7 W
dat = pbc[1:312,]
7 o% K. T9 O' ~/ u/ V! Q/ P' Zdat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]
$ K3 X' N) E& V) J( t" U7 L6 Z1 j
# Z$ }4 Y$ I, D. \str(dat) # 数据长这样1 Z$ q6 S$ u) V) w2 S
1' j' ?- o: E! Z$ [9 p- R
## 'data.frame': 232 obs. of  20 variables:
2 j" Z) F# S; `7 g3 a. z, o" u##  $ id      : int  1 2 3 4 6 8 9 10 11 12 ...
: c7 B; r6 [: r##  $ time    : int  400 4500 1012 1925 2503 2466 2400 51 3762 304 ..., |; P8 v4 y6 @
##  $ status  : int  2 0 2 2 2 2 2 2 2 2 ..., M" k. k  M* B6 w5 f8 z2 u. k3 m
##  $ trt     : int  1 1 1 1 2 2 1 2 2 2 ...
5 `! D7 u. T% Y! s7 J6 p##  $ age     : num  58.8 56.4 70.1 54.7 66.3 ...
( `; N6 r  B6 ~$ S& D9 P+ J( G##  $ sex     : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...! U# l  t# _8 T) }7 A
##  $ ascites : int  1 0 0 0 0 0 0 1 0 0 ...
1 n: p7 J1 b6 m9 r! `, k3 L2 s# u##  $ hepato  : int  1 1 0 1 1 0 0 0 1 0 ...  d7 ^6 t4 z* \1 f: A
##  $ spiders : int  1 1 0 1 0 0 1 1 1 1 ..., g& Y3 O) l$ Q) v; L
##  $ edema   : num  1 0 0.5 0.5 0 0 0 1 0 0 ...) {2 m8 C! j5 E% L" E' A# N& d! K
##  $ bili    : num  14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...; ~0 u8 ]+ u7 @0 r
##  $ chol    : int  261 302 176 244 248 280 562 200 259 236 ...
: Q% i# @3 D* d# v$ y5 m6 Y##  $ albumin : num  2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...
' \7 }9 p  D/ Z##  $ copper  : int  156 54 210 64 50 52 79 140 46 94 ...
5 N7 ~7 H5 w$ \3 d' h##  $ alk.phos: num  1718 7395 516 6122 944 ...
0 \& P' _' Q2 H: {# G+ R2 e0 c##  $ ast     : num  137.9 113.5 96.1 60.6 93 ...8 t, ]+ M/ J( y7 h4 M6 ^' I
##  $ trig    : int  172 88 55 92 63 189 88 143 79 95 ...4 s. H7 e3 H7 x- e
##  $ platelet: int  190 221 151 183 NA 373 251 302 258 71 ...
% k5 l1 o" e, K' u- P9 A, {0 ~##  $ protime : num  12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...2 R1 }: [: ~$ X1 x( a( R$ C, O6 i
##  $ stage   : int  4 3 4 4 3 3 2 4 4 4 ...& B$ a) {8 v  Z+ `; H

+ ^& w0 m) ~  T. S- u1 N1
6 J: Z1 q6 w  I3 `3 T( m9 U% qdim(dat) # 232 20
/ k# m9 G0 N* @6 u1% K0 ~6 }2 V( w( x
## [1] 232  20& q0 C, Z) ~& a/ x5 p( k
1
3 H+ Y- v  Z7 i! e* Y9 _* [然后就是准备计算NRI所需要的各个参数。
: R% _+ ~" I, k1 I6 F! K; G) G6 B4 T' W) l% @/ p
# 定义结局事件，0是存活，1是死亡
: V' c+ z! }6 a1 @' mevent = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)
, {* }' k3 A3 Z/ J5 @7 M0 P* @4 \6 w: _+ H
# 两个只由预测变量组成的矩阵3 I/ ~  t5 S- r! n- e, u8 O4 C
z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))5 Y- i) `, {; l% B
z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime))). s/ w% E3 y' c4 X- {

: t, V% q# S* x# 建立2个模型
8 h' }  |5 l. Q0 D0 }mstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)7 n3 ?: C  |. r
mnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)
8 g- R% c7 j: j7 W& L, C) R# T8 g) y
# 取出模型预测概率& ^1 I/ C, D+ b' {6 F* n2 W
p.std = mstd$fitted.values6 p( ?( O/ s' r7 ?! b
p.new = mnew$fitted.values& P3 g: ]. r4 s" p+ |5 [( I
+ ?' X, d' r' R: K$ `, W
11 h- S8 N4 l0 ]4 Y! `
然后就是计算NRI，对于二分类变量，使用nribin()函数，这个函数提供了3种参数使用组合，任选一种都可以计算出来（结果一样），以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。
* U: F4 K4 j1 N, w+ Q
! `# s$ h3 w1 Z7 v) l  X# 这3种方法算出来都是一样的结果
/ E% K3 ^% l2 K0 @) a" v; _3 z( o3 @4 }$ W
# 两个模型
1 ~/ @9 v0 K7 ~. nnribin(mdl.std = mstd, mdl.new = mnew, ! f: P: {! U0 G; {7 E: K+ ]
       cut = c(0.3,0.7), - m! z4 q+ o/ u- y
       niter = 500,   \+ f2 B: C) ]7 O/ k
       updown = 'category')
8 p3 J% ]$ x4 g2 b. i
( H5 a0 O- Y' `( N& V' E# 结果变量 + 两个只有预测变量的矩阵
7 X& p  }) j' D4 y$ @nribin(event = event, z.std = z.std, z.new = z.new, % t& p, h$ {1 C6 [) o3 ~* X
       cut = c(0.3,0.7),
8 l) R! t% R$ c  f       niter = 500, 9 g" y8 p6 P# _0 [: o' o
       updown = 'category')
/ [$ S% M7 B! w: Z, c3 N; A6 p8 U" {* f& X/ J$ j4 m  e) f
## 结果变量 + 两个模型得到的预测概率
' M# M) Q$ B: u! j# wnribin(event = event, p.std = p.std, p.new = p.new, & s* e! e3 R. m( \8 {8 R% {
       cut = c(0.3,0.7),
% b* f3 }7 J' z" K2 z1 N       niter = 500,
8 B3 ]6 O) {# a       updown = 'category')
7 n4 ^  y% v+ ^
" i0 S: o, ]! P0 k  C- I1
8 Z% X0 B, W" o* f其中，cut是判断风险高低的阈值，我们使用了0.3,0.7，代表0-0.3是低风险，0.3-0.7是中风险，0.7-1是高风险，这个阈值是自己设置的，大家根据经验或者文献设置即可。3 o2 |7 o$ u- x- O' n

, |1 Z9 L$ ^5 q2 P' w' Sniter是使用bootstrap法进行重抽样的次数，默认是1000，大家可以自己设置。; b% ~3 T' b% f% n' s1 S
  u! p' {% X; U+ a$ V" o3 H
updown参数，当设置为category时，表示低、中、高风险这种方式；当设置为diff时，此时cut的取值只能设置1个，比如设置0.2，即表示当新模型预测的风险和旧模型相差20%时，认为是重新分类。5 ~: B$ c/ Q" o0 {! `0 C+ }& P

4 R5 V6 G+ z; d上面的代码运行后结果是这样的：
* g; J8 h3 ]% r4 R' Z8 W% f
# B+ }* ~- Y, z1 Z' k. ^UP and DOWN calculation:5 x/ ?( `5 v8 Z6 p
  #of total, case, and control subjects at t0:  232 88 144/ ~, V' a9 M8 H5 o

6 d# N* {" m2 H/ A( G+ V  Reclassification Table for all subjects:$ ~0 r. O4 c6 l9 y' }! V
        New6 o, @: e1 a; x$ g5 F; {: `0 y: d4 C0 G
Standard < 0.3 < 0.7 >= 0.7
7 Q! z9 k: G2 d: _  < 0.3    135     4      0+ O2 v7 A" f# z  y# C1 l
  < 0.7      1    31      4
& K' K! I( B! L. P' K6 Z  >= 0.7     0     2     55
$ T4 D% G( E" _" d( v
  C, P  B* D* o) `  Reclassification Table for case:! n! b) h( X* f: n1 Q" o* Y- o
        New' {% R3 w4 `8 F1 i  z: x8 n* c( j
Standard < 0.3 < 0.7 >= 0.76 _) g0 F& T) |! B$ R/ V
  < 0.3     14     0      0; Y. }* R- F8 ?" {( }: o1 F7 I
  < 0.7      0    18      3
. D; d; P) e( H  >= 0.7     0     1     52
) {* C5 ~: J! p" Z( P9 H8 f/ t) T7 U) t8 \4 k
  Reclassification Table for control:9 J) Y* S, A5 j! L* j  x! V
        New
& E) t3 D3 f& a) e: D1 d- z+ o/ SStandard < 0.3 < 0.7 >= 0.7
( }, E* w1 i6 H* p$ D0 h  < 0.3    121     4      0
$ b( p" r  o2 b+ a: \/ K) F7 n; R  < 0.7      1    13      1
. Z2 L  d+ T' ?1 G; X1 ]* L  >= 0.7     0     1      3
3 z. ^' v0 n4 R0 o+ o! k& }0 _2 p  S. R: l* g" E& ?$ k& _
NRI estimation:1 N- D* P, B; P) C' G
Point estimates:
2 [- q# u* u7 a1 ~5 n# z) B; w                  Estimate
8 e& h: d0 T- r% nNRI            0.0018939391 V. I+ R$ L# C) b- {$ x! k: K6 }
NRI+           0.022727273' n# {: C2 |- S2 D1 S& U
NRI-          -0.020833333; \& c8 T5 R' h( x4 |) B# u7 P
Pr(Up|Case)    0.0340909099 ?$ F& N5 x8 X& I  `4 x, m
Pr(Down|Case)  0.011363636
) ^0 @3 x( c! B  F/ R! SPr(Down|Ctrl)  0.013888889# U: n8 [& ^7 ^( @0 B% @6 A& i* [' w
Pr(Up|Ctrl)    0.0347222226 U3 U1 D& c" U7 [
6 o2 o+ H' x- {- Z
Now in bootstrap.., b; T: b, a8 ^& d

5 q+ [0 Q/ _, N6 I" IPoint & Interval estimates:. H4 y2 {, m6 F- M
                  Estimate   Std.Error        Lower       Upper8 P5 F2 \8 x/ B  H: k2 @. ^
NRI            0.001893939 0.027816095 -0.053995513 0.055354449
, B2 f1 d+ {: S9 M) G& v$ s( ZNRI+           0.022727273 0.021564394 -0.019801980 0.065789474
% i$ t% w0 O8 H$ I- xNRI-          -0.020833333 0.017312438 -0.058823529 0.007518797
, ~7 c2 D( A" KPr(Up|Case)    0.034090909 0.019007629  0.000000000 0.072164948: \* D) H0 I$ w4 z$ }3 w9 p
Pr(Down|Case)  0.011363636 0.010924271  0.000000000 0.039603960
6 q) |( W, V$ r& h2 a) _" q) xPr(Down|Ctrl)  0.013888889 0.009334685  0.000000000 0.035211268
+ e- H9 C8 D( v* @, ^" S& gPr(Up|Ctrl)    0.034722222 0.014716046  0.006993007 0.066176471( H) Y1 B7 ]) S) A3 T0 J! ^1 y

" ^1 l- B! I9 Y* T" M7 I1
: p: O2 J; c8 g0 r+ c$ {( }& d首先是3个混淆矩阵，第一个是全体的，第2个是case（结局为1）组的，第3个是control（结局为2）组的，有了这3个矩阵，我们可以自己计算净重分类指数。
; I) L' k4 Q7 C% z1 y2 t, `$ m8 e
看case组：
% B: G4 f! g! C" I  U
# ^" W+ w: G' b5 v净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273
/ t+ M  f- j3 J# r3 q! }) D* s3 @) w2 Z6 L( l
再看control组：5 w( l4 e& b; O5 a% N; f% z
$ Q0 n% A0 U0 _
净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333" b5 A; O2 |6 u4 e
# l7 [- _0 J) ^* `$ |) `
相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.000315657, c4 G2 d4 V. m% {2 \% n9 \$ p

; I$ Z7 n# K1 `  ^  m/ C再往下是不做bootstrap时得到的估计值，其中NRI就是绝对净重分类指数，NRI+是case组的净重分类指数，NRI-是control组的净重分类指数（和我们计算的一样哦），最后是做了500次bootstrap后得到的估计值，并且有标准误和可信区间。0 e, J  ]' x" Y" k6 }
7 X& y' E, R1 p: u1 ~9 _0 c
最后还会得到一张图：
6 x" D- D% Q2 a9 L* w( e1 ^6 r0 x4 Q+ P6 w' E: x) Y4 d
这张图中的虚线对应的坐标，就是我们在cut中设置的阈值，这张图对应的是上面结果中的第一个混淆矩阵，反应的是总体的情况，case是结果为1的组，也就是发生结局的组，control是结果为0的组，也就是未发生结局的组。
+ ?1 v- ]$ \# [; S
$ G8 ]- H+ q3 m9 T) JP值没有直接给出，但是可以自己计算。0 T$ T4 ~- H5 `$ Z% N& S# [

( p% l& `; }. \' m1 J6 T# }# 计算P值) [7 t% q3 ~- R6 b0 N
z <- abs(0.001893939/0.027816095)1 a3 J; |( N) h: k. P
p <- (1 - pnorm(z))*25 A! x& M* l$ ~  l$ [& ?* X8 X$ X# j! D
p( s2 c! q" w+ M  K/ d  O; e' R# T% e
1
( k7 T' Q0 h  |## [1] 0.94571571 e- [4 F  K" @3 b. Z4 ~
1
) A2 T; E9 k5 l. [% rPredictABEL包  U7 @& X& a$ {3 H3 D9 m4 X
#install.packages("PredictABEL") #安装R包
$ L9 L" V+ u, [library(PredictABEL)  9 {# z1 Q# m0 X' U, `! z& o

7 Q; \9 g5 J* p* f4 _5 [# 取出模型预测概率，这个包只能用预测概率计算$ c( S4 w9 T3 H2 A2 K0 X+ |- z
p.std = mstd$fitted.values
' {4 E6 K: }6 X9 L0 x! @" a% c! Rp.new = mnew$fitted.values 4 ~8 P. \. ~3 @* L5 _1 [4 O
1- l" D' s( C+ k
然后就是计算NRI：
* j5 w* _1 J4 W* x9 K+ u( h# |
dat$event <- event
7 a5 z  C8 |' r1 P# \3 L6 l6 K7 C! d2 F3 e7 u1 a- D* n2 Y
reclassification(data = dat,- H/ k% j: `3 o. f5 Y# C
                 cOutcome = 21, # 结果变量在哪一列9 D' ]- M$ D" P6 n! e
                 predrisk1 = p.std,
4 }: Z+ l) ~" S6 ]                 predrisk2 = p.new,
. o9 ], \" T  y8 `$ ]! N                 cutoff = c(0,0.3,0.7,1)
' S6 R1 r$ O3 n, Q                 )/ h! k6 C6 k$ E% e; q+ ]
1# R& |9 o1 _  @5 E8 o; B
##  _________________________________________* p. ]: e) F- I5 X
##  ) V9 N) l( z' D, S
##      Reclassification table    & w. z( L  W! ]
##  _________________________________________# s5 f0 k8 E% z3 z  G7 E
## 9 r$ v( U6 e! O4 T
##  Outcome: absent 7 `* M6 M+ {0 x% Z" t1 r: y
##   , U3 g# n7 G7 b6 r
##              Updated Model: b! Q& }6 e3 |5 V/ \3 Z6 n6 Z
## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
3 J) ^' e" `1 A4 f( P##     [0,0.3)       121         4       0               3, |) a: D' k4 _) O. M7 Y( ?, J6 Z
##     [0.3,0.7)       1        13       1              13
% f' \9 n* E( Q4 O" [  n##     [0.7,1]         0         1       3              25; e$ }4 C6 [5 P: N' E9 |
## / ]# n& v; _- s2 K& I
##  
% ]& O# ~, r  z1 M! Q' m' `##  Outcome: present
& a" w7 ^) X+ H7 ], x##   , H( n4 ^  H$ n$ G- l
##              Updated Model
! \7 Q' m  o" X0 o5 n* c& g## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
0 n( f* ]% Z. R4 b##     [0,0.3)        14         0       0               0! p$ i" J# E5 ^: S" f
##     [0.3,0.7)       0        18       3              14
( K- {1 U6 w) {2 k8 I. {* u##     [0.7,1]         0         1      52               2# y! [3 H5 y  _$ S8 ]. a5 o. a" u7 o
## ' _& Y1 ~' d5 n- v- ~: G' r+ o
##  
9 }4 z. B3 V+ m3 J% }: p! d! L5 \##  Combined Data 4 t, Z2 ]/ O( U7 @: |
##   4 U1 K' t- e# O
##              Updated Model3 b7 r* V' A! g- c; l& @* m7 h' Y
## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified- l. w) X. ]( v! b7 @" H& {
##     [0,0.3)       135         4       0               3
3 Q  v. \* S' s4 u/ g##     [0.3,0.7)       1        31       4              140 m  ?, g/ l- ?5 S% \  O8 n
##     [0.7,1]         0         2      55               4
. l4 G+ {  v% X4 Y2 f, e##  _________________________________________$ Q. R+ z3 J; o/ a5 U
## + n) B% \  ]# c8 ^3 ]0 g
##  NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806
- @+ w; C  P: G4 v7 I) o5 Z##  NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048 + M  {  `4 j' H4 h$ d9 n* @9 W  Z
##  IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.283964 @& T- f1 v0 _* W/ ^2 P7 k
/ J" d  J- ]6 N/ M/ u; {# `
1$ O! V3 ?# z5 J  }. |4 w
结果得到的是相加净重分类指数，还给出了IDI和P值。两个包算是各有优劣吧，大家可以自由选择。( g. W& k" w' c; T# R

7 M( Z/ n+ d& R+ M  b& ~9 r生存分析的NRI- D6 E; `9 ^9 s" N  D
还是使用survival包中的pbc数据集用于演示，这次要构建cox回归模型，因此我们要使用time这一列了。
( @* F( b9 d. E' l
/ n! d- \0 n. o7 |: Pnricens包- m8 Y4 ~$ Y  i# q- a. }' l
library(nricens)0 L3 O9 x+ B2 |& q% K
library(survival)
; X6 }/ D' ^: l: c0 G( F* `7 o9 s5 J3 H4 _% r0 {# r
dat <- pbc[1:312,]
" {7 p  V7 X9 z3 H- Ldat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着，1表示死亡6 r8 e7 k* z- o( i- b
1( X9 L* }3 d, S6 |
然后准备所需参数：
$ H$ {6 x' `4 f* k! b5 v
- ~2 w* @, w& T/ a% i+ L# 两个只由预测变量组成的矩阵0 A) c# ^5 x. f
z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))3 D5 G6 k' q: N6 t) t
z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))6 b1 k9 [7 z4 E

( `3 m7 x5 _5 ?" c4 t/ X# 建立2个cox模型& l$ d  T4 f: Q9 ~) i
mstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE)
; A0 F' y. a7 n" I# k5 C$ cmnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE)
. W  m, e) f+ j: z3 ]9 N- X1 i" z/ N
# 计算在2000天的模型预测概率，这一步不要也行，看你使用哪些参数
0 P# Y' [6 x7 ~+ Q8 Cp.std <- get.risk.coxph(mstd, t0=2000)6 G0 X/ N/ [6 r% l0 o' }
p.new <- get.risk.coxph(mnew, t0=2000)
! T' o% k8 E9 t, [1% R1 O0 }: k4 s" W
计算NRI：
: G# w! x6 b6 [3 |, B
' f& T6 d: c' I9 Pnricens(mdl.std= mstd, mdl.new = mnew, 7 l% p( R5 t3 d  W; p0 ?
        t0 = 2000, + i: I* \; l" O- G8 D# H
        cut = c(0.3, 0.7),$ C3 M* @) i/ u/ {
        niter = 1000,
4 G  k. t6 S$ S" E; B  x  x        updown = 'category')
" ?, X, P: R6 {- v% O, t9 P1 K0 w( q
+ p( _* y3 e" ]UP and DOWN calculation:* P: v- J3 ?; J
  #of total, case, and control subjects at t0:  312 88 144
# c( @9 C9 L; A
+ v: Z& C+ z! a; x2 f' s8 q  Reclassification Table for all subjects:
# @+ z6 Z$ E# E& n        New( S7 {8 X' J+ q3 U+ N: A" O$ u
Standard < 0.3 < 0.7 >= 0.70 d  j1 ^% U0 A9 f
  < 0.3    202     7      0. q4 e+ a6 u" |# V1 @& A
  < 0.7     13    53      6
) b7 }1 l( a* x! y6 T+ {  >= 0.7     0     0     31
! g4 H7 Y" v& \7 S# T% W: u; Q$ b, O5 t
  Reclassification Table for case:
. O% y) }( `& W% G        New( J, m6 y' E, q* d% \) b& D+ ?
Standard < 0.3 < 0.7 >= 0.74 z( o/ p" A2 Y+ U: b
  < 0.3     19     3      0
% g) M& W" C# @) b  T5 B+ }  < 0.7      3    32      4
2 F9 T& t& `$ H4 g, R2 Y) v' Q* g  >= 0.7     0     0     27
+ {0 r' s5 X% a2 i- n0 K9 A) o$ Z# }4 c3 D: e; _
  Reclassification Table for control:, G* c+ c6 s, ?  U. ?% W9 ?
        New; C$ o1 n7 x( P) y0 T9 l
Standard < 0.3 < 0.7 >= 0.7
( L5 `% j8 Z0 W  L  < 0.3    126     3      0; p/ X+ D9 B7 v  L: C4 ^! A
  < 0.7      5     7      2
2 d) S. v3 y0 Z# ]0 ]- e2 @  >= 0.7     0     0      1' s5 H2 r+ g5 X) {) E& B3 `0 q2 T2 @* m

, T) K0 g+ \3 w9 u) J3 b, PNRI estimation by KM estimator:; h! k: v# p, E- z& _5 a
9 a" h; h" m) x; r: |# V* P
Point estimates:
. n& H' g9 Z: l* b                Estimate
7 B( t9 b" x' ]7 n* NNRI           0.05377635
0 N$ w$ ?7 \7 y7 p1 W1 i# n6 RNRI+          0.03748660# t" o# M" c4 q* S8 {8 M
NRI-          0.01628974
! r" J- U7 x% i. g/ }- h9 EPr(Up|Case)   0.07708938& l+ r8 q6 A. [) v# `2 n1 w
Pr(Down|Case) 0.03960278
/ J! L5 ~$ b$ @/ zPr(Down|Ctrl) 0.04256352
" I- |/ I/ ]  u1 A4 w/ a7 a. sPr(Up|Ctrl)   0.02627378
- x  f8 J! r- f
, B3 H2 z3 \# l( r- A& y+ ~Now in bootstrap..
6 u" x) m/ j8 W" v& ^2 D& S' d& W# j# Q7 Y+ r
Point & Interval estimates:0 a# s5 d. ~) t0 @$ z& [9 `3 x* n6 ^
                Estimate        Lower      Upper+ _( A. H: r) J4 V1 h6 U
NRI           0.05377635 -0.082230381 0.16058172+ t* f9 O' I9 @) Q* o- x
NRI+          0.03748660 -0.084245197 0.13231776' h* D. w  D, P; s3 P2 X! ~1 c* h
NRI-          0.01628974 -0.030861213 0.06753616
, h9 j8 f1 b0 C8 b: @+ p8 ~$ [; EPr(Up|Case)   0.07708938  0.000000000 0.191022910 E6 Q: q2 D) X+ j5 t8 o! i
Pr(Down|Case) 0.03960278  0.000000000 0.15236016
7 n, D& R' Q+ @+ A2 a- zPr(Down|Ctrl) 0.04256352  0.004671535 0.09863170& T8 N% @  [5 k3 |' i& O
Pr(Up|Ctrl)   0.02627378  0.006400463 0.05998424
8 O& x9 g1 u# G, d0 q& ?. B; i
% g1 ]/ x7 p7 X7 g1
: p! i. p. j6 l2 B/ ^; W
: h% z, V# A! U: ^, C2 uSnipaste_2022-05-20_21-49-387 e! t, a  v  A4 e
结果的解读和logistic的一模一样。6 q5 K8 S" c; L; A! Y" Y5 f' {

$ k: ~; W% Z- S7 nsurvNRI包
- o# G# E8 B3 Q! x2 D# 安装R包
* I9 j7 L% O; Sdevtools::install_github("mdbrown/survNRI")
% @/ z3 v& F% c, R, H: Q15 C' V2 L" s  @- Z" E
加载R包并使用，还是用上面的pbc数据集。" y6 T, R2 W( y5 ~

. z2 q) _/ M$ y1 i: @& Y& Rlibrary(survNRI)) r" o, f1 Z0 b7 h. A' ~
16 Z" l1 l9 C( J: T( K5 C- G
## Loading required package: MASS. T! g! u+ v8 V0 d
1
& \+ q' u: J- _# Jlibrary(survival)6 F% h% z2 C5 K

. T( a+ c' ^  l# 使用部分数据8 T: |9 E$ w/ z  s; q, @6 I
dat <- pbc[1:312,]
. p' I+ H0 Z* ]4 Wdat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着，1表示死亡: A: q3 p% M- ?" p  y* K; D
  ~" P0 o: b; w# h
res <- survNRI(time  = "time", event = "status",
. {: e! e* E& x8 o        model1 = c("age", "bili", "albumin"), # 模型1的自变量
' ~+ ^# |& X' y- m7 L1 O6 Q5 }        model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量
: m& e& J/ A( r3 Y& d        data = dat,
- x; H* F( Y' a( n7 G0 R        predict.time = 2000, # 预测的时间点/ Z3 r# c. Z; g. s0 j
        method = "all",
  k0 x* a8 p* k1 q* y        bootMethod = "normal",  
4 b7 Z9 Q: b8 `        bootstraps = 500,
1 b7 `! X7 r! V. Z/ E" W' M        alpha = .05)0 C0 D, Q% n  u  X
# D! e& y% I2 K* E" ^3 A$ j; v" _
1! j$ B. R/ B- P2 m5 Y
查看结果，$estimates给出了不同组的NRI以及总的NRI，包括了使用不同方法（KM/IPW/SmoothIPW/SEM/Combined）得到的结果；$CI给出了可信区间。/ b0 N$ t( ?9 Z
; N, O4 B: d+ W  ~
res
+ A- \0 f  y) {% N! n16 K. w" r$ N$ m2 y' G
## $estimates
: Y0 H3 f3 z8 U2 `! X, h##            NRI.event NRI.nonevent       NRI% n. t& F4 Y4 n2 O7 R
## KM        0.20445422    0.3187408 0.5231951  l" ^. y9 |( |
## IPW       0.22424434    0.3273544 0.5515987
, P" i" l0 ^& |* [+ w## SmoothIPW 0.19645006    0.3144263 0.5108763
/ z% W2 }/ t' \) N# a& _- \## SEM       0.07478611    0.2632127 0.33799884 t) i, u; G1 W( I2 M  T1 k
## Combined  0.19633867    0.3143794 0.5107181
- V) U: R" u- l1 w, ~  H##
% h& ?2 u  Q- _1 v7 {2 P## $CI
8 W. }& D: ^& l, C## $CI$NRI.event
3 D, u. o% M! Q+ |##                     KM         IPW   SmoothIPW        SEM   Combined6 `) b2 I4 H4 W! w1 C% @$ L* M
## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723# J7 R( M) B" o; N+ [  W
## upperbound  0.44806768  0.47033936  0.44014214  0.2658309  0.4400496) ?& j. S# N$ w$ P
## & P; u4 U+ p/ V9 }2 d% C
## $CI$NRI.nonevent
0 s! |. ~4 L& y- N$ w7 N2 w1 k% a8 X# c# R##                   KM       IPW SmoothIPW        SEM  Combined
+ F( Z1 D8 A2 K; C6 d## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.12864262 |3 r# J8 N( Q: f. Q2 N
## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.6964549
7 w9 l& h2 B4 e; b# |# |# d4 d5 E8 W## 3 B' ^. m3 T' x; U. v9 k7 H: J
## $CI$NRI2 Y& D7 [: N& Z1 H9 b* Z
##                     KM         IPW   SmoothIPW         SEM    Combined1 z( Y6 w$ ~( \& L+ l) \* Y& x
## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409
4 ?0 e7 v8 [7 @$ K) V$ x4 V## upperbound  0.89306122  0.92464359  0.87970125  0.64253510  0.87953153
( ?9 w1 u3 ?/ o: k4 j##
% C) f8 p8 K+ \, K- Z. I1 ~( w## ) V9 x" ?8 H& F& q- ^
## $bootMethod4 z4 N( o7 v& G5 }6 G( c
## [1] "normal"1 I+ D& T0 q1 q
## , R0 H3 `# o/ g
## $predict.time
6 z* Q6 m2 F( o1 a$ _## [1] 2000
+ J: J6 N9 ~0 z2 S##
7 W3 K: {4 q2 j( r& m## $alpha
1 R1 ^9 Z: C& G2 v; ^1 V+ A' O: J## [1] 0.05
& a5 g/ W; q3 [$ _; j## 5 b0 @( P2 o; w, o
## attr(,"class")* _" a" i5 W6 c& b: y- s8 Y1 z
## [1] "survNRI"
% i1 _4 ~, q# p' g. t. R. P0 @* A
9 b/ O' V8 N! I- h1 k4 Y1
: }" K& A( O, L5 m! VOK，这就是NRI的计算，除此之外，随机森林、决策树、lasso回归、SVM等，这些模型，都是可以计算的NRI的，后面会继续介绍。大家如果有问题欢迎在评论区留言。
# m$ ^% c: O! V) e' r8 Y
3 s& A* a: H1 k( L8 Z3 `" y% x本文首发于公众号：医学和生信笔记
7 M4 O+ G) I! F8 e! X
/ v1 R4 s  P. ~0 L9 ~2 ~) N0 o“ 医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
' R$ x- _. U: @8 w本文由 mdnice 多平台发布
# ]% Q; d: m: y, K# m4 e/ F' r————————————————
% [8 p6 |, k5 T/ R3 T0 _版权声明：本文为CSDN博主「阿越就是我」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
' [& e* B' H' r  T( v& Q2 o% L原文链接：https://blog.csdn.net/Ayue0616/article/details/126768006) l+ ~  h5 K+ c" H. S. Q* o

/ T! j' a: `5 k% h: E- L
( Y( F) H+ ^: x! f; ^- P



欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5