QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 3053|回复: 0
打印 上一主题 下一主题

[其他资源] 净重新分类指数NRI的计算

[复制链接]
字体大小: 正常 放大
杨利霞        

5273

主题

82

听众

17万

积分

  • TA的每日心情
    开心
    2021-8-11 17:59
  • 签到天数: 17 天

    [LV.4]偶尔看看III

    网络挑战赛参赛者

    网络挑战赛参赛者

    自我介绍
    本人女,毕业于内蒙古科技大学,担任文职专业,毕业专业英语。

    群组2018美赛大象算法课程

    群组2018美赛护航培训课程

    群组2019年 数学中国站长建

    群组2019年数据分析师课程

    群组2018年大象老师国赛优

    跳转到指定楼层
    1#
    发表于 2022-9-12 18:43 |只看该作者 |倒序浏览
    |招呼Ta 关注Ta

    7 ?+ T1 F3 f% b8 R( B净重新分类指数NRI的计算
    % V* P5 o* ]1 ^9 A1 i, }; E' G“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
    4 ?+ U, H4 L% M9 Q/ n# p" INRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!
    ' Z+ y) ~0 f7 ^
    ; z% ^' c& L7 }4 p# S$ N在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens包,PredictABEL可以计算logistic模型的净重分类指数,survNRI可以计算cox模型的净重分类指数。
    7 h# @3 P& H! S: r1 |# N
    4 W) ~/ ^) Z; ilogistic的NRI
    ) N: s9 G7 q, e; Tnricens包! |" J: U3 a, q* G
    PredictABEL包$ i4 E+ t$ W7 p: G' a/ C2 ~
    生存分析的NRI( Y& R; h3 C8 h+ d
    nricens包
    ( @4 W4 p# R/ w$ C+ G  JsurvNRI包5 _9 O. I* K! I6 a. P" y1 Y, A
    logistic的NRI
    / |8 g2 ]1 N) b( V3 Tnricens包
    9 a  H" k9 k  ^& H3 [0 A1 P#install.packages("nricens") # 安装R包, t) N0 }7 V0 |
    library(nricens)
    / k" P& `8 c1 x. `' I2 r1+ e! D0 E9 H( t
    ## Loading required package: survival5 w) l, b: s7 B. ]: A- v
    1
    / o0 R) _& M0 b0 M2 M使用survival包中的pbc数据集用于演示,这是一份关于原发性硬化性胆管炎的数据,其实是一份用于生存分析的数据,是有时间变量的,但是这里我们用于演示logistic回归,只要不使用time这一列就可以了。; S! M. @2 r' [5 _
    7 q. V7 H- G5 r" o& U: Q' }9 L
    library(survival)1 N# ?- F+ U  \
    # i( u. d! i0 r% O4 t: F  `: O% j
    # 只使用部分数据" @1 h6 Y* V, k6 M2 p0 q; J
    dat = pbc[1:312,]
    4 e. G' d1 _5 g; g+ k! n8 Edat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]1 J3 K4 i1 L; F* P) w

    - E8 V0 s4 j* [* gstr(dat) # 数据长这样1 N6 Y# w8 E' V9 V9 S$ c
    1/ g/ I! x' n  j7 K2 V  C
    ## 'data.frame': 232 obs. of  20 variables:
    ' e+ T; d# N0 `4 U+ @2 W##  $ id      : int  1 2 3 4 6 8 9 10 11 12 ...
    4 F0 ]; {4 j2 n- W* R3 W9 h##  $ time    : int  400 4500 1012 1925 2503 2466 2400 51 3762 304 ...4 _* K+ k* |: Z+ z* K/ j8 Z
    ##  $ status  : int  2 0 2 2 2 2 2 2 2 2 ..., r: R5 X: j' w* z+ h
    ##  $ trt     : int  1 1 1 1 2 2 1 2 2 2 ...
    7 D) ~) R+ u! J5 V2 T##  $ age     : num  58.8 56.4 70.1 54.7 66.3 ...  b+ x3 r. C% w: F# t
    ##  $ sex     : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
    $ c/ e. I$ m3 j! q! W( F##  $ ascites : int  1 0 0 0 0 0 0 1 0 0 ...
    * i3 A# E3 Y+ a. n( }# R3 y, w##  $ hepato  : int  1 1 0 1 1 0 0 0 1 0 ...
    5 [, C0 j$ Q- E# C* R- L: ]* i##  $ spiders : int  1 1 0 1 0 0 1 1 1 1 ..." C0 z  S7 P7 |6 N# h0 w
    ##  $ edema   : num  1 0 0.5 0.5 0 0 0 1 0 0 ...
    ; k5 Q9 B( ^5 Q0 c! O! T9 E9 d) M##  $ bili    : num  14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...: l! P0 y+ N4 t
    ##  $ chol    : int  261 302 176 244 248 280 562 200 259 236 ...! C$ k' O0 X- b2 F2 j- q3 x; t$ h
    ##  $ albumin : num  2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...
    9 v9 M, a4 S' Y; B* v##  $ copper  : int  156 54 210 64 50 52 79 140 46 94 ...+ }9 x+ e9 Y, C4 E
    ##  $ alk.phos: num  1718 7395 516 6122 944 ...
    ) f9 @, q+ M0 P6 c- y. K% I##  $ ast     : num  137.9 113.5 96.1 60.6 93 ...
    ! h! D! \- `! k" T) b6 g+ ?##  $ trig    : int  172 88 55 92 63 189 88 143 79 95 .... a7 L9 l2 u% G0 i6 X
    ##  $ platelet: int  190 221 151 183 NA 373 251 302 258 71 ...
    5 j1 I! w; j6 r' L1 S3 ^##  $ protime : num  12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...; V. F9 s. k* @. l* o2 H, t: T
    ##  $ stage   : int  4 3 4 4 3 3 2 4 4 4 .../ B, M6 H; n6 j1 P) {9 I
    " v+ i, H, X) `9 d: e
    13 i9 J7 I$ M3 F: W- C
    dim(dat) # 232 201 _7 b  c- ^, N
    1
      e. H% [; Y; ~## [1] 232  20
    % T0 e/ y& y: s" [, z: u3 M/ F1" b! V- ^- [( x, V
    然后就是准备计算NRI所需要的各个参数。/ U  E8 W- w+ U! W5 Z- I. {. E

    5 E2 b6 m. G: l+ R0 w  g# 定义结局事件,0是存活,1是死亡
    3 f) R- ~6 S" J0 P5 Eevent = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)6 A# ^4 ~! X4 @. v  w7 s- S

    . n& J( q8 L6 ~' j# }  \$ B# 两个只由预测变量组成的矩阵; e; P0 T1 N1 k0 X+ N6 H0 B
    z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))( O- V9 |9 X0 [
    z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))% U% p  P$ F( C9 ]. G
    7 z. B; W6 n! N$ F2 W* ]
    # 建立2个模型$ v  \- I: m7 A+ B
    mstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)9 j: q9 L' v8 ]! [( \
    mnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)  `2 M* S# ^2 A0 {! g

    $ a$ {# n; K5 n# 取出模型预测概率
    . w; R: c1 H5 ~& s* n( i& {p.std = mstd$fitted.values2 Q$ ~6 p& l. e# r
    p.new = mnew$fitted.values
    ( X( n- l2 a) n- g7 p1 H5 L: r9 H
    1
    ' \2 v) U0 x9 [. f然后就是计算NRI,对于二分类变量,使用nribin()函数,这个函数提供了3种参数使用组合,任选一种都可以计算出来(结果一样),以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。; f  @/ G6 I! ]4 }) R7 t* b7 g% J

    # A, o' s" f1 _7 U6 w7 s: W# 这3种方法算出来都是一样的结果( x5 Y/ E# r7 T! C& U

    7 A* ?9 X' i, ~" l& V9 O# 两个模型0 k- L" s9 w8 ?4 D7 F
    nribin(mdl.std = mstd, mdl.new = mnew,
    1 @( x- U! A6 j  F3 \       cut = c(0.3,0.7),
    ! y9 D' N8 b  s* k8 F       niter = 500,
    0 G, M0 ?* _9 e' P       updown = 'category')' C2 A' I4 ^) V* l% V" `+ Q
    5 V9 `4 M% Z3 I2 p0 q; V
    # 结果变量 + 两个只有预测变量的矩阵
    : `& a1 {/ M! W& w* K$ enribin(event = event, z.std = z.std, z.new = z.new, ) i# d6 |2 j3 T: `4 ]" B7 e( ?* S
           cut = c(0.3,0.7), # M% X2 H9 x+ b& q7 y
           niter = 500, ) ]" g  w! Z6 p( A5 L
           updown = 'category')
    1 w+ j- ]  O& Q9 [; {4 X
    8 U$ `" Q& G- `* A/ J+ a## 结果变量 + 两个模型得到的预测概率
    % @8 ]) ^, F9 e! Q, inribin(event = event, p.std = p.std, p.new = p.new, . I9 f7 v! r6 [3 C
           cut = c(0.3,0.7),
    * M) j# M3 {' Y1 m; t  a% ~       niter = 500,
      f) S( l5 ]8 s( P! D- R       updown = 'category')
    . u& X  J" X, F/ p9 Y. T: ^* V% F
    . S! ]# S/ q- n1
    * K9 X8 J  j9 r其中,cut是判断风险高低的阈值,我们使用了0.3,0.7,代表0-0.3是低风险,0.3-0.7是中风险,0.7-1是高风险,这个阈值是自己设置的,大家根据经验或者文献设置即可。
    " p( L8 K  r" x4 t
    ; b3 R( B( @* e3 J: ^niter是使用bootstrap法进行重抽样的次数,默认是1000,大家可以自己设置。' O' k/ j1 n  r" v. x
    , S5 q" o$ T* M% X/ X2 r) F4 E
    updown参数,当设置为category时,表示低、中、高风险这种方式;当设置为diff时,此时cut的取值只能设置1个,比如设置0.2,即表示当新模型预测的风险和旧模型相差20%时,认为是重新分类。
    & A6 @0 l1 G; G3 z. U$ w, Z" l2 q- R+ V& N0 v
    上面的代码运行后结果是这样的:
    % T# R* B1 \4 s5 ]7 b. ?- ]) Z, F9 O) y& i/ Q# j7 a8 J$ V4 N
    UP and DOWN calculation:
    0 O2 @* @% E  U9 m/ u0 T7 A  #of total, case, and control subjects at t0:  232 88 144
    & i4 W! C% X) b4 E( n
    ; O0 e6 j9 u( p4 k5 q  T8 T7 ~1 _3 u5 I  Reclassification Table for all subjects:
    9 i) \1 m" H' F        New  E. R6 s: H7 V
    Standard < 0.3 < 0.7 >= 0.7* l+ B% \& ~" `% k+ z
      < 0.3    135     4      0: |) a: z! b7 g. t# t+ M
      < 0.7      1    31      49 f  J/ M& N" M9 P; P+ S5 k
      >= 0.7     0     2     55
    / O) W; ^7 Q) q' Y0 o' ~$ h" q" c
    * B/ G* \5 U- e7 ]& o$ i  Reclassification Table for case:# C9 X( C6 Z/ m) w( X  R7 O; S/ |% W
            New3 N$ i$ }8 t0 w. G% ?8 N
    Standard < 0.3 < 0.7 >= 0.7
    1 u6 H# p( W2 q2 o' o$ m/ a; J  < 0.3     14     0      0
    " R% G; s) D* \& w  < 0.7      0    18      3
    ' x# z( R: T, c! }/ a& W  >= 0.7     0     1     52
    2 v3 d$ ~; N7 c' v2 F+ E# t9 w
    2 _/ i1 c5 J- ?' L# f% \  Reclassification Table for control:' p+ J. W, k% U+ y
            New
    0 l- S9 m* S; t" O0 l% v1 n9 E( O% T- pStandard < 0.3 < 0.7 >= 0.7, D/ X1 d$ ^7 ~
      < 0.3    121     4      0" w5 V1 Z/ D7 a. Z. @6 G
      < 0.7      1    13      1
      ]0 q, n% b3 z' W- I2 y0 \- P  >= 0.7     0     1      34 G8 v# g0 H$ G3 \
    * [4 G% M: n) C9 `! k
    NRI estimation:
    % |6 c) g3 R0 z  Q+ APoint estimates:
    & v( F3 h2 N/ R6 q6 e" s- o                  Estimate" }! ^! _3 i' v/ B) g1 _
    NRI            0.001893939
    : [0 F& l5 u# Z* gNRI+           0.022727273' Z6 w5 |/ M6 e  o1 R
    NRI-          -0.020833333
    6 \; e. Y1 u: x1 p$ R& FPr(Up|Case)    0.0340909095 B% G6 s5 z8 @' _" B5 q
    Pr(Down|Case)  0.011363636
    : ^' E: z5 s8 W3 l0 d# _Pr(Down|Ctrl)  0.013888889
    ( K4 s: b3 j6 q* R/ V' jPr(Up|Ctrl)    0.034722222
    $ m5 F' }. W; e: e
    1 |* M+ E( L5 g+ L; Y( B3 CNow in bootstrap..
    ( d/ l! o# d1 O' \+ C/ E
    ; y* D9 l' k- x# N7 \Point & Interval estimates:$ O- ~1 \% n' P6 Y8 K
                      Estimate   Std.Error        Lower       Upper
      K+ k# U' _5 U" A# nNRI            0.001893939 0.027816095 -0.053995513 0.055354449
    7 Q6 ?$ A0 {* E( B( y, r2 qNRI+           0.022727273 0.021564394 -0.019801980 0.0657894741 X" S3 c6 z% A: g
    NRI-          -0.020833333 0.017312438 -0.058823529 0.007518797
    : w0 `8 e  k0 X4 }+ A. f! T) @Pr(Up|Case)    0.034090909 0.019007629  0.000000000 0.072164948
    ! e# j  [+ l0 X* a9 \Pr(Down|Case)  0.011363636 0.010924271  0.000000000 0.039603960! Y4 j! i( [6 m* Z! ^- b
    Pr(Down|Ctrl)  0.013888889 0.009334685  0.000000000 0.035211268
    & X4 {. R1 J3 L( k$ }Pr(Up|Ctrl)    0.034722222 0.014716046  0.006993007 0.066176471
      X4 h/ C, j) l' U6 q# b* P
    2 h3 H7 b$ z/ q8 y3 G/ d12 p: J" ]9 K& ]& v
    首先是3个混淆矩阵,第一个是全体的,第2个是case(结局为1)组的,第3个是control(结局为2)组的,有了这3个矩阵,我们可以自己计算净重分类指数。
    1 `9 M) G2 K1 B1 M' ?, y- I: Y# w5 e9 K, @* P
    看case组:9 \3 [# W8 M8 F  L

    , ~# Q4 l8 n5 b9 |7 a净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273
    - i1 w/ ]+ }& E) a8 h+ g5 a- {( ?
    再看control组:+ c4 ~% K0 k/ M& _

    1 g0 t& P. n3 @净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333$ g) o5 w* {: f/ a: ]
    # U5 C! X8 v) n( _
    相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.000315657
    " J' y0 E& z1 _  M4 F
    5 Q: [0 C, G. {, J& K3 c, O3 y再往下是不做bootstrap时得到的估计值,其中NRI就是绝对净重分类指数,NRI+是case组的净重分类指数,NRI-是control组的净重分类指数(和我们计算的一样哦),最后是做了500次bootstrap后得到的估计值,并且有标准误和可信区间。
    # n9 _6 o3 ^. f6 H% _$ u; J0 A  X* _; L% R; r% j9 M
    最后还会得到一张图:! u/ j' g6 m! @+ x# H; p: _
    & ~, b: B- D  S4 Y
    这张图中的虚线对应的坐标,就是我们在cut中设置的阈值,这张图对应的是上面结果中的第一个混淆矩阵,反应的是总体的情况,case是结果为1的组,也就是发生结局的组,control是结果为0的组,也就是未发生结局的组。
    $ N2 N* y6 g) @3 R# I8 a2 k1 X9 u
    # N4 V* T* E7 w# h  K) vP值没有直接给出,但是可以自己计算。
    , `+ ~. c1 T% m
    # j9 q. z5 `0 @" h4 a0 T# 计算P值# s  [# @5 Y0 O6 }# T& O' e9 X
    z <- abs(0.001893939/0.027816095)
    1 I. |& K8 N+ Y* Z2 W1 M( G: J# @p <- (1 - pnorm(z))*2
    + m# C  I9 t& _p6 z5 X8 E) w- x$ a
    1
    7 l$ h$ _8 U& H: C## [1] 0.9457157
    7 H( Z$ Q3 h$ L& i9 X15 Y* R5 Y) H' Q6 O. K8 J! ?
    PredictABEL包) o) U. \4 I9 S9 D/ u, P
    #install.packages("PredictABEL") #安装R包
    / V* m  ?( H. l$ Flibrary(PredictABEL)  
    ) ~8 T9 G0 j# M4 a7 n
    # G. `% z1 n3 j5 [; s  R6 _# 取出模型预测概率,这个包只能用预测概率计算+ [% |: g3 x& P! P
    p.std = mstd$fitted.values
    9 i" T, N- G6 `# t5 p1 r5 Qp.new = mnew$fitted.values
    , y, [. ^  y7 S3 K, E9 Y* E& F17 l9 ^# F) M5 b. p" @
    然后就是计算NRI:
    3 W7 x% w; H1 y3 z& z( V2 u6 R3 o, z3 _, Q6 o9 r: S- p3 `2 L! U
    dat$event <- event( P9 k( R* h4 E. F3 p" R& |' G: J

    , \% S1 z, t1 Z' \, y$ preclassification(data = dat,
    1 X) G* p1 \/ w1 t) m# g$ B                 cOutcome = 21, # 结果变量在哪一列
    4 k  J" ?' T8 W- P                 predrisk1 = p.std,& T/ |6 y1 l& Y& A/ `
                     predrisk2 = p.new,) G5 G- I% v" G+ H: J5 v
                     cutoff = c(0,0.3,0.7,1)& ~7 H+ [  N9 W8 o, X" l
                     )3 B1 h# S) p9 j0 d0 b% J; w
    1
      y1 e% R* x! Q##  _________________________________________
    # a# l, I0 b' x6 U# q6 l# L##    Y( f& z# w4 D* B. `( b
    ##      Reclassification table   
    6 X* n) n0 p/ j: b( I##  _________________________________________" K& e" K7 h( q8 u! d) x
    ## - J- U9 q; g; C) Z
    ##  Outcome: absent 1 E( M" S* c6 ~6 d  @& f
    ##   
    ! g; K1 h2 Q0 L) o##              Updated Model
    % G+ u6 X- X7 p' o7 Z## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified) E  R7 i4 Z  Q3 |( `
    ##     [0,0.3)       121         4       0               3
    : k/ L, E1 O$ Q& R( v4 {9 f##     [0.3,0.7)       1        13       1              13
    - l2 @/ G* }, ], ?  M##     [0.7,1]         0         1       3              258 Z) l: y  q% ~  N' _7 A% y
    ## & o" c2 r5 V) i" y% H* g% M
    ##  6 ?% B0 s5 \+ x1 ?7 @3 i8 Q
    ##  Outcome: present 0 \4 E, g. Z0 M: S
    ##   3 R1 F* Z' ]3 t  x9 W0 _( P9 M
    ##              Updated Model
    . ?  C; d1 y" u2 {% c## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
    2 o# C, w& P, A( L: ?/ w##     [0,0.3)        14         0       0               0# V& B* y. q  l& g  G: W
    ##     [0.3,0.7)       0        18       3              14) N: X2 S6 x: ^# V
    ##     [0.7,1]         0         1      52               20 A' Y# f6 L2 Q4 r
    ## 9 T2 J/ u+ M& J$ }$ `
    ##  
    $ R( }2 l8 F2 I0 w##  Combined Data - t# `* n2 Y  I: F- |
    ##   
    6 d3 S3 J+ H* m. L3 l+ t##              Updated Model
    5 l; w8 `/ ^% d2 g$ o2 s/ z8 S+ h## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified+ _- \5 `  b  W( a" l' W
    ##     [0,0.3)       135         4       0               3
    " T; t- K3 t6 p##     [0.3,0.7)       1        31       4              14
    & \+ C4 G2 z! _##     [0.7,1]         0         2      55               4" R) ^  Q$ [- Z: [7 `; D/ t
    ##  _________________________________________
    ; @+ x& M0 F& P) x, V2 L" Q## . M  K4 P) a/ w& R
    ##  NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806 ( a* }1 K7 d& A
    ##  NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048 ( B* c  M& }" S- W6 o  f7 |+ @
    ##  IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.28396
    : [! D7 q; e6 }" L& K! f3 e! Q+ u) R4 u& y
    1  ~' W8 j2 ?4 f$ @; |$ f
    结果得到的是相加净重分类指数,还给出了IDI和P值。两个包算是各有优劣吧,大家可以自由选择。: L7 U7 T( l1 p. U4 y

      b; D! I- _$ p6 Y8 l生存分析的NRI  x/ X! A4 L4 p% V0 Y$ c
    还是使用survival包中的pbc数据集用于演示,这次要构建cox回归模型,因此我们要使用time这一列了。
    $ c* F# x1 G* I9 ~. v. Q& i7 d- ]# l# l
    nricens包+ q4 G: ^; K: ]8 p
    library(nricens)
    + [3 j+ N& ]5 I3 X: G+ E/ Glibrary(survival)( a4 U, ^5 `/ q4 ]) q. X; s! v  l

    * p( ^1 c8 s, d) v9 O1 jdat <- pbc[1:312,]
    . x2 g5 q  T! |! N+ \5 Bdat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡
    * r& t2 y$ p! P% |( d6 n7 J1
    & N3 i. ]) L+ P) X然后准备所需参数:% k: S' i% o; X4 d) N9 Y* z& K

    & B$ b' c' ^0 f- h, j7 {# 两个只由预测变量组成的矩阵! |( j% o  w5 S% y$ _+ B2 f
    z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
    * A9 i! m- h. B" H7 U6 iz.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))& V0 `6 N0 M1 N0 |
    7 k% H: _. _- {0 I/ q
    # 建立2个cox模型0 a4 a; v1 b, X+ T0 c" h" q
    mstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE)9 t$ B* z& s; L; r  R4 m1 U
    mnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE)
    8 P! _4 V1 H1 \. g+ Y" B$ B1 m: a' Z$ c
    # 计算在2000天的模型预测概率,这一步不要也行,看你使用哪些参数: h4 f& X5 P* H( w' \. ~% m
    p.std <- get.risk.coxph(mstd, t0=2000)
    ( F4 z$ z7 P/ K( ^p.new <- get.risk.coxph(mnew, t0=2000)  H% V# v+ o2 D6 X5 \
    1
    , }; T9 h2 v1 N7 j计算NRI:" j" ^# L/ Y. M5 {; `$ p

    " {: @+ ^+ ^. S  F' R$ f. W% a9 ]& fnricens(mdl.std= mstd, mdl.new = mnew,
    : L7 U- q! |% ]2 t2 X/ x        t0 = 2000,
    # \7 R* G8 U) c6 o( A& ^& r2 ^5 G        cut = c(0.3, 0.7),7 r6 \( r. O- M
            niter = 1000, " [6 d# O: i0 c: Z
            updown = 'category')3 r( v# `) `4 Z1 O( d

    4 Y8 t) n" I2 i4 a1 L, D! dUP and DOWN calculation:" P1 H4 J% z8 e4 z0 H: B$ I
      #of total, case, and control subjects at t0:  312 88 1446 d/ B- s8 M0 T; K

    5 K# C3 m5 s" _) K% q  Reclassification Table for all subjects:
    ) g. B. S/ Z! j3 C' l% Y0 |        New
    6 V4 q6 l1 B1 E! g. UStandard < 0.3 < 0.7 >= 0.7
    . Q: K1 R8 e  Q7 Q% _" E9 d  < 0.3    202     7      0
    7 u8 c2 Z! Q& g9 x. I5 z  < 0.7     13    53      6
    ( D& H& b' u3 q+ `6 r  >= 0.7     0     0     31: ?3 r$ O/ G( a2 u- i
    ! L2 B9 i  u8 E4 i: W, H
      Reclassification Table for case:
    1 }+ a- Q5 x1 R: A' _1 ]        New
    + l, X% E  ^9 c. NStandard < 0.3 < 0.7 >= 0.7/ j0 A( j  f7 B
      < 0.3     19     3      00 E' |4 i" ~4 m/ I* E% Z6 N8 T
      < 0.7      3    32      4, |2 O% V$ D1 H) j
      >= 0.7     0     0     27# i. K# }( @3 f4 }* Q: @
    * V. w! q$ W7 |- h, ~1 }; u, V. \2 Q
      Reclassification Table for control:/ ^8 d' K* F5 c) m- r
            New
    8 D. n9 q9 C- Q  l# n! L' k* [Standard < 0.3 < 0.7 >= 0.7
    ! x5 ^1 C3 h$ T3 B$ ~" d: w' H  < 0.3    126     3      08 z  Z$ V6 i/ h! o3 u
      < 0.7      5     7      2
    8 G: g& ^& U$ O# C  >= 0.7     0     0      1
    ) m+ |% Q7 _- N2 Q, |! X+ _; G! j4 ^- ?  i  u( T
    NRI estimation by KM estimator:
    ! q, S3 X; G' C  V7 o- ]
    + g6 `# X) v% }; _3 Y3 B5 ^Point estimates:  ~/ h. l+ w( a7 o
                    Estimate9 o! \8 H( u  O, ]8 ]# d- z
    NRI           0.05377635
    1 M4 c+ M# I1 L7 n3 u, {: F, JNRI+          0.03748660. u1 z+ Z' |3 j3 m5 ?' \& U
    NRI-          0.016289746 N9 u0 F2 S. H, Y) B9 b9 G
    Pr(Up|Case)   0.07708938: X$ {) G1 j) d, g* S5 l0 ?
    Pr(Down|Case) 0.03960278: L& q* h/ n5 b! E, n" O( ^: P2 O
    Pr(Down|Ctrl) 0.04256352
    : z' r+ P9 z( M' L5 TPr(Up|Ctrl)   0.02627378
    # r% h+ ^( W( X4 X
    ; q8 N" x3 y# x; h% hNow in bootstrap..4 v7 |) M# G/ Q) t# `- t$ Q

    ' G# u1 r+ S* y* PPoint & Interval estimates:
    " h7 A* j6 z3 a# _! ]5 y                Estimate        Lower      Upper
    " v8 P$ `6 |& O5 gNRI           0.05377635 -0.082230381 0.16058172
    / j! A& ]. g! [NRI+          0.03748660 -0.084245197 0.13231776
    . t' r9 i" D( @, d/ @1 iNRI-          0.01628974 -0.030861213 0.067536168 f: y2 w; o  c
    Pr(Up|Case)   0.07708938  0.000000000 0.191022911 b$ A2 A/ Q; `" f9 p
    Pr(Down|Case) 0.03960278  0.000000000 0.15236016+ F% H: v. h. e
    Pr(Down|Ctrl) 0.04256352  0.004671535 0.09863170
    * s- \1 {0 F* {Pr(Up|Ctrl)   0.02627378  0.006400463 0.05998424
    ; _  ^9 a" _# S  Y& U
    % ]( @* m! z& T" X+ r1- I# ^0 {! A  i7 {' ^6 |$ h

    ( D( m* S& {% {) C$ {; N( b; \Snipaste_2022-05-20_21-49-38
    5 W! Z, l1 S" C4 k. N7 `7 P结果的解读和logistic的一模一样。0 G$ {* S3 {% ]  p2 I, Q

    ( M$ O( I& @- A; S* ssurvNRI包( F. L' I; F; L6 U' g) v
    # 安装R包. P* Q3 N. i: d: k+ y) t% V
    devtools::install_github("mdbrown/survNRI")) W9 E. z7 _, t7 ], C
    14 t8 k! C, x! W7 \# L/ i
    加载R包并使用,还是用上面的pbc数据集。
    * j) Y! A0 b" L- \' c
    0 m3 \+ P$ c( H" M! [( mlibrary(survNRI)
    0 @+ e& q0 N/ ~1 Q" p1, h7 W0 n7 q& Y1 ~# q$ G: d
    ## Loading required package: MASS, q. [/ h# B. J% J
    1
    $ f" G3 A9 P6 l" k) y" zlibrary(survival)+ m4 m1 ~% m; a% M& N/ r' \! d: z

    , O5 {; d+ p% j! C5 P) D8 w+ X5 [5 f# 使用部分数据
    * e5 ?2 z; V: ~2 G& gdat <- pbc[1:312,]
      v5 P* Z; b& S9 Y7 ?6 r7 Vdat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡1 q: q! C3 Q# g6 c

    . N5 J; w  E) p; J# G7 Eres <- survNRI(time  = "time", event = "status", - u9 n* d1 T# ?( m
            model1 = c("age", "bili", "albumin"), # 模型1的自变量* w& Y" N9 I  [4 _4 o9 D
            model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量
    " K" `7 H  o# k( w7 S. p3 ^6 t1 u        data = dat,
    1 X5 y3 z3 ^5 p% u$ ^$ ?% x9 m$ r6 B        predict.time = 2000, # 预测的时间点5 u! M) X0 w6 V  D2 f
            method = "all",
    5 F- i+ i- U( i! C& m8 P& p+ B        bootMethod = "normal",  
    2 J7 m! m' \! p+ N( f3 C        bootstraps = 500, 1 g1 k3 @: d# O# J/ x: J# [% P
            alpha = .05)
    1 R% g. F/ B1 b, v$ A0 m4 F
    2 O+ `" ~, m6 E% W5 H8 E1
    # p/ y, N% B" R' }/ i/ m查看结果,$estimates给出了不同组的NRI以及总的NRI,包括了使用不同方法(KM/IPW/SmoothIPW/SEM/Combined)得到的结果;$CI给出了可信区间。( `, @9 U; v$ b! X0 {) h

    # f- o7 W0 ]/ f1 X3 {res' H/ r% q) z' q0 r' y5 V: Y' W
    1
    1 ^6 j8 S- o7 E# l* O9 h## $estimates
    % \1 f" D2 g2 c! \) U  `  g0 E" F##            NRI.event NRI.nonevent       NRI  G+ B9 i  g7 o3 x
    ## KM        0.20445422    0.3187408 0.5231951
    + w! e! f, S! j; `/ s, s$ W* S## IPW       0.22424434    0.3273544 0.55159875 P# x; q/ S- d9 v2 c* q3 `
    ## SmoothIPW 0.19645006    0.3144263 0.5108763
    1 w2 u) L) {; M% i% W## SEM       0.07478611    0.2632127 0.3379988
    & D4 T' Y' C, @## Combined  0.19633867    0.3143794 0.51071817 @# `, l: I3 S& Q% r& j
    ## 8 y  [$ e8 I( c6 R! T( j
    ## $CI
    4 @. q' B/ H  |; J## $CI$NRI.event" z0 e! ^% ^7 o+ Z+ W4 d8 z
    ##                     KM         IPW   SmoothIPW        SEM   Combined* K) r; ], a+ a5 k
    ## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.04737239 |7 w9 m* F8 }, b) M3 Q0 Z( I6 K
    ## upperbound  0.44806768  0.47033936  0.44014214  0.2658309  0.4400496
    ! `' T5 r, w$ n  \! O##
    1 K1 H4 }% V& p& M2 T## $CI$NRI.nonevent
    / Q) H- d5 r; N  \3 e##                   KM       IPW SmoothIPW        SEM  Combined2 \9 u/ |% Z% W
    ## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.1286426
    1 e% |( w( q( z# y2 ~8 Z9 r## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.69645490 m1 ~$ W/ Z) U7 ]% |) q0 ^: K# I
    ## 9 j+ E( w! y% ?) Y0 q" R
    ## $CI$NRI
    3 L( T% F! _9 T2 g) q% V$ n##                     KM         IPW   SmoothIPW         SEM    Combined
    8 m8 Q4 Y8 {, N9 P' @; [7 S## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409& W) d# o% y0 B- I: P  f
    ## upperbound  0.89306122  0.92464359  0.87970125  0.64253510  0.87953153' f% d8 P/ l6 R
    ## : j- Q$ s- V/ {  l; D) j
    ## - r' f- X( c7 I
    ## $bootMethod
    6 l2 V: L3 a( b* f, j+ M' N5 {& _## [1] "normal"
    ; h8 j# u3 w! S# z# |## 2 o: g! U. Z0 Q/ Q
    ## $predict.time  k0 M" W1 j  ?
    ## [1] 2000
    , m+ d1 [* }( _6 H## 0 y/ R% s( E% x6 J$ m
    ## $alpha
    # v5 b- N" F1 U3 j3 R## [1] 0.05
    7 [# d9 n* Q% i/ |## ! L4 ^2 X( Q& J$ i& N  M2 t0 b/ G. x* Z
    ## attr(,"class")
    . b* P. D% v8 [! \## [1] "survNRI"
    9 u; w" J9 p( Y3 R, d4 S) d9 h) ?2 P) D
    1% @, E  t# U3 E' Q2 H, F0 u8 z9 [7 ~
    OK,这就是NRI的计算,除此之外,随机森林、决策树、lasso回归、SVM等,这些模型,都是可以计算的NRI的,后面会继续介绍。大家如果有问题欢迎在评论区留言。* E& x! i1 S' @% A7 y
    . @: r* ]3 H6 `
    本文首发于公众号:医学和生信笔记
    . R" }% y( u, n0 H& e3 k: p* z: r9 H. }( x7 C$ h
    “ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
    % f& F$ H- y# ~; O3 `本文由 mdnice 多平台发布/ q! i% t. @5 @' [! {8 H- o; r
    ————————————————, Y. `- d1 D; D1 @. d
    版权声明:本文为CSDN博主「阿越就是我」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。: e- Z4 U; B& h+ @% r9 F5 Y( Q
    原文链接:https://blog.csdn.net/Ayue0616/article/details/126768006
    9 C  Q4 k; O# l9 i
    : d! d4 ]: i# {& ~
    + L4 {4 s( d; `% v. Q9 Q
    zan
    转播转播0 分享淘帖0 分享分享0 收藏收藏0 支持支持0 反对反对0 微信微信
    您需要登录后才可以回帖 登录 | 注册地址

    qq
    收缩
    • 电话咨询

    • 04714969085
    fastpost

    关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

    手机版|Archiver| |繁體中文 手机客户端  

    蒙公网安备 15010502000194号

    Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

    GMT+8, 2026-6-12 09:24 , Processed in 0.459884 second(s), 50 queries .

    回顶部