QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 3054|回复: 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 {) X+ e5 J! t净重新分类指数NRI的计算
    - x, n+ s  \* U6 d- g. U9 u“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。* x, ]; y* |" Y$ `. u( }" |
    NRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!- O! |) f% m; c  e* j
    6 T0 t& |" p# J0 Q( J( v
    在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens包,PredictABEL可以计算logistic模型的净重分类指数,survNRI可以计算cox模型的净重分类指数。" a: u! M. I4 @9 A5 i) ^' e

    7 J7 W+ a/ V) t/ c( ^  F6 s9 Ulogistic的NRI3 i) r1 i' j% M1 J" S% O: g
    nricens包( p4 i9 m# X  ~2 L
    PredictABEL包4 c4 J- C& u2 y: W- Y
    生存分析的NRI4 {7 l) Y6 G) F+ y( ?' |/ y- \
    nricens包
    0 n/ y' S; R. @4 m/ csurvNRI包- X1 L! k, @/ b0 h
    logistic的NRI
    3 Z: R. H# v- Z3 r0 ~, q0 anricens包
    ! D" A% a5 ~$ X#install.packages("nricens") # 安装R包
    8 E: J: ~9 i0 [0 M5 b2 I+ B: G7 Dlibrary(nricens)& W) P, t" i6 g+ H* {) Y
    1
    3 \, a  Q/ w8 {( M! p  c$ {9 r## Loading required package: survival9 W! O3 ?1 N. g; z  |
    14 ~% X2 T* h1 f  `
    使用survival包中的pbc数据集用于演示,这是一份关于原发性硬化性胆管炎的数据,其实是一份用于生存分析的数据,是有时间变量的,但是这里我们用于演示logistic回归,只要不使用time这一列就可以了。
    2 E. U6 C6 A* M
    " U: ^. e  p$ |library(survival)
    + Q  y: P/ [8 [) U% {, _
    , E  e3 Q9 E. p$ R$ f# 只使用部分数据
    5 _/ O6 U; |1 t8 y7 \6 B6 m! Tdat = pbc[1:312,]
    ( G9 W5 B# ~. \; d; k. I3 V. Vdat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]
    . `2 ^5 C# _& x4 O+ i6 }- r" ~
    : F4 ]8 G9 I/ E$ z* Astr(dat) # 数据长这样& Q5 }8 ^' l' P, Q) {
    1
    / d4 E' q7 _* }* m+ u+ v7 {" j## 'data.frame': 232 obs. of  20 variables:: x) |# a  P( r
    ##  $ id      : int  1 2 3 4 6 8 9 10 11 12 ...
    % V6 m1 n5 D8 F7 B: V4 {##  $ time    : int  400 4500 1012 1925 2503 2466 2400 51 3762 304 ...
    $ H! C2 q$ T& N0 W##  $ status  : int  2 0 2 2 2 2 2 2 2 2 .../ j% V; V0 H3 u: n
    ##  $ trt     : int  1 1 1 1 2 2 1 2 2 2 ...; S- X) O/ ?0 j/ w8 q
    ##  $ age     : num  58.8 56.4 70.1 54.7 66.3 ...
    4 z: Q/ Y7 [& c7 S; h##  $ sex     : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
    ' q) ~# C8 r8 ]) H4 J##  $ ascites : int  1 0 0 0 0 0 0 1 0 0 ...! r3 ^* j! Q* ^/ k3 f
    ##  $ hepato  : int  1 1 0 1 1 0 0 0 1 0 ...
      p3 G. ]2 A- ?  Y" j" Y9 q##  $ spiders : int  1 1 0 1 0 0 1 1 1 1 ...
    % d8 S5 a3 F/ S##  $ edema   : num  1 0 0.5 0.5 0 0 0 1 0 0 ...! z# I1 w% d5 H! D  D* o9 a: _. r
    ##  $ bili    : num  14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...6 I0 H% W3 K' s- a  q% c9 c
    ##  $ chol    : int  261 302 176 244 248 280 562 200 259 236 ...7 K! s5 k2 k  c. X0 @0 v% a
    ##  $ albumin : num  2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...6 i+ q  D# I. Q% k' _
    ##  $ copper  : int  156 54 210 64 50 52 79 140 46 94 ...
    2 z0 W% e# E8 M( o6 T+ @##  $ alk.phos: num  1718 7395 516 6122 944 ..., F& W1 c6 P+ a4 I
    ##  $ ast     : num  137.9 113.5 96.1 60.6 93 ...
    1 _$ t; A) W% q$ h( Z, `##  $ trig    : int  172 88 55 92 63 189 88 143 79 95 ...
    & V# @3 Z# m/ \. |5 J9 f##  $ platelet: int  190 221 151 183 NA 373 251 302 258 71 ...
    : u, f$ N, z9 U0 i& J' {##  $ protime : num  12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...
    3 T3 m3 K" x8 C4 b6 p##  $ stage   : int  4 3 4 4 3 3 2 4 4 4 ...4 z  k5 {* t/ D( V3 c* i
    / H" y0 o2 {8 B6 n+ h2 h
    1; a8 p& k8 E6 Z& }% {
    dim(dat) # 232 20
    ! [3 R& Z6 r) u8 _) i& U3 K1  E0 V6 ^3 j. o' N5 h! D
    ## [1] 232  20) f/ j- i$ D. k% j
    1
    / s* Z0 d& t" K5 x# [然后就是准备计算NRI所需要的各个参数。! Y8 o6 C) Y: Q* l" D9 B2 E
    4 Z% |0 [5 p. B6 l4 m7 U% ^/ p) ]! E
    # 定义结局事件,0是存活,1是死亡
    , T: O+ J* j. [1 S4 [" fevent = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)7 d/ A& i" @2 _( m0 b  I

    1 t" i- X. f) X% T' ]# 两个只由预测变量组成的矩阵
    , k3 C- ^7 t, |% D% Q. X# i- {z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
    % I. `7 i$ v1 W7 W/ k" [4 mz.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))4 B1 }$ @; B% i" }/ |

    % `1 q0 @" E8 t( x, F# o2 k, A# }# 建立2个模型9 m, A2 _4 P" W" i/ D& c/ P" y
    mstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)* X. g6 F5 X8 W2 x/ {" l" y
    mnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)% _% H- v# i7 Z% F

    7 }+ |3 R6 C) A2 I8 g# 取出模型预测概率
    1 t& Y, t8 Q) D" v, R' Fp.std = mstd$fitted.values0 n5 F/ b: D9 u3 t7 h$ y% D' _
    p.new = mnew$fitted.values
    $ r6 b0 B- `0 O: P4 i. {- ]7 `
    - q# H4 r7 d7 c* a1
    ; x* S; s/ ]. q然后就是计算NRI,对于二分类变量,使用nribin()函数,这个函数提供了3种参数使用组合,任选一种都可以计算出来(结果一样),以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。
    & l# U$ h7 Q7 ?2 Q% o$ l) O' X  B# |# d8 p7 S- j  u% o
    # 这3种方法算出来都是一样的结果: M1 b' \1 e) @) ?9 {8 `1 ?

    % Q- t0 ]- o! \) g( g/ Z# 两个模型
    : Z* @# I5 y$ onribin(mdl.std = mstd, mdl.new = mnew, 2 g6 a/ t0 y4 z
           cut = c(0.3,0.7),
    ! C6 P) N" R" L6 N$ n! m       niter = 500,
    . ]+ i' L+ G- i       updown = 'category')6 T" G( b' B# S5 B. J4 b0 A7 o
    % j" Q  l2 i; ]' K# n4 r. H' [; x
    # 结果变量 + 两个只有预测变量的矩阵& p& e, t% l6 P$ d7 I/ g
    nribin(event = event, z.std = z.std, z.new = z.new, & D7 n1 e' z6 d* S; v
           cut = c(0.3,0.7),
    & Y( x" x5 B8 A6 d4 [5 J2 V: z5 c       niter = 500, - s! k0 `; w' Y
           updown = 'category')
    # ]+ ^9 F8 q$ j$ I" ]
    - O9 A0 m$ W+ n, I  i( {' H9 z, ]## 结果变量 + 两个模型得到的预测概率
    ' L# M5 c! t+ T* E" T3 U8 pnribin(event = event, p.std = p.std, p.new = p.new,
    : j7 ]$ Z0 j: _8 u- N! ~2 p9 _) o; m       cut = c(0.3,0.7),
    : ]8 ?+ B; p2 a1 ~: }# H/ K+ B       niter = 500, 4 O2 E/ i6 s9 r- _
           updown = 'category')7 o/ C9 h: p& H- W. @  [
    5 F1 E8 j  ~+ @  z7 U5 b
    1# k! _) d( H: H
    其中,cut是判断风险高低的阈值,我们使用了0.3,0.7,代表0-0.3是低风险,0.3-0.7是中风险,0.7-1是高风险,这个阈值是自己设置的,大家根据经验或者文献设置即可。3 [; O* M7 v2 n0 p

    $ l: H  j# f: ?; zniter是使用bootstrap法进行重抽样的次数,默认是1000,大家可以自己设置。
    $ R8 c# V; f+ }2 @7 s
    3 w5 `+ B( _4 Y9 T; ^updown参数,当设置为category时,表示低、中、高风险这种方式;当设置为diff时,此时cut的取值只能设置1个,比如设置0.2,即表示当新模型预测的风险和旧模型相差20%时,认为是重新分类。
    4 r' t; s/ _; t; w5 b  ?+ w% T) H# r4 _; X
    上面的代码运行后结果是这样的:
    & E* o) f6 M, J. Q/ J) j9 B
    , c6 S. g3 a* V% MUP and DOWN calculation:
    1 P; I8 ]0 T8 a/ q4 R) G  #of total, case, and control subjects at t0:  232 88 1443 U/ c4 m1 o1 w) @2 b. k
    + p/ I2 E& Y3 m8 G& F  P4 |
      Reclassification Table for all subjects:- Y1 Z4 Q# k7 ?0 r/ u
            New
    $ e( Q  l2 |& {- f4 @% r' bStandard < 0.3 < 0.7 >= 0.72 Z; z5 Y: R1 K$ x) I& M. E
      < 0.3    135     4      0) }# M9 S  G; m9 u
      < 0.7      1    31      4
    " A* Q; U6 R. Z( y! ]4 W- i  >= 0.7     0     2     552 y& R2 b3 T, p& P$ b* `9 z
    4 {5 X' w$ }& g& D! [
      Reclassification Table for case:$ u" Z1 s, q: Z$ U6 q3 c0 s7 y
            New
    ! ^# K: }. f% A, u2 J" `Standard < 0.3 < 0.7 >= 0.7" x  K- e$ ~, P7 }! Q3 I
      < 0.3     14     0      0
    ) n# Z4 C; c) b  < 0.7      0    18      3; b- A" S4 m; \: f8 f9 V/ ]. A
      >= 0.7     0     1     523 e: I# {8 k1 l$ E: _0 ^+ v; I1 g
    0 H1 T/ \0 k3 B$ c2 f* s
      Reclassification Table for control:7 c, s% N) H* t1 f/ a
            New( M9 {2 _  v! g5 a1 }5 U. y" s
    Standard < 0.3 < 0.7 >= 0.7( n  N7 f: ~# y; U
      < 0.3    121     4      0* M# c5 X+ X2 x( Z
      < 0.7      1    13      1# |( t6 Y( Z) B7 G
      >= 0.7     0     1      3" W& y7 O/ ^, m+ ^

    4 c/ @7 t8 r, A  c5 X6 ?NRI estimation:; D$ K* S, F; Q# V  P
    Point estimates:/ p" H* D5 @& |# d
                      Estimate
      s8 M5 \3 @& {2 j  @NRI            0.0018939394 s% J/ q" u1 p& L6 S
    NRI+           0.022727273
    0 k! a9 p: h5 |7 P3 A' [: hNRI-          -0.0208333331 O0 L; |9 |* \" k. I$ X9 x
    Pr(Up|Case)    0.034090909; j5 d0 E" I4 b* k% g/ k
    Pr(Down|Case)  0.011363636
    * p5 h8 k8 x" G3 b3 tPr(Down|Ctrl)  0.0138888895 }* |" k- a" r9 U( M
    Pr(Up|Ctrl)    0.0347222226 c8 z* y! I# N1 }. `

    , l5 |: r- @. a' O% A! |) PNow in bootstrap..2 L) L% p9 D& K4 Y+ P
    ; F, V: }) `% j; K8 U
    Point & Interval estimates:
    . V  t7 e9 D2 [  l! f                  Estimate   Std.Error        Lower       Upper
    ( p- R" L$ k8 E* hNRI            0.001893939 0.027816095 -0.053995513 0.055354449
    ; I- {7 k* F! ?7 e8 w/ DNRI+           0.022727273 0.021564394 -0.019801980 0.065789474
    ; H8 n& ~/ @% K9 A" U; k: Z8 }NRI-          -0.020833333 0.017312438 -0.058823529 0.007518797
    " \8 z6 F  }5 NPr(Up|Case)    0.034090909 0.019007629  0.000000000 0.072164948
    . B8 h5 m" I  k, \9 o0 K' @, }Pr(Down|Case)  0.011363636 0.010924271  0.000000000 0.039603960
    4 T( Y6 s* D6 q7 nPr(Down|Ctrl)  0.013888889 0.009334685  0.000000000 0.035211268, ]4 E' q$ n5 |& G
    Pr(Up|Ctrl)    0.034722222 0.014716046  0.006993007 0.0661764719 l$ ~: m4 y* |. ^8 X

      A1 O6 P- U+ D0 Z% P% L1
    " Z9 E2 m% o  n) `8 `首先是3个混淆矩阵,第一个是全体的,第2个是case(结局为1)组的,第3个是control(结局为2)组的,有了这3个矩阵,我们可以自己计算净重分类指数。; w% f) s3 w2 U6 s  `- @

    8 W' j4 o7 q+ K7 h0 t看case组:
    , O+ k3 n9 l3 q
    ( T. i- ^6 C4 l( E- _5 U" q& a净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273- [' X% {: x; P
    + l4 N7 [* A4 Z; {# J* ~" w
    再看control组:8 D9 Q- [+ h' P9 D
    9 L4 A7 I" q% u+ S1 s1 [
    净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333
    % j; D' o; l4 T" v- `$ ^% B# s1 |, R& B, ?
    相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.0003156572 C+ k8 ^" Q1 O& K3 T# W$ Y1 W% y

    ! e7 p0 K. ?7 }, b2 C: n7 b再往下是不做bootstrap时得到的估计值,其中NRI就是绝对净重分类指数,NRI+是case组的净重分类指数,NRI-是control组的净重分类指数(和我们计算的一样哦),最后是做了500次bootstrap后得到的估计值,并且有标准误和可信区间。
    % q: S  }$ Q  i$ R1 V+ p
    & a% W7 e3 D: f2 G  o* F$ _0 j最后还会得到一张图:
    - X. R6 y- `3 K9 P& V" \) w# Z. f2 P- X( U- v
    这张图中的虚线对应的坐标,就是我们在cut中设置的阈值,这张图对应的是上面结果中的第一个混淆矩阵,反应的是总体的情况,case是结果为1的组,也就是发生结局的组,control是结果为0的组,也就是未发生结局的组。! x. e- m6 w' t( f6 F

    " y. C5 ]2 F# z0 p& GP值没有直接给出,但是可以自己计算。+ @: j9 b3 |' U' X5 c' ?

    0 ^; S& l* T+ V1 j3 w# 计算P值$ `: P  R' Z$ a0 Z( @, T) [
    z <- abs(0.001893939/0.027816095)
    ! d; ^5 S0 O$ Jp <- (1 - pnorm(z))*24 b1 l4 B' Y* U! V2 Y
    p
    9 m8 _, L1 x) ^1; K: `5 z( \) z: v% m& \7 I
    ## [1] 0.9457157
    & c. J1 v1 g* A3 ?4 C7 C1
    2 t2 d* \& k( X" _4 ?# F& CPredictABEL包
    ' I% b  N( Z- }" L3 q0 |#install.packages("PredictABEL") #安装R包
    ' I8 W" j  l1 D3 D! Ulibrary(PredictABEL)  1 {; p0 U& x* N9 e0 ^
    / s' m2 ^7 }! e: Y/ R9 D6 t
    # 取出模型预测概率,这个包只能用预测概率计算
    " ]% K. x& x3 S& Hp.std = mstd$fitted.values' o# U# T* w2 o1 _+ p0 i/ a. d$ i
    p.new = mnew$fitted.values * ~6 S, D4 n9 k8 h/ s
    1
    5 W5 ]  h* J3 b' l然后就是计算NRI:
    . k( h+ z+ m  p/ v, v- Q' v. x* s- r% b" F
    dat$event <- event
    4 w1 i: D5 K& ^+ o4 {
    3 Z( T' A- P' g, w1 t6 A4 Treclassification(data = dat,
    4 n% F4 r/ r- K* k) A6 _# `                 cOutcome = 21, # 结果变量在哪一列( o& Q  u9 j! s
                     predrisk1 = p.std,
    + P4 m* Y& r8 n                 predrisk2 = p.new,: ]8 V4 a; i. \( \0 G
                     cutoff = c(0,0.3,0.7,1)0 Q7 p' \( U7 e2 \1 A
                     )* o! Q) {' @# F# Q0 x3 Z; H9 Q
    18 K* p* Q$ S7 A8 ]: `" \
    ##  _________________________________________0 |; E" x6 `) [4 S7 I& {/ X
    ##  ; v" L; ]! ]! e$ y% F( w
    ##      Reclassification table   
    . Q1 M0 ?1 V2 F1 |2 f- r. ?' N. p##  _________________________________________" A2 J! E7 v, c. i- H% l% U$ i
    ## 6 [& p5 G& D+ s8 ~5 c1 }
    ##  Outcome: absent
    8 _0 P3 b3 w! ?! Q& c! f##   
    & ~  X- [; O- u& b3 A1 `##              Updated Model
    1 s1 W3 m5 N/ ?  V' D; Z## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
    : L% V  ], k; P1 W4 o##     [0,0.3)       121         4       0               3
    + S  H) r& S" }' t0 ^6 R. h##     [0.3,0.7)       1        13       1              13
    / |1 S  q; J" k! F+ b4 E; t1 R/ ]##     [0.7,1]         0         1       3              25
    9 S1 b/ t2 f6 {/ h4 a  S## 1 h" |( U9 ^! j: k2 q0 p
    ##  2 C5 ]' {- x: p( M# Z$ R9 S9 |
    ##  Outcome: present " M/ a* ^1 G6 E3 C6 f( Y
    ##   
    9 ~; {" l2 N$ ?6 D! b& t6 i0 a% ^##              Updated Model' ?* ^* i0 e8 X8 ^7 n: H/ K
    ## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
    * n2 M- G" W  D* C1 w! k/ c##     [0,0.3)        14         0       0               0
    + T  r- l$ p& v: \1 s' q8 s##     [0.3,0.7)       0        18       3              14
    # r# v/ @7 B. b3 N& t##     [0.7,1]         0         1      52               2
    7 d8 m/ K4 ^* X( c; s* P3 Z2 x## 4 r& ]5 s% F# b6 _- P" L8 R
    ##  
    & l7 j0 }" O4 [7 T( e- L; m##  Combined Data
    4 m4 Q3 L) e6 u+ O7 y, K; r( n##   
    3 _: Z: n1 V/ L% W5 @) N##              Updated Model& m6 E" i! y* l( J- {
    ## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified1 n) R0 \& G* f/ {# h. m1 |- F0 h
    ##     [0,0.3)       135         4       0               3
    # F0 ?2 O0 h2 o& J; A  O3 @##     [0.3,0.7)       1        31       4              148 D. l& Z% O+ J, A2 O& L) X
    ##     [0.7,1]         0         2      55               4
    ' U" F7 B; D" p' X8 t- o##  _________________________________________
    4 @" g& E1 N0 v0 W## / E9 p% g2 l1 U8 @
    ##  NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806
    ) r; r. e& e4 [" b7 a##  NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048 3 H* Q. P4 g! M8 O' T$ ?
    ##  IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.28396: ?! ~% y7 C+ h2 J
    + n8 X3 \1 t5 L1 ~( c% T  b
    1
    5 A, S, V9 A8 _3 X' @: l) p( Y" w结果得到的是相加净重分类指数,还给出了IDI和P值。两个包算是各有优劣吧,大家可以自由选择。! z% e6 G9 E  s6 }
    $ P( {2 f) E( P9 l6 Q& G' `
    生存分析的NRI- l/ ~" w; I6 p9 F
    还是使用survival包中的pbc数据集用于演示,这次要构建cox回归模型,因此我们要使用time这一列了。$ e0 s4 _! u# e& i& f4 \
    : \- Y; F: {, I0 K  k. h# ^3 m
    nricens包
      O: J7 a( u: U% k4 vlibrary(nricens)
    % o* Z  R  j0 k! P1 {' Z6 ~library(survival)
    + J+ J( f" Y, B
    $ Q+ t4 e9 X! k5 Edat <- pbc[1:312,]
    & D0 g- X  d2 e% {' I0 {dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡
    3 C1 N, B2 _$ J12 @6 g9 T- ]! n
    然后准备所需参数:- C# d1 V4 F* H2 v& L1 O; e8 L& z

    + q# C0 J: Y2 D# 两个只由预测变量组成的矩阵
    7 j" e" R7 a9 V) {7 `z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))4 P5 X9 I7 \! Z: w. r
    z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime))), h. ~6 T' v8 g3 y% V( q

    & f" _8 N: ], m9 z. j# 建立2个cox模型
    / @6 O2 A0 {$ j' x3 M; F" H8 o( gmstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE)
    2 _( k9 I- e: x8 w/ Imnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE), H3 ]" s3 q: p

    8 ?* Y! m$ z1 E$ A; t0 ?4 v1 f# 计算在2000天的模型预测概率,这一步不要也行,看你使用哪些参数
    . w2 j9 n3 j; F/ }9 r( f6 R' Xp.std <- get.risk.coxph(mstd, t0=2000)" i8 W' j7 X% G- @7 x  ~9 ~
    p.new <- get.risk.coxph(mnew, t0=2000); I  P2 u6 B/ T% U! Y; \5 L0 ~
    1
    0 [1 {3 ]$ H# _4 U: j计算NRI:
    * u4 H# H' ~% b( b" _& {5 i  i; M6 }# g
    nricens(mdl.std= mstd, mdl.new = mnew,
    # Q: V8 u4 I! p+ X; W" t5 i+ A        t0 = 2000, : F; ?* W3 h8 n2 m' K- I' \
            cut = c(0.3, 0.7),
    ) q5 T( x1 Z9 C; E: P* X% m* \        niter = 1000, # U2 o) N. G3 T$ @) c
            updown = 'category')
    ' A- S' S/ U' z+ t2 {; }" k) R' Z
    # Q; T, E. G2 M) UUP and DOWN calculation:
    - F) d* V) }) G  #of total, case, and control subjects at t0:  312 88 144  y2 Q5 T. n9 @' o* D! j& i7 S( j

    9 I3 w% \% G$ @, N6 w  Reclassification Table for all subjects:  @) b' {0 O. G5 c& E0 a
            New
    ! @) y( `4 w# J8 ]6 UStandard < 0.3 < 0.7 >= 0.7# c1 Y9 g& n3 Y' I+ a
      < 0.3    202     7      0, b7 C" j2 _5 j+ g, |% ~
      < 0.7     13    53      6. O  n9 `. q  _% _4 ?6 G8 p4 j. z1 C3 d
      >= 0.7     0     0     312 [" x9 p/ F! W/ C8 ?' r0 B% P
    6 A2 s: g6 z  X4 X6 p
      Reclassification Table for case:- g/ l+ H1 C6 D8 ^! a, x2 e, X
            New
    4 G9 \  E1 K' {Standard < 0.3 < 0.7 >= 0.7
    7 n, n6 a+ h* Q: Z# P& d2 s5 D$ x  < 0.3     19     3      06 H/ q" ^( L6 Q3 f( o
      < 0.7      3    32      4; Z0 j; d# l3 e5 _; Z: Q  Q7 {
      >= 0.7     0     0     27
    ; N5 C# ?$ d; c% r+ E% x& A! ~' N- d9 V6 s! u3 t7 o$ v/ ^
      Reclassification Table for control:
    $ W) l  p: h  I. \  s( u6 c        New
    7 |) T- Y1 h# |+ B# x! i: `Standard < 0.3 < 0.7 >= 0.7( b7 }, P% K" y. p! s- G  b6 D
      < 0.3    126     3      0
    5 S0 z: k! v& f/ r  v  < 0.7      5     7      2) i1 ?6 I/ g$ e2 a: s
      >= 0.7     0     0      1- E5 ]+ ^$ r- ~

    & U0 L6 Z9 }8 f& H! i6 c2 gNRI estimation by KM estimator:# L' J0 Z% B1 R4 l' @+ w" [
    2 s4 @9 K) b% H0 Y! b7 K: R
    Point estimates:4 U, r! `  D4 w9 E; S3 r/ |
                    Estimate7 e8 m# I/ _. Q% [0 g7 B
    NRI           0.05377635) _3 _: D; n- M7 Y0 f* G  [2 \
    NRI+          0.03748660, I) ?* z+ @6 w  ]
    NRI-          0.01628974
    0 h* m+ x2 i6 s) e2 aPr(Up|Case)   0.07708938
    4 N' s# c) z: ]* jPr(Down|Case) 0.03960278
    2 s3 t, N$ J  Q" h) sPr(Down|Ctrl) 0.04256352  Q7 O# X' N. O  C
    Pr(Up|Ctrl)   0.02627378
    6 k% X8 N' F& |  a: B; ?# q+ W; C5 X$ z6 d9 G0 T" `
    Now in bootstrap..  T& h$ G5 k1 S- \

    " |! i# n  G3 S( v" OPoint & Interval estimates:
    2 J/ Y$ @7 m" ~8 b; {# G7 C                Estimate        Lower      Upper; @0 w2 Z9 J0 z, @( Y
    NRI           0.05377635 -0.082230381 0.16058172
    : n  o+ V1 J" U4 V- CNRI+          0.03748660 -0.084245197 0.13231776
    0 _( k* t, L: u6 KNRI-          0.01628974 -0.030861213 0.06753616
      U1 e8 l( n9 r* V% j+ }Pr(Up|Case)   0.07708938  0.000000000 0.191022916 F% D" q* L( v' [7 i
    Pr(Down|Case) 0.03960278  0.000000000 0.15236016) G8 [& d  B6 R; V
    Pr(Down|Ctrl) 0.04256352  0.004671535 0.09863170
    + `, k* p6 M  x% K9 ZPr(Up|Ctrl)   0.02627378  0.006400463 0.05998424
    9 v) |- F$ Z9 A' ~! r' y) X& ~
    ) P# h) _  }& z3 ?2 @4 Z9 f' x1
    4 g9 X) X9 a4 G# K4 t* Y  F( G8 |4 ^- T% C; Q1 i5 z& E
    Snipaste_2022-05-20_21-49-38
    & g: O; u6 l$ t# R" X结果的解读和logistic的一模一样。6 k6 X6 A* s' R

    4 I% U- B/ M3 e& |2 q9 {0 gsurvNRI包
    / a8 B" m: }) z  u0 l1 e8 [: `# 安装R包
    $ ?, s! Q+ N/ [( Q8 z  A! L3 edevtools::install_github("mdbrown/survNRI")
    ) [) k9 X  G# u0 u0 J' g1 w- [/ g3 V( N1& Z5 Y  K* v2 p- ?3 z- u' K
    加载R包并使用,还是用上面的pbc数据集。
    ! e+ g4 M3 |3 X2 [' {
    2 e! K4 V9 Z+ n) C9 q6 O  Xlibrary(survNRI)
    % F- Y* y1 \) d1/ G0 M+ [+ K6 T; @
    ## Loading required package: MASS) E# K0 ]2 `: l0 F
    1
    1 e7 Z5 X# Q6 _3 o5 k7 blibrary(survival)
    ; J5 Y9 {/ ?' j( |, L
      K6 J8 Y1 M7 m# L" v. _. f  Z: h# 使用部分数据
    ! v2 f3 m7 L/ ^/ fdat <- pbc[1:312,]. a7 R2 m3 x$ s2 q
    dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡" g: i5 S: ?* o; R# _' ]* [; F) ^9 t3 w

    * m. u+ q2 Q3 y; Sres <- survNRI(time  = "time", event = "status", 0 j: [; H% }# [8 `/ U+ S$ K7 _
            model1 = c("age", "bili", "albumin"), # 模型1的自变量
    ; f' A0 O, C" j; p' D0 P# V* y        model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量
    ( A6 P4 Z# ~+ q+ A        data = dat,
      ~  N  {( \- r) e! N2 c$ i        predict.time = 2000, # 预测的时间点3 z9 ^; z$ m  v/ X; w0 A; r) s8 C' z
            method = "all",   r* A9 G' [8 t: ~$ s. R
            bootMethod = "normal",  " ?2 z1 Y/ n  D( m" n6 w
            bootstraps = 500,
    0 C" T1 ^- r$ J. K2 C5 [        alpha = .05)
    , ]: q' k1 b& d( R4 z1 G: c
    ) w! K4 g6 j; g. x  y) ?5 B1
    7 g& X; I, s1 C7 `6 Q: |查看结果,$estimates给出了不同组的NRI以及总的NRI,包括了使用不同方法(KM/IPW/SmoothIPW/SEM/Combined)得到的结果;$CI给出了可信区间。6 ^5 H5 C0 W) H, ?
    8 h) n  y* G5 D7 u" `
    res8 |- U% T2 z; h5 o. o, p+ t
    1
    5 j' I' U0 `2 g## $estimates
    , o8 O( |& K  j2 f, k4 k  p* |" P##            NRI.event NRI.nonevent       NRI) k& @1 Z: Y2 [1 \, ^
    ## KM        0.20445422    0.3187408 0.5231951. B& H) R6 G# D, G
    ## IPW       0.22424434    0.3273544 0.5515987
    , y: J# S: J$ z: T6 E: Q## SmoothIPW 0.19645006    0.3144263 0.51087635 L: o0 i1 d7 z( ^0 f4 z/ ]
    ## SEM       0.07478611    0.2632127 0.33799881 P) [' l( j6 Y) V- z
    ## Combined  0.19633867    0.3143794 0.5107181
    % P9 o0 m; I7 n, e% [/ ~## 7 x9 T, ?& c* t5 w
    ## $CI
    1 s# e5 M# C  u, l* ^; l2 p## $CI$NRI.event' R; m- q  ?' m+ Z, g
    ##                     KM         IPW   SmoothIPW        SEM   Combined
    , ]  U. w+ b: `% o) K1 Q& e: F## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723
    0 B+ t# l( n! ?( R  j## upperbound  0.44806768  0.47033936  0.44014214  0.2658309  0.4400496
    9 [$ A0 o" V: R  a7 k$ ]7 J## ! H0 b, h( c/ B- L  b% E
    ## $CI$NRI.nonevent& N5 p: L" K, _  Y( Q
    ##                   KM       IPW SmoothIPW        SEM  Combined( w! \* I* z  ~! T
    ## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.1286426- p. y1 e# K* Q! f  p' L0 `0 m
    ## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.6964549, `7 F5 I! B) h# Q0 `
    ##
    % @6 G) g& h7 M* C  H## $CI$NRI+ l: l/ N/ t3 k) N' u% ?8 R! o
    ##                     KM         IPW   SmoothIPW         SEM    Combined
    ! _- S* x, b( S4 X; N; {## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409
      n4 U  V1 Y0 [" |. M# s## upperbound  0.89306122  0.92464359  0.87970125  0.64253510  0.879531538 \& U) R  J  ]' v% t
    ##
    * L0 u$ j4 L' N& }  n1 T6 A##
    $ }- ^  \7 c. c: H# h# L3 D& N6 w## $bootMethod1 [% t& n6 m4 k
    ## [1] "normal"
    + w# [& U) K; P& }& ~5 E5 A## & [9 x! a5 M& I. ?! i8 Q
    ## $predict.time) d4 Y- A# i9 U8 [) y2 ~! }/ y
    ## [1] 2000/ y- F1 B1 e0 k3 I4 T
    ##
    . u$ v; G9 y7 ?## $alpha5 o- N1 F1 B# c! n6 z( j5 k4 y) k
    ## [1] 0.05
    ' I; A, U+ l( X, b/ ?8 q! O## ! N0 x2 x. ?$ y' J" a! _2 |/ {
    ## attr(,"class")
    / |% N. ^0 m$ a; l## [1] "survNRI"
    6 B' ~* B9 J% V! e2 G" E! ^* I5 s! c& _! E; Y9 N( B
    12 v5 u* O! l. W) @1 `/ l
    OK,这就是NRI的计算,除此之外,随机森林、决策树、lasso回归、SVM等,这些模型,都是可以计算的NRI的,后面会继续介绍。大家如果有问题欢迎在评论区留言。4 H& d3 h( {3 F/ _9 f6 N

    2 w# F3 B; O5 o2 O本文首发于公众号:医学和生信笔记
    : \2 k- G4 L$ V) D. N/ W; h; z+ a: z# l' x; L
    “ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
    9 b. _+ l% U' z( \3 @+ ]6 U本文由 mdnice 多平台发布+ z2 f& j$ C6 n- m" E
    ————————————————
    & p/ h# _' Z5 l' O( s$ i版权声明:本文为CSDN博主「阿越就是我」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。$ B! d) \! ~4 U
    原文链接:https://blog.csdn.net/Ayue0616/article/details/126768006' j2 U2 b" H; \6 e
    9 s! Y- g3 j) E9 }: G" W1 @6 P. b
    / p3 I- F" e0 Y4 R1 r! {% ]
    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 10:55 , Processed in 0.445602 second(s), 51 queries .

    回顶部