- 在线时间
- 1630 小时
- 最后登录
- 2024-1-29
- 注册时间
- 2017-5-16
- 听众数
- 81
- 收听数
- 1
- 能力
- 120 分
- 体力
- 541607 点
- 威望
- 12 点
- 阅读权限
- 255
- 积分
- 167860
- 相册
- 1
- 日志
- 0
- 记录
- 0
- 帖子
- 5324
- 主题
- 5250
- 精华
- 18
- 分享
- 0
- 好友
- 163
TA的每日心情 | 开心 2021-8-11 17:59 |
---|
签到天数: 17 天 [LV.4]偶尔看看III 网络挑战赛参赛者 网络挑战赛参赛者 - 自我介绍
- 本人女,毕业于内蒙古科技大学,担任文职专业,毕业专业英语。
群组: 2018美赛大象算法课程 群组: 2018美赛护航培训课程 群组: 2019年 数学中国站长建 群组: 2019年数据分析师课程 群组: 2018年大象老师国赛优 |
# S% D; r6 z3 [7 p5 [7 {
净重新分类指数NRI的计算
; } F9 v0 F" U. s“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
* }: f$ ?0 Q1 I: JNRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!
/ U' N# ` k6 d# w6 i* l; o/ f/ D0 l4 ^$ _0 x. T
在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens包,PredictABEL可以计算logistic模型的净重分类指数,survNRI可以计算cox模型的净重分类指数。$ A5 ?1 U& k0 w% y* i5 i; _' s
* M' `" p) C4 L; B$ V+ w
logistic的NRI
7 \$ q9 S) F, n& r( H) p% [nricens包0 d0 T# X- \4 O% E, u5 P5 e
PredictABEL包) n1 k( D, o7 E9 l
生存分析的NRI3 J E2 ?1 Z4 k9 F- d1 k
nricens包
" e: s8 `* b" U1 q4 XsurvNRI包
' \4 I$ \1 K) ulogistic的NRI/ ~ f8 M! W+ P6 {: a d
nricens包
2 r. N( `9 m! V#install.packages("nricens") # 安装R包' a" P% L$ [4 p( M8 b# c
library(nricens)7 m. w, U* u7 b- ?" D& u) l
1- b- x/ T% o7 s; f) o
## Loading required package: survival
2 _0 l% t ?3 q6 ?1! q; F/ S2 h- K( l1 O/ z! F
使用survival包中的pbc数据集用于演示,这是一份关于原发性硬化性胆管炎的数据,其实是一份用于生存分析的数据,是有时间变量的,但是这里我们用于演示logistic回归,只要不使用time这一列就可以了。8 i: l) g0 E5 z
: @) m' d) h) z% F" L2 ?2 S9 [library(survival); m. _: k- O4 f+ `, ?" R2 A
1 y# O; o% s4 w/ F( g
# 只使用部分数据
; j0 m! y. l4 m4 N+ Mdat = pbc[1:312,] 6 L3 X7 [+ `6 e
dat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]& d7 o$ V' k! |
D& P+ x% [. ]. a/ {* u: m! _; c4 v
str(dat) # 数据长这样
( ?: B; }! S$ `: r) g C1: t) _! s' {0 [0 e q
## 'data.frame': 232 obs. of 20 variables:$ K" g ?( S, o5 J& c0 |
## $ id : int 1 2 3 4 6 8 9 10 11 12 ...$ [2 Y4 b7 V) W: e0 d) s7 l/ n
## $ time : int 400 4500 1012 1925 2503 2466 2400 51 3762 304 ...
9 u! ^" |6 e$ s## $ status : int 2 0 2 2 2 2 2 2 2 2 ...
! l: q, |* n5 Q5 N; _## $ trt : int 1 1 1 1 2 2 1 2 2 2 ...
7 d# n& e3 y, l, O$ m( B## $ age : num 58.8 56.4 70.1 54.7 66.3 ...
$ ^) q- ~2 ~8 l3 H1 M- R4 C6 o' O## $ sex : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
4 S* [- e& U, t& B/ x: |## $ ascites : int 1 0 0 0 0 0 0 1 0 0 ...2 W% ~& f! [1 A1 J% t8 ~
## $ hepato : int 1 1 0 1 1 0 0 0 1 0 ...
1 |$ s" o/ K6 B4 q+ x0 s## $ spiders : int 1 1 0 1 0 0 1 1 1 1 ..., ~( \0 w1 A) g* r) S- [. X
## $ edema : num 1 0 0.5 0.5 0 0 0 1 0 0 ...
9 i1 q! V4 z9 W" n/ \8 }## $ bili : num 14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...
6 c9 U6 ?* F, [: s( n9 k## $ chol : int 261 302 176 244 248 280 562 200 259 236 .... z2 y9 E+ z' b3 s* L- {- P
## $ albumin : num 2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...% `/ a0 U3 p* r6 j$ Y ~
## $ copper : int 156 54 210 64 50 52 79 140 46 94 ...
x5 I* Z+ Y' q* z+ ^7 O## $ alk.phos: num 1718 7395 516 6122 944 ...
/ A2 V/ |, \" r; ]0 }, Z& b## $ ast : num 137.9 113.5 96.1 60.6 93 ...; h5 a4 b; g& h( [3 q0 q
## $ trig : int 172 88 55 92 63 189 88 143 79 95 ...% n8 T% e0 u: Z& E# P& {
## $ platelet: int 190 221 151 183 NA 373 251 302 258 71 ...
+ \) ^5 h0 x& R, f3 S## $ protime : num 12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...
4 d! L! ^8 m# {( o( p% L## $ stage : int 4 3 4 4 3 3 2 4 4 4 ...; q$ z) n0 n- S' U5 j1 g
- }+ V0 }' \" Z: N; i/ N2 r6 t11 H1 M3 g4 ]% t( A$ Q
dim(dat) # 232 20
9 k( @1 @( z4 A- Y% _) ~13 s+ {& | Y2 L! U
## [1] 232 20
7 U9 F1 k+ B9 d" R1: Y& u0 |3 O8 ^. R! x
然后就是准备计算NRI所需要的各个参数。7 L) U& V6 w" U
, A0 x5 r! v" c7 K+ N1 s
# 定义结局事件,0是存活,1是死亡
$ E4 t7 T6 Q* x6 y+ Kevent = ifelse(dat$time < 2000 & dat$status == 2, 1, 0) k1 z% }' S; C( B; {. t
; e( ?+ C; C) C! \+ i2 a9 v. ]
# 两个只由预测变量组成的矩阵
" ^( J5 `% I: Z( Uz.std = as.matrix(subset(dat, select = c(age, bili, albumin)))* y2 K( Y [4 G/ ^9 ~( m8 ?
z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))% A% P x# ^0 I: C
8 \+ Z" Q' E' J4 A! E4 M3 `# 建立2个模型
" a# Q8 |' i* K! g9 c" emstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)
' T- f; T2 [7 A d T! ^mnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)
3 a v; f3 Z4 g5 W# A9 X# W" T- o% {! w T# _
# 取出模型预测概率% }- N/ F1 D$ z% W* s
p.std = mstd$fitted.values
0 e+ h6 W4 b( c. F$ Z6 _p.new = mnew$fitted.values
0 h1 b8 a7 r: _9 ?3 l" F4 T3 [$ r
) u6 L+ L" |, i. Z1! u$ w0 c4 r' W* i S
然后就是计算NRI,对于二分类变量,使用nribin()函数,这个函数提供了3种参数使用组合,任选一种都可以计算出来(结果一样),以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。& ~9 z# [! G; _+ p
: G' `$ f# G/ |$ L0 H j1 T
# 这3种方法算出来都是一样的结果
8 g: `9 y" Y8 W, h( o) M# v4 m: X
# 两个模型! P9 U# C3 F& l4 N) K/ i
nribin(mdl.std = mstd, mdl.new = mnew,
' |; |+ k& \; N3 o4 n! I cut = c(0.3,0.7),
) K6 ^8 ~/ H: A niter = 500, 9 b$ T! T" q1 n% B' x( T0 y
updown = 'category')
: V3 t# z }+ @7 V" q: n
$ N2 `: T8 `$ H3 g# }0 @ c9 \3 t# 结果变量 + 两个只有预测变量的矩阵* z H( q. U# n/ z; t0 G
nribin(event = event, z.std = z.std, z.new = z.new,
; d' R2 W+ v2 | p3 w5 O cut = c(0.3,0.7),
* B9 G) y3 i5 H niter = 500,
1 L1 v6 b/ w4 L, Z updown = 'category')9 P$ {. x' `& |: U& w# f/ ^
1 x# w1 U& z, W## 结果变量 + 两个模型得到的预测概率' ^5 F" q4 z2 {: p' B% \
nribin(event = event, p.std = p.std, p.new = p.new, 8 G3 L& u* W% ^2 }% I5 T6 }. }
cut = c(0.3,0.7), 8 |2 ~ F z' F1 X2 X) P
niter = 500,
% t* j8 @2 T. a! m9 }5 ` updown = 'category')
3 G g: Q& {% k8 Y `6 g/ q' _3 }6 u5 U3 [, F. p) S# B
1) r8 H* V0 `3 V: w
其中,cut是判断风险高低的阈值,我们使用了0.3,0.7,代表0-0.3是低风险,0.3-0.7是中风险,0.7-1是高风险,这个阈值是自己设置的,大家根据经验或者文献设置即可。
6 h. p$ J. \0 h0 z- D) O
5 w& }, D- O5 F7 K6 ]niter是使用bootstrap法进行重抽样的次数,默认是1000,大家可以自己设置。# e W# `2 k$ n. f
3 r5 u/ Z9 T0 X g: y
updown参数,当设置为category时,表示低、中、高风险这种方式;当设置为diff时,此时cut的取值只能设置1个,比如设置0.2,即表示当新模型预测的风险和旧模型相差20%时,认为是重新分类。
, q5 B) u. i* u0 V" w" t; ^% q2 K. \
上面的代码运行后结果是这样的:- d- O+ I# u* i* J6 W& N
/ z* S7 c) n- E7 \$ T6 y: w
UP and DOWN calculation:1 w- @9 _# I) @* |
#of total, case, and control subjects at t0: 232 88 144
( I$ I) ?- \, w+ O" M6 D( l- ?# c+ a' F8 Y7 `
Reclassification Table for all subjects:$ Z1 P/ b, G. t
New
3 t5 w0 |9 ~5 I% \: m" e3 [6 BStandard < 0.3 < 0.7 >= 0.7, e& x+ s! |7 |# r
< 0.3 135 4 00 k$ e4 y/ E4 o
< 0.7 1 31 46 S! S7 t5 C/ m
>= 0.7 0 2 557 b: y: j& Q8 ]
. w# X# |5 v& b$ o1 \: ^ Reclassification Table for case:
1 c- {. H3 y6 A9 X, S; T New& R/ v1 Q/ p4 |% g
Standard < 0.3 < 0.7 >= 0.7
( O: V* M( ?9 V! n9 Z& D, {: } < 0.3 14 0 0
2 t; B2 k8 x! O1 e4 M9 Z < 0.7 0 18 3& c* I: ^' ^2 v* ~+ W1 j/ G4 |
>= 0.7 0 1 52
. Q3 C! z9 P% a7 I( Q8 @/ t
6 }& l. v, A2 w' \: X Reclassification Table for control:5 Q( z( g' r; y' g2 \) p
New5 v P: F: X$ k7 o6 i( Z# g
Standard < 0.3 < 0.7 >= 0.7
& E7 G- r C+ @ < 0.3 121 4 0
& d/ F. @( s3 U- R5 D < 0.7 1 13 1: Z+ y4 k$ I+ W1 C0 v8 e8 @
>= 0.7 0 1 3! x9 `% }, ]' s0 h( T
u7 z% V% D) S2 vNRI estimation:- K8 }! L+ F8 x* {
Point estimates:4 K, n! r5 p( V5 ?8 Q3 Z( U
Estimate
/ s/ Z! m. Q0 Y! B+ uNRI 0.001893939# d$ T, N; M n
NRI+ 0.0227272738 ?# H- ]1 C) [/ N# Q
NRI- -0.020833333# h/ p" ~. R n0 q
Pr(Up|Case) 0.034090909
- Z# Y9 J, I3 L( V/ KPr(Down|Case) 0.011363636
) u. A! F0 t% Y) G% q: m, i YPr(Down|Ctrl) 0.0138888897 p* ?2 A3 e }; S5 V
Pr(Up|Ctrl) 0.034722222
5 e% W, |* n8 n. ~5 [( G' q. M
+ c* V, w5 [0 rNow in bootstrap..% X {7 C8 M9 W
# p* h, ]4 |, r/ R( IPoint & Interval estimates:+ a9 E0 v- b4 [7 Y8 F( H7 e
Estimate Std.Error Lower Upper
@% {. [2 T" \) [$ V/ oNRI 0.001893939 0.027816095 -0.053995513 0.055354449
% K# v& R5 T; ]; pNRI+ 0.022727273 0.021564394 -0.019801980 0.065789474
. p+ y( I6 c. C# `+ P5 yNRI- -0.020833333 0.017312438 -0.058823529 0.007518797
9 j: Y0 [( W7 b8 e3 ~Pr(Up|Case) 0.034090909 0.019007629 0.000000000 0.072164948
* z( [1 ]! U/ n# C3 [" O% j/ r! r3 wPr(Down|Case) 0.011363636 0.010924271 0.000000000 0.039603960! }" N. e) W8 k. @% p/ N& h* i
Pr(Down|Ctrl) 0.013888889 0.009334685 0.000000000 0.035211268& g9 i; X& D Q/ p
Pr(Up|Ctrl) 0.034722222 0.014716046 0.006993007 0.066176471' \% J* y( @& v0 i2 e* [4 X
: b/ Z- \9 h# O) }- a
1
. ?% \7 H D7 c首先是3个混淆矩阵,第一个是全体的,第2个是case(结局为1)组的,第3个是control(结局为2)组的,有了这3个矩阵,我们可以自己计算净重分类指数。: G0 B$ }1 S6 w: b
! k0 E1 ^4 N8 o8 @ x; b3 [看case组:' ~3 D$ V( ~# R5 s2 x8 x" ~1 v
$ u8 P' @: Y O1 r- c( z( X
净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273
8 \$ M7 W) L+ M+ A0 Z7 K$ ^! M5 x
8 X8 c& k, B6 h- Y0 J: H再看control组:
# L" w7 x7 {9 p( a# P6 V3 Z* c
1 ~ x* ^& k3 ^* C净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333$ }/ ~0 N1 a' ?* L, ?' d
" V4 k& ]9 E: l$ m相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.000315657
5 Z0 T! a0 |! u: O) y
# k6 |" j& m+ Q3 K# ?& A再往下是不做bootstrap时得到的估计值,其中NRI就是绝对净重分类指数,NRI+是case组的净重分类指数,NRI-是control组的净重分类指数(和我们计算的一样哦),最后是做了500次bootstrap后得到的估计值,并且有标准误和可信区间。9 r* x5 z3 _0 a# T( @- }
! [6 l) _: X! Y0 h" I
最后还会得到一张图:/ c8 T* ]: j1 r5 n1 C
3 P' P4 ~5 A) q) t8 E: e% g1 }这张图中的虚线对应的坐标,就是我们在cut中设置的阈值,这张图对应的是上面结果中的第一个混淆矩阵,反应的是总体的情况,case是结果为1的组,也就是发生结局的组,control是结果为0的组,也就是未发生结局的组。: L6 \( e% q) d
' p! @7 |: h* }/ s1 E4 K+ {P值没有直接给出,但是可以自己计算。& {% d* m3 m! w+ g }( y _8 |
. H" c/ g) j8 t
# 计算P值
( ]: X2 R# G; W9 Tz <- abs(0.001893939/0.027816095)3 Y- j$ I$ S' I' `5 {. `/ X; w4 ]! O
p <- (1 - pnorm(z))*2% y7 c. T3 x& G& j# O9 ^
p
- M) Y1 h E+ H A9 s1
; n/ a. R' q8 q$ k## [1] 0.9457157
* H" m7 P, O; Z: O# E1
# y6 D- N7 M& W# w; u/ tPredictABEL包& i7 T, P: _8 I/ e! A
#install.packages("PredictABEL") #安装R包
$ q* g& A( u, s6 ~2 M+ ilibrary(PredictABEL)
& y4 K% E; i. N5 G; F. t
% v3 Q+ a- F) j3 I3 _# 取出模型预测概率,这个包只能用预测概率计算
+ i+ q7 f7 l2 v* j' Qp.std = mstd$fitted.values5 T/ N$ R3 V- V* o0 n! d* B
p.new = mnew$fitted.values
" W5 ]3 l$ G* z z7 |5 X1 N3 |4 E5 F( w; x
然后就是计算NRI:% q- z e% O) C. N6 `: ?
* F8 |5 k4 I- J S) p
dat$event <- event
3 I# \0 ~ P$ f' v' ^
) \# {% m% P# h& zreclassification(data = dat,* _6 ?+ |# G* ^* }+ O) Z+ G+ }
cOutcome = 21, # 结果变量在哪一列! J/ L% d4 Z6 f; k
predrisk1 = p.std,% q! {: ~. s7 E' M# v
predrisk2 = p.new,
9 @) x N* C; u$ N+ L7 D cutoff = c(0,0.3,0.7,1)
& ~' o$ n; F, e, `) [8 J: @ )
8 m* s( \& B' x- e) ~' O* E2 m9 H$ ?1( W$ R- R& D9 U
## _________________________________________
/ t% A- F1 z/ Q2 p2 F3 r' J##
& S- R/ T; Q' ^0 U4 I) p" \## Reclassification table $ u# _+ h) l0 L4 b R3 N0 @
## _________________________________________ y! S6 B# ~; k4 w, [' p* _( n3 F' Y
##
: |- `* G& g8 o- G2 r( \## Outcome: absent
' P4 n8 z* R4 s+ r## 5 m' P8 }4 X+ Y# V
## Updated Model
4 P) Y( V$ g% ?/ X& g: X## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified2 b4 K4 ]4 B- o; |4 _9 F
## [0,0.3) 121 4 0 3
4 a' ]# w8 Z5 t& S- l# z## [0.3,0.7) 1 13 1 13
# X) c+ ]4 U3 f8 s## [0.7,1] 0 1 3 25
4 `% j/ R1 y# n: U## - [% e$ b+ ?5 h C
##
9 B$ F, B" U0 h2 a* B## Outcome: present 5 T- f+ y. i) F8 T4 m: `
## : @2 D9 F" ^ z$ Q3 c1 Y3 w
## Updated Model
9 |: D S, c8 r$ p1 ^6 S8 W4 g## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified" q+ U' i2 L8 T) Z
## [0,0.3) 14 0 0 0
& H' d% W1 z* K e! X+ E## [0.3,0.7) 0 18 3 14
0 s) v- Z9 u/ A7 V0 A# g: l, ]## [0.7,1] 0 1 52 2
8 [7 O5 G+ w* [6 Y' |! y$ v##
6 Q# f/ [3 e, ]7 A' p0 Y! r) G## 6 @$ _* ]8 O1 A: `/ x5 b% ~
## Combined Data 9 }* f9 V3 g+ ?9 T$ S: `, f
## - z# n$ Y0 C! g
## Updated Model9 u( t7 |" L. X
## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
' A; R' i0 X; m1 u8 h* {## [0,0.3) 135 4 0 33 d6 p5 L* H$ P# F
## [0.3,0.7) 1 31 4 149 w7 y# e! V6 \/ G* }
## [0.7,1] 0 2 55 43 s4 w9 h- x: B" K2 ?3 J
## _________________________________________5 ? X+ Z5 B1 n' I- W
## ' w" W7 x2 d g
## NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806
6 q6 N X% l9 _5 I## NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048
6 ?# p9 ?# h3 N## IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.28396
) g9 a" `& y4 c* p) ?/ w
( p3 R/ U2 L1 @8 k. d1/ Y9 x: F8 r+ A& D7 b$ G) F
结果得到的是相加净重分类指数,还给出了IDI和P值。两个包算是各有优劣吧,大家可以自由选择。' M4 _7 [* j9 r* C1 p
7 C+ R& h* B& @8 E生存分析的NRI
4 S2 S* t0 X! v' T( y3 @还是使用survival包中的pbc数据集用于演示,这次要构建cox回归模型,因此我们要使用time这一列了。
2 @; K' i$ ~% i+ P, [
4 X9 l+ A9 l2 c( _% C8 q+ Onricens包! t+ P# e) g/ J9 v4 K; T
library(nricens)& y4 Y0 \3 |! h& l# o
library(survival)
/ w9 b' i' R8 }, B U9 @, V
# }6 @, l, F2 c8 D$ Fdat <- pbc[1:312,]3 T5 |- `! k) q; C3 `6 s+ ]
dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡
; U. a( N2 u( K! I1
# t' c1 ?) f9 f* G5 }6 C. [然后准备所需参数:
4 x$ E/ K. ~( v+ f+ a
, d- _9 u; n! J( V" O x x4 I# 两个只由预测变量组成的矩阵
; v4 |' x$ r& S" p) rz.std = as.matrix(subset(dat, select = c(age, bili, albumin)))1 L* L7 b( \9 D
z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))
/ h/ v0 X2 N! J( T V" J8 y
) @9 Z) M# Y4 p2 ?# 建立2个cox模型
+ u5 [. {; F- F1 r/ ^% u$ O) Dmstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE); l! C8 F+ F7 u; |' m. [& g& a
mnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE)
% o# U* A6 B/ U$ h( j$ b
8 a) T) U4 P& {) s- N" }% V# 计算在2000天的模型预测概率,这一步不要也行,看你使用哪些参数
" n5 p6 B: J* Tp.std <- get.risk.coxph(mstd, t0=2000)
) H+ Y: w; [9 `. ?( Qp.new <- get.risk.coxph(mnew, t0=2000)
7 @# n+ q' j s( T. e3 T! [2 z1
% t2 O# j) N0 D! n+ K0 J2 _1 e+ {计算NRI:7 j6 U* b0 m$ j& |
; o8 l9 q( x3 ^ `
nricens(mdl.std= mstd, mdl.new = mnew,
! |$ Y8 J3 n# Z0 ? t0 = 2000,
# k. F& ?& }1 m! a' S+ x cut = c(0.3, 0.7),
; u# Z) \2 B3 r9 p) { niter = 1000, ! C: j2 Q: o4 ]1 U+ r% V8 y8 l
updown = 'category')
3 J' Y) l: [, B4 k! i4 f2 _# p: Q* Q& Q3 e3 N; q% M
UP and DOWN calculation:9 \ A }0 [, P5 T" e1 l' {
#of total, case, and control subjects at t0: 312 88 144* ]8 L. g/ ]% D/ n
9 J$ b7 i1 H1 y5 }8 \% N
Reclassification Table for all subjects:
2 l) W: a1 g( P4 p5 O8 M; W- Y New
0 ?6 P% ~. a( d1 zStandard < 0.3 < 0.7 >= 0.7
' W+ x2 q9 N/ @: H' L6 P7 e* ] < 0.3 202 7 0) G9 C1 r5 a$ N5 i8 x0 u% \
< 0.7 13 53 6
* x# p' _+ @- n% F; n* I' z3 ] >= 0.7 0 0 318 R1 ], J# N; V, `6 F( U) v
4 h0 P% M& U! {( G Reclassification Table for case:
7 u( `9 O3 e% ~' o- v1 x s; b New
" D$ B. }3 d+ S7 DStandard < 0.3 < 0.7 >= 0.7
# E! G; \6 B: K3 G < 0.3 19 3 0
0 C: d/ n" r( G- O! F" N < 0.7 3 32 4! v; ]$ A$ }- A, }7 G* W/ F
>= 0.7 0 0 278 x3 G# L" _+ V; V5 ?- I8 n
. V* X. T5 _- \, b& j7 f Reclassification Table for control:) O4 n6 r6 r$ c' O0 {2 u! V5 O% Q( y
New* N; t7 A( E1 ?2 E
Standard < 0.3 < 0.7 >= 0.71 z) l, G6 ^" P; I+ Y5 k
< 0.3 126 3 0
" e3 T3 `- l" {0 M < 0.7 5 7 2. l. S }" ?4 |0 `# l0 p
>= 0.7 0 0 1
+ ?: h1 P: l2 h; e7 U9 l% Q% e; @" M* P5 k( D
NRI estimation by KM estimator:- @3 a$ O: e& x" t) _" B7 f
& D7 H: r4 f8 f4 \Point estimates:
! e8 z" Y9 k: x0 T$ Y! c Estimate
l, g/ U( c4 _% WNRI 0.053776357 X) _( M+ \: P
NRI+ 0.037486605 B. {$ ?- Z" e) U, e; G
NRI- 0.01628974
2 s. a# `" D) j! _# RPr(Up|Case) 0.07708938
2 j$ k/ F- C, e7 y! L! L2 IPr(Down|Case) 0.03960278
% E9 @" _/ f7 p8 L3 ?+ p4 ^) `5 @Pr(Down|Ctrl) 0.04256352& U5 |% H5 g+ c, O, u: C
Pr(Up|Ctrl) 0.02627378
. T$ g& h2 s+ Q2 n# L' \0 h2 Q5 u- M; O4 p6 T5 ]" p- v0 j( z
Now in bootstrap.." _* l8 C6 ?# q: y" |8 N, w! h
+ c" g) @) E; Z0 k' {Point & Interval estimates:1 Y. q" L% k6 ~, v
Estimate Lower Upper6 ^0 }/ r: n. n: r
NRI 0.05377635 -0.082230381 0.16058172) O5 h' f9 S' H' D
NRI+ 0.03748660 -0.084245197 0.13231776
5 }5 b" W5 k* Y/ m: N' p- B, ^! XNRI- 0.01628974 -0.030861213 0.06753616
" n0 F2 P; v. x$ d ], u2 O: QPr(Up|Case) 0.07708938 0.000000000 0.19102291
6 p. p! M$ P3 O2 r4 ]9 CPr(Down|Case) 0.03960278 0.000000000 0.15236016- c* f+ R/ A4 c2 x" ?: ]# n8 Z
Pr(Down|Ctrl) 0.04256352 0.004671535 0.09863170
F7 m* t0 M! `2 QPr(Up|Ctrl) 0.02627378 0.006400463 0.05998424
( E" S- ~. `% ^% g
: J" ]/ q. W! K8 n% d" b1 `- {9 |7 u; d
X8 h6 ]" {& g: O
Snipaste_2022-05-20_21-49-38
3 I% u! c9 q) P( |" I结果的解读和logistic的一模一样。4 u3 B# s0 d4 A+ A# U! A
# {) n0 B8 b5 N
survNRI包; {4 H- R0 y3 Y1 Y1 {' X. P! r
# 安装R包7 x$ y- v! p# P
devtools::install_github("mdbrown/survNRI")
3 |/ e2 p- o' G" k1 t$ ~% u9 }) F0 X1
3 D1 L+ |7 U2 Z, C# `! i加载R包并使用,还是用上面的pbc数据集。$ I: {7 j& Z, d. S, S7 t" B
& n9 y: o7 `& ~9 i2 d4 F' tlibrary(survNRI)% o- O) f* l) c( } H& i
1
/ u. A8 G4 l n## Loading required package: MASS
$ U4 H4 N8 c Z' g0 D( O! d5 Q& M19 O% K# U& q( C( S% ?
library(survival) _- m6 `# g3 I( Q7 a r/ }+ C
3 o7 q$ `- P3 \& h! F) Q2 A6 x+ {# 使用部分数据1 [( [& e6 a" x0 e
dat <- pbc[1:312,]1 i7 L) Y7 X5 ]: L0 j
dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡/ r Y7 d: F# o
( G! }; ^% u* _& b) ures <- survNRI(time = "time", event = "status",
2 k7 H: j' K& x: o. o: y0 m model1 = c("age", "bili", "albumin"), # 模型1的自变量
! Y" h& a/ }& o- }+ K0 x3 O model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量
" @: U* G- V( S# D4 _7 G data = dat,
6 M3 Z7 f' Y8 n0 r+ ] predict.time = 2000, # 预测的时间点
4 s: B4 B6 d3 _" c+ w method = "all",
5 e5 a V# i# }/ p" L% r bootMethod = "normal",
6 P8 w/ o5 y5 r/ [/ Z& k8 O bootstraps = 500,
( I3 y) @; k' B: V2 G4 a4 _ alpha = .05)$ l+ w0 z; }3 m( t% l) h
. d' m! K& |+ M, |$ Y; R! M1# W! G0 F& f. u! I7 ~6 |9 \' u/ r' R
查看结果,$estimates给出了不同组的NRI以及总的NRI,包括了使用不同方法(KM/IPW/SmoothIPW/SEM/Combined)得到的结果;$CI给出了可信区间。
) ~ L, L8 m- P% Z! F% ]7 F, u! {( ?. M {4 @; N: ?
res1 v, `# c* A+ _. l+ K
18 M/ ]7 X& H* \2 N
## $estimates9 r" k3 k: m$ C: n) p& \; n0 ^/ b
## NRI.event NRI.nonevent NRI
/ a* |4 a9 a2 i5 r7 X## KM 0.20445422 0.3187408 0.5231951- }0 I. r* Y5 S0 a l1 E6 m# w* }
## IPW 0.22424434 0.3273544 0.55159872 I2 _. O8 W; ]6 l4 p7 _ W# \
## SmoothIPW 0.19645006 0.3144263 0.5108763
$ r, w* f9 K2 V9 s, k8 d4 I## SEM 0.07478611 0.2632127 0.3379988
/ h% Q' f5 I9 `- G. T' m3 T& H, V' ?## Combined 0.19633867 0.3143794 0.51071812 I4 o: s, D [- [) ^* m) n+ j# n8 ?
## . m; M6 T* y. Z3 d D
## $CI2 M& w9 _: j% n9 m4 s5 w
## $CI$NRI.event
; u# M* k; @9 H+ i## KM IPW SmoothIPW SEM Combined8 j- K* j; V; b* I. q! n
## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723& E s5 b7 I, h& R4 Q9 [' ?8 n
## upperbound 0.44806768 0.47033936 0.44014214 0.2658309 0.4400496
/ u) u' `0 T" f6 P. |) ?, W##
5 K% j6 D: d) `! D7 [% k5 K## $CI$NRI.nonevent& x5 w2 q4 O' \7 r) Q
## KM IPW SmoothIPW SEM Combined9 g: o7 ~# @% v
## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.1286426
% p: P0 r7 w2 D0 @2 G. d## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.6964549' t) ^+ s8 K; Z
##
, M# a8 g( p8 {. p1 ~## $CI$NRI3 a9 u$ e+ I* b s% ~
## KM IPW SmoothIPW SEM Combined$ t, n2 a3 ]; r4 E# ^
## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409
8 V. S8 j! b6 l5 T1 B6 c## upperbound 0.89306122 0.92464359 0.87970125 0.64253510 0.87953153
# _3 y: D; [1 `' E# F8 Q##
, p1 I. z5 n3 h##
) d( g" A7 i; I. O& N% C## $bootMethod
) L$ P) r9 ]3 t8 w) M## [1] "normal"
$ m; A3 a8 G5 z0 S## 6 q/ q/ R+ F7 V+ e8 z8 Q. x1 K
## $predict.time
: W* X& u& k+ }! i1 Y* @## [1] 20004 o ]7 B. B# t3 e, p1 B7 W% |+ }
##
, T8 m3 ?, R' M## $alpha
% F7 c4 P) V7 o- `* ?6 Y2 M## [1] 0.05
. m( X; K8 h. Y# G0 D## $ T, } w/ `5 |/ I
## attr(,"class")
& U' R; y3 o" ~) b& U: X## [1] "survNRI"
! z! b; ]0 r4 X: z2 ^3 A$ `
5 I9 e4 {6 o n# w/ V1
3 n. \+ a. L( l) MOK,这就是NRI的计算,除此之外,随机森林、决策树、lasso回归、SVM等,这些模型,都是可以计算的NRI的,后面会继续介绍。大家如果有问题欢迎在评论区留言。) w# D) b4 U# d4 R9 G' k
' H5 c6 R% f) c. a, {
本文首发于公众号:医学和生信笔记% v- P2 q4 N7 S# W
, |' }9 F9 r6 s* k8 j2 d4 S
“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。2 n' `7 o& k+ h
本文由 mdnice 多平台发布6 U: O- o" z8 c8 q& j! o
————————————————; d4 Z( k, m! b
版权声明:本文为CSDN博主「阿越就是我」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。+ Y- f4 X2 ^# i, _9 z% ]* q1 q) h4 k
原文链接:https://blog.csdn.net/Ayue0616/article/details/126768006
! _1 X% X: Y+ H# P+ R1 l1 S6 c) {1 B3 E7 a( E! C
+ @* |$ G1 e8 b7 _7 R5 K, ~+ U* R9 T9 D |
zan
|