QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 3056|回复: 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

    9 g" w1 R# W8 i净重新分类指数NRI的计算
    ( J( s" D; T; B7 O# Z“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。) H! i0 l' _8 [8 c8 I
    NRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!
    ( c! H5 a. a) q5 O' d" P
    ; t/ W  r. n$ i" U4 q0 C4 ]在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens包,PredictABEL可以计算logistic模型的净重分类指数,survNRI可以计算cox模型的净重分类指数。
    5 J+ S% i2 Z. R) T% ~- v, O& v2 @6 u9 w
    logistic的NRI
    0 f% W5 Y# \* Y1 ?+ f' h9 M) ~, q2 C9 Cnricens包
    + P; Z# _4 h! @4 iPredictABEL包
    $ [- q2 p& Z. l生存分析的NRI
    1 k$ Z3 x4 P. q( P9 t) z! Nnricens包: A7 }/ u( J% \# K$ S3 U
    survNRI包$ m% d8 `' c* F
    logistic的NRI1 h( p" k: F' ]6 R% B
    nricens包! e! L; z1 Z: p: K9 _. W/ X
    #install.packages("nricens") # 安装R包- P3 m4 h( Q& b0 `# P
    library(nricens)
    4 o- g/ X! V/ s: o6 Q1  U% ?( u$ B+ }" P, M8 s
    ## Loading required package: survival
    6 b/ H1 D' S- e& z  N: S1
    $ l0 Y) }. I& j+ d使用survival包中的pbc数据集用于演示,这是一份关于原发性硬化性胆管炎的数据,其实是一份用于生存分析的数据,是有时间变量的,但是这里我们用于演示logistic回归,只要不使用time这一列就可以了。
    0 W' F: j/ `& A% m' W
    " H2 `0 V& L7 U% [) |library(survival)0 m- E1 l7 }% X; A/ b3 S: J

      a' ?1 c. F& I  q+ b  B# 只使用部分数据
    & t# u, n: R6 c& N2 o; @dat = pbc[1:312,]
    3 {) A! F% N4 V9 f9 xdat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]) a. h! f" h1 U3 E; K/ p+ ~8 I5 E

    6 S" J9 S) u$ c. Mstr(dat) # 数据长这样( }3 `* {1 |6 B5 p  V
    1) z9 h- B5 z$ S
    ## 'data.frame': 232 obs. of  20 variables:# G4 T: `/ H' Q( @5 P9 |1 C( H! P- Z
    ##  $ id      : int  1 2 3 4 6 8 9 10 11 12 ...- e' U, I; J/ K8 Z8 X( P
    ##  $ time    : int  400 4500 1012 1925 2503 2466 2400 51 3762 304 ...
    ; S* ]+ \4 R& M, R4 W( _##  $ status  : int  2 0 2 2 2 2 2 2 2 2 ...
    9 w# K9 W* R' l( q$ C##  $ trt     : int  1 1 1 1 2 2 1 2 2 2 ...
    ' K6 @  i) S) B; p$ }2 k+ T4 s##  $ age     : num  58.8 56.4 70.1 54.7 66.3 ...
    ' g1 J4 R2 K+ X##  $ sex     : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
    ) W$ ^# l6 I6 B& e' q##  $ ascites : int  1 0 0 0 0 0 0 1 0 0 ...& c! L3 p8 h5 E; `
    ##  $ hepato  : int  1 1 0 1 1 0 0 0 1 0 ...- y2 Z8 e" E, W8 {4 H7 t. j
    ##  $ spiders : int  1 1 0 1 0 0 1 1 1 1 ...) G# `8 v* U5 S- A. Z4 \8 m
    ##  $ edema   : num  1 0 0.5 0.5 0 0 0 1 0 0 ...9 l# C. X5 O! C2 u  A
    ##  $ bili    : num  14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...
    8 ?+ S/ i, I  ~) a( J& s1 |0 J- ]##  $ chol    : int  261 302 176 244 248 280 562 200 259 236 ...: [6 C: M1 a2 }
    ##  $ albumin : num  2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...
    ) N$ v- r1 k9 |9 e##  $ copper  : int  156 54 210 64 50 52 79 140 46 94 ...0 b" \# R" k% s" ^
    ##  $ alk.phos: num  1718 7395 516 6122 944 ...
    # y4 {& Z9 Y8 I/ F##  $ ast     : num  137.9 113.5 96.1 60.6 93 ...
    8 P5 i1 K0 H% o* U##  $ trig    : int  172 88 55 92 63 189 88 143 79 95 ...
    - u' Q& L0 u2 i) i##  $ platelet: int  190 221 151 183 NA 373 251 302 258 71 ...
    0 I; `: w  [8 I3 L+ Y; I##  $ protime : num  12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...7 v% a* i% B7 V* }
    ##  $ stage   : int  4 3 4 4 3 3 2 4 4 4 ...
    5 R! ^/ U2 ~0 r8 }; v+ R! x- a* s7 S
    13 o) d6 m3 {* G0 {& C# J. d
    dim(dat) # 232 200 y4 Z( _7 _! o/ X1 L* [
    1% K8 I$ K. ^6 E: C) x1 y
    ## [1] 232  20. y; T+ {3 ?% `
    1
    8 J8 g: g. u0 k/ F然后就是准备计算NRI所需要的各个参数。+ Q, y4 F$ t5 C. \) E! M6 O

    2 g7 @: C8 E/ D; n$ I# 定义结局事件,0是存活,1是死亡9 G" U- S+ h3 y8 W
    event = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)
    7 B7 @5 g& l3 q/ {  f/ r' l1 v
    # 两个只由预测变量组成的矩阵+ @3 y8 @7 r* X) Q5 U0 ?
    z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))% H5 S# s! V. c1 X# G, h, C
    z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))6 S" v0 g: n8 H9 S6 Z% N

    ! D% X! P, s/ d- _% U  `# 建立2个模型
    ! H  ?1 f" ~2 H: S  M2 `8 Dmstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)
    2 h0 m. }2 Z% w* r$ U# Smnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)" f$ o: c1 H, [5 L/ x
    0 i# S) t* g' _. E. k: q
    # 取出模型预测概率: L' v/ ?/ d1 U5 B7 \% Z) ?2 v
    p.std = mstd$fitted.values0 d& }2 w& ~* g$ [5 k+ _. _( o% L
    p.new = mnew$fitted.values) X0 \3 W# Y/ H" @9 t9 ^% [

    : e, Y# z1 s9 R' c) j18 D' i2 T3 Q3 |* W" r& Q
    然后就是计算NRI,对于二分类变量,使用nribin()函数,这个函数提供了3种参数使用组合,任选一种都可以计算出来(结果一样),以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。' v3 {+ E* z8 S# y, s: n

    : a4 u& ]4 Z% S1 K# 这3种方法算出来都是一样的结果
    ; \3 \1 J! a6 ?+ x' V" A1 ?& W
    " Z+ M0 r% j0 L+ p- \- I7 `, R# 两个模型7 P6 g' S' V+ S
    nribin(mdl.std = mstd, mdl.new = mnew, & k( W# `! H. j/ y! p  Z; r: g) ?
           cut = c(0.3,0.7), : E1 z5 I8 f, H. ^6 `7 G7 a
           niter = 500,
    ) i; }' s  v1 }( u       updown = 'category')
    ( C  D# _5 Q, u$ H8 a- Q( e3 t9 g% F& a! [! g
    # 结果变量 + 两个只有预测变量的矩阵
    " m" I0 v  X9 E7 D2 Unribin(event = event, z.std = z.std, z.new = z.new, 2 x' u$ t2 Y2 W, S
           cut = c(0.3,0.7),
    # \# C, i# d/ O/ S       niter = 500, " @0 C, f9 h. u! K
           updown = 'category')9 d, X+ h5 B$ ~7 x0 C! O
    ' t. n9 @+ ^4 _# F3 L$ ?9 V$ {
    ## 结果变量 + 两个模型得到的预测概率  h: ?7 V; B5 z6 p4 |( e% |
    nribin(event = event, p.std = p.std, p.new = p.new, 4 {9 K) V( Q. z# P
           cut = c(0.3,0.7), $ s% Z; s% A8 a# x
           niter = 500, ( V: l4 U& I! m, L2 B
           updown = 'category')& M8 j* h9 K6 T' `

    ; q9 |* W. v/ A2 S8 c1
    0 L% [$ }* X0 B2 a: s1 r其中,cut是判断风险高低的阈值,我们使用了0.3,0.7,代表0-0.3是低风险,0.3-0.7是中风险,0.7-1是高风险,这个阈值是自己设置的,大家根据经验或者文献设置即可。5 y" t7 I% k9 N" K+ y1 H, H
    + W6 B1 S9 [8 ]. V/ U9 H
    niter是使用bootstrap法进行重抽样的次数,默认是1000,大家可以自己设置。
    $ p6 o( ~4 p- k( g6 ^+ ]6 g3 a, D. _
    + X0 m' G1 m0 d3 l0 V9 wupdown参数,当设置为category时,表示低、中、高风险这种方式;当设置为diff时,此时cut的取值只能设置1个,比如设置0.2,即表示当新模型预测的风险和旧模型相差20%时,认为是重新分类。4 S$ a9 ]# H) I2 q5 D
    7 G. y; f* \! w
    上面的代码运行后结果是这样的:: _! e' m, q. T( n$ ~0 n, Q

    . v# \0 b0 a- M( I# T, EUP and DOWN calculation:) Y4 m% Q, |1 _% N
      #of total, case, and control subjects at t0:  232 88 144
      `) q, i. Z% W! o& Z( f; ~9 c) F3 G6 _
      Reclassification Table for all subjects:! O5 O  F* j( h0 i  w) n
            New
    5 ]/ |3 `5 s* R, T& l+ ^Standard < 0.3 < 0.7 >= 0.7
    . L7 `0 [7 c2 P, x# l5 h4 N  < 0.3    135     4      0  E6 f+ T! l- \# u
      < 0.7      1    31      4
    # R8 E7 p9 c5 O  >= 0.7     0     2     558 b( F0 f; l0 M; c
    ( I6 M& L9 F/ s, G- a" G
      Reclassification Table for case:& Y; s/ @, |* f% ]* e9 X% ~1 r/ I
            New
    " z7 {" p7 K2 J0 [Standard < 0.3 < 0.7 >= 0.7
    + i! d& r3 v* a1 R  < 0.3     14     0      04 m: U4 ?0 w- D6 f8 g* }3 b. D' N
      < 0.7      0    18      33 O: _# r6 }+ @) T! L
      >= 0.7     0     1     52
    & t7 D+ W- N( e. K' Z9 n/ z& D3 }3 T% N# M( q' V
      Reclassification Table for control:5 u: A, u/ x# K8 ]0 C
            New/ k5 v) [; P6 V3 N
    Standard < 0.3 < 0.7 >= 0.7
    ; h7 s( g7 H/ }) g+ Q% ]  < 0.3    121     4      0
    , e' j! l5 H& _9 ?  < 0.7      1    13      1
    " y5 X- S; }$ Y0 j  >= 0.7     0     1      3
    0 q1 z8 m  t$ N* ?, F  `  ~- H7 f" R8 @+ b; E
    NRI estimation:
    ' l& p: W6 Y% V. YPoint estimates:# b7 S) @, H' e
                      Estimate8 v; P' b6 n" |6 B5 {
    NRI            0.0018939390 |) B! F: n9 n4 L% Q
    NRI+           0.022727273! Q3 q9 u# N3 S. ?; F
    NRI-          -0.020833333
    * X% i5 N7 U" n8 _- {Pr(Up|Case)    0.034090909
    ) X! d. y! @6 q( oPr(Down|Case)  0.011363636/ {  U' z/ v1 D5 [  h4 w
    Pr(Down|Ctrl)  0.0138888898 R6 L6 h0 A8 m6 e2 C3 {5 Z7 |; D
    Pr(Up|Ctrl)    0.034722222
      G6 N7 M/ f2 r  Z7 E8 J5 l- ]1 ?5 x( a* Y0 t2 Q
    Now in bootstrap..5 ?, S3 p4 x2 D( a7 R; }/ x0 G( u( W
    : n7 z0 }2 O; C  B
    Point & Interval estimates:
    5 |/ m/ ?; I. R! i# \/ d. m5 Q                  Estimate   Std.Error        Lower       Upper
    - \1 G4 R' V) s, _; yNRI            0.001893939 0.027816095 -0.053995513 0.0553544494 {! W- i& ?: k" ?+ C" b8 j9 O
    NRI+           0.022727273 0.021564394 -0.019801980 0.065789474' r* M" J# u6 b- W
    NRI-          -0.020833333 0.017312438 -0.058823529 0.007518797
    ; t" x- g! [8 m  |Pr(Up|Case)    0.034090909 0.019007629  0.000000000 0.072164948
    * {0 Z/ P5 ^2 Y9 p: `2 ^2 VPr(Down|Case)  0.011363636 0.010924271  0.000000000 0.039603960! Y2 J0 @+ Z4 n' j0 h8 r
    Pr(Down|Ctrl)  0.013888889 0.009334685  0.000000000 0.035211268
    / Q3 d4 S  Y7 n( FPr(Up|Ctrl)    0.034722222 0.014716046  0.006993007 0.066176471& V( a7 L) F( c- V5 ?
    # f8 y8 P" Q3 U  g3 d  @+ W
    1( h& R4 Q( z& V0 g) g0 {0 X; T
    首先是3个混淆矩阵,第一个是全体的,第2个是case(结局为1)组的,第3个是control(结局为2)组的,有了这3个矩阵,我们可以自己计算净重分类指数。
    , R4 E4 K5 m2 P/ p+ T
    ; p6 r/ `; r3 r4 b8 ?看case组:
    4 _: G, b3 ]8 [6 |1 k( m0 M  T  f- k( g; u, U4 N" u
    净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273
      Q1 G! ^. Z) G/ H, P& m' ~$ t: n7 N: h! e7 A
    再看control组:
    1 Q  R2 }8 _: L7 b4 \  {. T4 d/ O, }+ W( \! E' L
    净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333
    + _! N1 R9 S# X3 F0 M) G& w$ O- _" Y1 X1 e$ D; Y% V  d
    相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.000315657+ ^3 {9 s) R0 ?: e4 K6 B- p5 R3 r

    6 [8 b9 ]" h- y+ I, C再往下是不做bootstrap时得到的估计值,其中NRI就是绝对净重分类指数,NRI+是case组的净重分类指数,NRI-是control组的净重分类指数(和我们计算的一样哦),最后是做了500次bootstrap后得到的估计值,并且有标准误和可信区间。
    + k2 x1 y6 n8 s( {. ~$ [1 t3 ?9 f! z* l: o  Z+ R
    最后还会得到一张图:* H# H/ L0 x4 Y/ y1 w2 ~1 Y
    6 q9 ^$ `2 ?" w  r
    这张图中的虚线对应的坐标,就是我们在cut中设置的阈值,这张图对应的是上面结果中的第一个混淆矩阵,反应的是总体的情况,case是结果为1的组,也就是发生结局的组,control是结果为0的组,也就是未发生结局的组。
    . |5 n/ ?0 K# c, L, [: G
    ( b4 o; ^) E: @6 C9 U% vP值没有直接给出,但是可以自己计算。7 F# ]9 U+ a. d% }. `
    4 H  B; _8 {+ L9 G; \/ O# [: O
    # 计算P值
    * S( A7 H; S' ~" n: `2 c7 S( j! sz <- abs(0.001893939/0.027816095)
    6 z1 q/ K" H( ]# k: ep <- (1 - pnorm(z))*2$ g3 z3 e( G, ~( U% _3 K
    p9 Q4 `" G5 P9 @, ~' s
    1
    ' @- v. E& w0 s## [1] 0.9457157
    6 G/ q4 Z" |4 m1 W  c* l15 L, C/ L* j/ x: N, W1 t
    PredictABEL包. J4 \5 k- p6 P) X7 }" P
    #install.packages("PredictABEL") #安装R包$ g$ Z* f, R: d0 r3 W
    library(PredictABEL)    A, x0 y4 E$ G2 X: U
    + _" m) k# K# z, M# P4 \
    # 取出模型预测概率,这个包只能用预测概率计算* N/ ~$ A) b, @- I8 p7 _  a% Y3 X
    p.std = mstd$fitted.values4 G$ Z& N' w& R/ Y$ y3 O
    p.new = mnew$fitted.values
    7 t9 @; _' I  s0 l/ V. u1& Q0 c& F  X2 _8 N, |- ]
    然后就是计算NRI:3 u/ V+ ]: l2 r

    ( o: O& x; C* J* U3 t* Fdat$event <- event
    5 b. p* X! G. n9 v
    5 L4 \' e: o5 c* {/ K" K/ qreclassification(data = dat,
    9 h4 r1 Z4 a* H& c9 u5 A# s                 cOutcome = 21, # 结果变量在哪一列
    6 v- U6 v, m: {. T: p: u" p6 @                 predrisk1 = p.std,
    & I( t5 m6 n0 R! s% ]  F) d6 D/ t                 predrisk2 = p.new,
    0 b8 ~: G- u9 \; D                 cutoff = c(0,0.3,0.7,1); \  o$ r" J6 J; z3 S
                     )
    7 b. R; z9 D0 G$ T! Y' c% c+ u; [5 o# {1% ^4 k. S- l6 T. `
    ##  _________________________________________
    / v& Y+ g5 e( \! W0 s8 |: @##  # e% o1 V6 z' m" k# Y* j
    ##      Reclassification table    ! k, I& v# B% Y* v1 f2 q2 M
    ##  _________________________________________
    . T  F  G& x4 Q( R##
    4 v. e: r/ W7 k# \# O) e9 B8 o1 t##  Outcome: absent
      p: w0 g; _/ l  ?0 X! J# w  ]$ m##   
    0 h9 R# d3 w' n) d3 }8 z: L##              Updated Model
    7 q% C! ?. E: e6 p! T- }## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified$ a% x3 Q& H/ B+ U
    ##     [0,0.3)       121         4       0               3
    8 e2 s7 F6 m+ j. B6 A4 _3 @##     [0.3,0.7)       1        13       1              13
    9 w9 N$ C5 _7 s! Z& w##     [0.7,1]         0         1       3              25( D7 V& o5 \/ l) d/ W2 ]- a; R+ r
    ##
    * |2 g5 W3 j9 L##  
    - c# `- x0 I' M+ r- r0 Q2 B##  Outcome: present
      @- A& J# G* R. f" H: K##   $ L4 S) B7 a; s  ?+ |
    ##              Updated Model
    - l! p- t: h! c9 l6 Q; n## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
    , ?# [" ^. j- t0 R' J3 W2 P4 K##     [0,0.3)        14         0       0               05 h/ W) N" j! P6 Z7 J
    ##     [0.3,0.7)       0        18       3              146 D6 V& |* k+ W& E& q
    ##     [0.7,1]         0         1      52               2
    % k- ]: X+ y! |* \" q##
    : C, A! Y0 h% X1 Q##  2 M* C. k, Z' W' o
    ##  Combined Data
    ' {+ {9 N5 h9 F/ z3 ]2 Y: ~##   / U8 I" {# x$ E& Q5 d* h" P
    ##              Updated Model
    3 X  q9 n8 r$ T9 j- I! D5 z% {# M2 u## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
    ) }. N- K6 I. z9 G% {9 V) k##     [0,0.3)       135         4       0               3
    ! I/ w3 g5 W% e7 g- n" F! T9 `##     [0.3,0.7)       1        31       4              14
    - e7 [0 b1 S% F9 _0 [. r##     [0.7,1]         0         2      55               4
      q$ E  x/ S% g4 _7 x3 p% ]##  _________________________________________5 h) k$ V( @8 r0 l7 u' L0 f/ ^$ }' g) u
    ## ! I- s' |3 {* i  ]6 r, Z% {* h
    ##  NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806
    # b8 `" x6 [" j5 r6 a2 B( P; G##  NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048 3 x* k3 }7 N$ F+ V
    ##  IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.28396
    2 W- H% |' _% M3 y1 m6 [0 d' K0 e: f6 e6 k5 K6 I5 [
    1  A5 A$ Q3 j. [4 c: Z
    结果得到的是相加净重分类指数,还给出了IDI和P值。两个包算是各有优劣吧,大家可以自由选择。1 S) P9 ]  `5 ?: D+ |

    5 A9 L* k  T; m& i8 M. ]' T  h/ o生存分析的NRI
    * P, J1 }( @6 I; \( Y' N" \" O+ y: L还是使用survival包中的pbc数据集用于演示,这次要构建cox回归模型,因此我们要使用time这一列了。
    ; H9 e: W' `: U: {! D9 q1 E0 I. b3 F+ A  m4 L. l* G( s. t
    nricens包
    , z% f! F* z+ o4 m8 Y( ?library(nricens)7 n3 r! G# H1 K# l
    library(survival)% Z# {, c5 T8 b; e. f, [

    ! _* n- \; }3 I/ \dat <- pbc[1:312,]
    0 F7 [6 s3 U- O% m$ X- ]dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡
    5 u7 G" {# M2 o; {! L' W0 l: S1
    % Q2 ~1 I# v' x然后准备所需参数:4 l5 g1 i7 v& p
    ( ?9 L2 M( K& |( y
    # 两个只由预测变量组成的矩阵
    + j, i, ~' _+ t  G+ I( w0 T; U  j8 tz.std = as.matrix(subset(dat, select = c(age, bili, albumin))); j  u  y* n- o4 V; l. p
    z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))! c7 y1 X5 P) J: z
    $ {( P* o( Z, B
    # 建立2个cox模型
    ; t3 S- A! m4 C  {mstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE)) x3 W+ q7 M: s8 S, L
    mnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE)
    ' A1 H6 k. ^1 F6 K( O
    $ q; x2 e3 N, w# 计算在2000天的模型预测概率,这一步不要也行,看你使用哪些参数
    7 _% ^/ _$ ]* N/ f; up.std <- get.risk.coxph(mstd, t0=2000)
    ! }% K" u" m5 ?; [2 ^/ mp.new <- get.risk.coxph(mnew, t0=2000)
    : }) M! U: n) C+ h1
    2 R5 d, _* b- w' l& c  l& S计算NRI:
    8 [3 }- _. E2 n( @
    9 |- H3 ]7 t' W0 o1 C) Dnricens(mdl.std= mstd, mdl.new = mnew,
    , z9 M) M/ e! r* S" V        t0 = 2000,
    1 R! F. c* x8 B- a7 l        cut = c(0.3, 0.7),5 f5 p3 b  o/ @
            niter = 1000,
    - Q0 Z3 e5 `; }        updown = 'category')
    # k# g  q4 ]; w- F# o- s; M+ e  S1 D- M
    UP and DOWN calculation:1 a7 x. q0 e/ Y! K& d# V0 e
      #of total, case, and control subjects at t0:  312 88 144
    - l, y7 N4 T; m% G
    4 x5 q  _, y  J) i6 N0 w  Reclassification Table for all subjects:, t7 P7 C. r! [! X8 U. j. h
            New
    / J2 }, N! g( ^  b* y' HStandard < 0.3 < 0.7 >= 0.77 b, k( P7 x2 {0 S* u! B
      < 0.3    202     7      0# ^6 K, s9 a& \, g1 `
      < 0.7     13    53      6+ i5 g; u) p8 I. w
      >= 0.7     0     0     31* m, }1 h- @- P9 X" D8 O
    . j$ @5 {! h4 @# S" I& s1 O9 \8 a, O
      Reclassification Table for case:3 ~0 h! p; l7 E8 Q1 O  F  D3 k6 N7 y
            New
    $ D2 ?$ d- \; Q0 q% @  j6 x1 MStandard < 0.3 < 0.7 >= 0.7
    9 h7 W; s+ `+ I# ?  < 0.3     19     3      0
    5 k/ D, Z  o3 A  n. Z  < 0.7      3    32      4
    6 N5 V, b% |1 u  >= 0.7     0     0     27
    ( v, r  x) n) X+ G
    2 a5 ~# T. [  P, F& }& c- n* i  Reclassification Table for control:! X) N! s( S! ~0 [% g
            New8 ~+ U: b8 m, }( w* g5 |( `! z
    Standard < 0.3 < 0.7 >= 0.79 w$ S% {9 k; ]1 L/ |+ {/ R
      < 0.3    126     3      0
    2 a! r( \; U7 y8 y- B6 Y  < 0.7      5     7      2
    5 j  g8 f' x$ S$ }4 A) m& @  >= 0.7     0     0      1
    - I0 d" \1 ^$ w
    2 j9 k4 i' w" ^# M4 M$ jNRI estimation by KM estimator:* \6 v/ H- I; y. l

    ! b1 v# x7 }% Y, ?. [6 xPoint estimates:
    ; `$ ?7 h) R9 J+ n/ m                Estimate2 R8 E9 C7 A$ n, @. s7 i- y5 g
    NRI           0.053776355 }5 C4 _. r; j6 j8 M
    NRI+          0.03748660( _! Q6 `( o  ^5 L
    NRI-          0.01628974# m7 M& S7 i) `% J" f2 p
    Pr(Up|Case)   0.07708938
    % u& E0 E# j, u. M! E( T7 X. QPr(Down|Case) 0.039602788 u0 O, _6 n+ Q: S+ C
    Pr(Down|Ctrl) 0.04256352
    $ L: r" R) c( |8 ]Pr(Up|Ctrl)   0.02627378+ B: A5 s8 l# w- B; |( s7 t4 Z

      U. i9 n% K- }% x) V. u2 b/ YNow in bootstrap..
    0 G% c9 a7 U) I$ t
    . n1 ^$ I% _9 \4 f( v0 sPoint & Interval estimates:! `8 [% R1 _$ L" y( h9 m
                    Estimate        Lower      Upper  L5 F! b9 W" ^2 p9 a; n5 q
    NRI           0.05377635 -0.082230381 0.16058172
    % w' E* E6 p% D7 LNRI+          0.03748660 -0.084245197 0.13231776" m2 \: P% y  O  }
    NRI-          0.01628974 -0.030861213 0.06753616
    , r6 S+ K8 t, [2 W9 C! J* WPr(Up|Case)   0.07708938  0.000000000 0.19102291* t% e/ T3 d. }; h# X' n- C
    Pr(Down|Case) 0.03960278  0.000000000 0.15236016
    ; j6 t1 ~" N, J4 ~Pr(Down|Ctrl) 0.04256352  0.004671535 0.09863170( D) V' W, J: G: @6 G: W
    Pr(Up|Ctrl)   0.02627378  0.006400463 0.059984240 ]+ P/ _3 u  w. g8 Y

    2 }+ P8 k& Y- i4 A/ A8 l2 m) [( E1
    0 ~2 K2 H( V0 q1 D# Y+ O7 f+ [5 u8 P. E0 h0 K0 Q: l" Y
    Snipaste_2022-05-20_21-49-38# J# G& z6 F" |. w# X, @: f* x
    结果的解读和logistic的一模一样。: X) Q8 v+ K' S
    / R2 [9 D) J. }
    survNRI包
    5 A; H9 w/ j) O3 a- I# 安装R包
    / B. x9 W1 N# Z7 R, B) h! L& j$ \devtools::install_github("mdbrown/survNRI")9 o4 V5 H  h. A& F
    1
    8 J' E0 ]; f) T( V加载R包并使用,还是用上面的pbc数据集。
    7 }! g$ E/ A  {7 C
    % ^, N6 S6 R" h( @library(survNRI)
    ( ~4 X4 p1 u4 A" I11 v- V* O, q# x' Y2 }! B8 S6 k4 D* ^
    ## Loading required package: MASS- C. s$ T% a; i
    1
    $ c$ |% M5 D) n5 j, [! {/ ~& \$ Ulibrary(survival)
    , U9 c8 n, C* u7 h8 d4 i
    3 x- z1 f8 e6 i) r) ^# 使用部分数据
    6 C% F2 m* W7 t9 ]% K) ddat <- pbc[1:312,]( D4 Y$ d% B* Y1 G5 A- Q' r
    dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡
      h$ `: z. T$ G) N5 }% |; V4 j' x/ f/ ^
    res <- survNRI(time  = "time", event = "status",
    , ]+ q& R3 t/ M' }        model1 = c("age", "bili", "albumin"), # 模型1的自变量7 N" ]' R; A$ \  p
            model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量* P( U# {& N- Z6 B+ k
            data = dat,
    ; U" }7 y2 i& s        predict.time = 2000, # 预测的时间点! v9 [: j1 B& a8 B4 n% ]9 u
            method = "all", * A6 G! T7 M" X% C
            bootMethod = "normal",  
    * o  R, q! |, {7 \, ~2 @( Z1 b        bootstraps = 500, 3 q* {& ]0 ?* g, \7 I& y
            alpha = .05)9 T2 C: I: H' @, j0 [2 }

    0 x& T* q5 G$ x/ ^6 b1
    % _! W. x3 f" x, a查看结果,$estimates给出了不同组的NRI以及总的NRI,包括了使用不同方法(KM/IPW/SmoothIPW/SEM/Combined)得到的结果;$CI给出了可信区间。# H3 T* T3 O8 A" `! L# A

    1 Y2 s/ C! C9 w* x9 b$ y7 Bres
    ' D; n" [- V7 C  ~2 s1
    . S' R1 m1 S# ~6 F/ p" e2 `' j## $estimates
    3 {2 H+ b% D" k) G##            NRI.event NRI.nonevent       NRI9 m+ d3 l* `* y# A* _
    ## KM        0.20445422    0.3187408 0.5231951! e* r- V% m6 y* V& F  ~& N
    ## IPW       0.22424434    0.3273544 0.5515987" C$ h) Z2 Z5 h
    ## SmoothIPW 0.19645006    0.3144263 0.5108763
    7 M5 O1 }* k- N' c6 a! x## SEM       0.07478611    0.2632127 0.3379988
    - T! n% h' [3 C+ P* c## Combined  0.19633867    0.3143794 0.5107181
    / y& x$ ?$ n4 e+ D4 `8 w##
    , W# ^! J3 z8 K9 T## $CI
    # c3 m7 c; ~, M6 b3 z) O* s% D. u## $CI$NRI.event
    / Z; _% j2 s# o: T##                     KM         IPW   SmoothIPW        SEM   Combined: f3 Y* B% r- c
    ## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723
    ' s; ?9 k) u4 p' y% F9 H  W3 k* Z## upperbound  0.44806768  0.47033936  0.44014214  0.2658309  0.4400496
    - c+ @9 U5 \0 Z" U## ; J. W5 K; x, V5 m/ F# ]
    ## $CI$NRI.nonevent! b- \, S/ ?3 h
    ##                   KM       IPW SmoothIPW        SEM  Combined
    * U. R5 @! ~! z' Q+ }- P% u# z. R## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.1286426
    ! {6 h  B& q4 ^# l! v0 ?0 F( w## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.6964549
    2 q  r( E* Z6 W& @) M, R1 Y## 7 F6 m; M: W7 w  |9 i! W+ @
    ## $CI$NRI
    8 o: c, X& ~4 d# O##                     KM         IPW   SmoothIPW         SEM    Combined4 S& {* \  a) J9 {' R0 ]
    ## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409
    ) {- G; O, |/ B/ _## upperbound  0.89306122  0.92464359  0.87970125  0.64253510  0.87953153
    % S4 N- A+ d7 `0 s! m4 E' f## 3 f- m( ^! ^: o* r4 O
    ## ( C" u# [; h( R/ ]3 x
    ## $bootMethod1 s! K: \  R1 s3 {9 F& c  k% b* ~
    ## [1] "normal": o: h! ~- W+ {
    ## . c# o5 n5 \; g* |/ V
    ## $predict.time
    % c+ O1 k. }3 u- Q9 t  V## [1] 20005 |9 I. T! K- R! d
    ## - l. g! X+ x5 A& V1 J& T/ e( n- c
    ## $alpha
    2 G, |0 }( q6 a& W% z  T0 Z" a## [1] 0.051 J8 C2 M5 \7 \' j8 t  R% {
    ## # K( [" L7 G7 ^) v% p
    ## attr(,"class")
    / C5 X: ^0 Y3 O8 U3 O## [1] "survNRI"/ S: ]) J: P4 {
    $ k2 f: L5 |- r. x$ s
    1$ i4 `6 r4 q: H( V& s6 ?
    OK,这就是NRI的计算,除此之外,随机森林、决策树、lasso回归、SVM等,这些模型,都是可以计算的NRI的,后面会继续介绍。大家如果有问题欢迎在评论区留言。! G% Y8 a5 I4 C" o5 L% V

    ( T. S' }- R5 [! L, j本文首发于公众号:医学和生信笔记) o! N8 q, r2 X# \5 I

    7 R( {$ V- f- E8 i! z“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
    / }' g* s; W  h本文由 mdnice 多平台发布
    + }, U; V: O) E9 O2 T————————————————9 ]) T1 p% I. s% ~8 S0 ], [
    版权声明:本文为CSDN博主「阿越就是我」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。) D7 \0 U7 f3 H, i4 q6 F1 C7 n: j
    原文链接:https://blog.csdn.net/Ayue0616/article/details/126768006
    % {# I" r6 Y  P0 O- ?4 |5 N4 {% v; w3 G' Q

      a7 S) I6 u9 d& ^, \1 K
    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 19:42 , Processed in 0.452495 second(s), 52 queries .

    回顶部