- 在线时间
- 1630 小时
- 最后登录
- 2024-1-29
- 注册时间
- 2017-5-16
- 听众数
- 81
- 收听数
- 1
- 能力
- 120 分
- 体力
- 543261 点
- 威望
- 12 点
- 阅读权限
- 255
- 积分
- 168357
- 相册
- 1
- 日志
- 0
- 记录
- 0
- 帖子
- 5324
- 主题
- 5250
- 精华
- 18
- 分享
- 0
- 好友
- 163
TA的每日心情 | 开心 2021-8-11 17:59 |
---|
签到天数: 17 天 [LV.4]偶尔看看III 网络挑战赛参赛者 网络挑战赛参赛者 - 自我介绍
- 本人女,毕业于内蒙古科技大学,担任文职专业,毕业专业英语。
群组: 2018美赛大象算法课程 群组: 2018美赛护航培训课程 群组: 2019年 数学中国站长建 群组: 2019年数据分析师课程 群组: 2018年大象老师国赛优 |
; v: H) v' j' E
净重新分类指数NRI的计算
4 q$ |9 B5 x! w7 D; \' }! B, n; t“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。 R! P* d' T8 @. s4 a' K
NRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!0 u, ^8 z8 i$ y5 }6 m
! D5 m; \; d+ z& ?$ f% P% d在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens包,PredictABEL可以计算logistic模型的净重分类指数,survNRI可以计算cox模型的净重分类指数。. t9 \% q7 |; a5 O) _3 B- I
0 z# b. {6 j6 }! k" d
logistic的NRI
/ R5 k( k7 e* O. w) ]: {3 X: A! s/ Qnricens包 x' Y' y+ s' Y! s
PredictABEL包
) V y1 S7 V& }$ t生存分析的NRI
7 `1 N; Q* M2 o: m% Dnricens包6 e/ B- M, Q% v" n% `* y; |& C
survNRI包
& C* Q* B N5 n; ]0 C& Mlogistic的NRI4 M0 h8 y$ z! ^( d! D, e+ U
nricens包. ]' g w: L7 V& R) ]" w9 p# r
#install.packages("nricens") # 安装R包4 J" A+ X3 ?8 K1 `7 ]# X
library(nricens)" P% P( n' O' W- O; Y" T
1* t1 c4 a! a7 C
## Loading required package: survival2 E8 Y. I: K6 K8 X% r
1
9 }8 w5 c( g( t: [使用survival包中的pbc数据集用于演示,这是一份关于原发性硬化性胆管炎的数据,其实是一份用于生存分析的数据,是有时间变量的,但是这里我们用于演示logistic回归,只要不使用time这一列就可以了。
0 L/ h3 g+ z/ Q6 f- ~ N9 {2 O5 ^4 c7 p. Q
library(survival)+ z7 z$ q8 N9 ]( w5 X. ^" }
* M% q+ a3 F' T; b, u( [9 q
# 只使用部分数据% P: U2 R! I) D F$ x5 z/ ]
dat = pbc[1:312,] ; x+ x4 v3 } d1 ?/ i) [
dat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]) V* i3 f" n' u' J. G& v0 x& D
& }. g6 c! g1 r2 r
str(dat) # 数据长这样: M/ s3 z! F. y. x; i2 R
1: C( e+ i P8 o1 o+ [
## 'data.frame': 232 obs. of 20 variables:
- O3 n2 u, @% ^% h## $ id : int 1 2 3 4 6 8 9 10 11 12 .... Z- V( f3 P V" x7 V; h; x
## $ time : int 400 4500 1012 1925 2503 2466 2400 51 3762 304 ...+ L! u% }0 k$ b7 G! ^% c. V$ K
## $ status : int 2 0 2 2 2 2 2 2 2 2 ...( `0 a% F3 U' b5 |2 l5 n4 j2 l- Z2 E
## $ trt : int 1 1 1 1 2 2 1 2 2 2 ...- U3 ^6 t) i+ N$ J: {/ k9 x
## $ age : num 58.8 56.4 70.1 54.7 66.3 ...
/ z% p Y6 u3 g0 [## $ sex : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
8 d( ]4 `( P* s. N## $ ascites : int 1 0 0 0 0 0 0 1 0 0 .../ v! M% Z' t5 f) z( x3 T! E
## $ hepato : int 1 1 0 1 1 0 0 0 1 0 ...
" a; [# Z: a% W/ Q) N## $ spiders : int 1 1 0 1 0 0 1 1 1 1 ...
9 g/ \/ S9 t/ E8 e+ ?: U& }## $ edema : num 1 0 0.5 0.5 0 0 0 1 0 0 ...) E+ p9 L- d2 d" q
## $ bili : num 14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...1 `% X! Z- c; b
## $ chol : int 261 302 176 244 248 280 562 200 259 236 ...( Y7 e" I3 Y# p6 S* A
## $ albumin : num 2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ... g) z- X# W6 G, r. b
## $ copper : int 156 54 210 64 50 52 79 140 46 94 ...
( N! g+ r9 K2 L2 d/ o## $ alk.phos: num 1718 7395 516 6122 944 ...4 Q% u6 s; f, C( B
## $ ast : num 137.9 113.5 96.1 60.6 93 ...+ h1 S+ O* R. Z3 C( f
## $ trig : int 172 88 55 92 63 189 88 143 79 95 ...
) g8 @3 z! o5 n3 c% F## $ platelet: int 190 221 151 183 NA 373 251 302 258 71 .... |( f* Z; Z+ T: p7 d; [; @
## $ protime : num 12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...0 n. j H* u5 x
## $ stage : int 4 3 4 4 3 3 2 4 4 4 ...
& n1 d7 H: D. }1 g3 o
+ ?4 t. t5 z8 s( P' D. Q1# }% ?/ c5 t* m( ]$ L
dim(dat) # 232 20
* x3 Z4 \" Q3 t# l1' R! l) R, T( o3 y C
## [1] 232 20
) \2 F3 J8 {# U0 F, \1 {* x% s1: {6 D2 I- A9 x }. W j
然后就是准备计算NRI所需要的各个参数。 M. I5 ~& ^3 s3 B/ F1 U8 q& @
9 r: W& L" v- F& Y; V
# 定义结局事件,0是存活,1是死亡
$ F- k/ G' l( N1 R: D3 Oevent = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)
8 d% X" w. H |) t
& b6 B d7 w' G0 j* C7 l# 两个只由预测变量组成的矩阵
7 u+ _) n0 K+ P5 Wz.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
; `! n: Y" d0 O Xz.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))
5 U; t- b' B4 `% N3 N5 m8 k
& c) o! H; n" i+ K# 建立2个模型3 v( s6 a9 ]" h7 w
mstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)7 C5 e7 C/ u: C4 |* `! e" M* \! H1 M" k
mnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)
b' c6 Z* p! g6 X, _; c* M4 F, X+ r. ^( \1 }
# 取出模型预测概率
% Z$ M& Z: C. [4 S# Ip.std = mstd$fitted.values0 h; l& i: G8 ~! W
p.new = mnew$fitted.values
9 b& r! p+ j4 W* p3 K ~3 r0 A. [( C1 O3 d( u9 j8 p& I+ h2 X
1
3 U9 E. [/ T4 m然后就是计算NRI,对于二分类变量,使用nribin()函数,这个函数提供了3种参数使用组合,任选一种都可以计算出来(结果一样),以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。+ M+ C7 V' i# A
3 A D$ c* n& M0 M: l/ W; z
# 这3种方法算出来都是一样的结果* o8 v6 J( c$ Q# h
- W) F& C5 Z* i7 { J# 两个模型
" l6 {$ H3 J; q6 i2 j1 [nribin(mdl.std = mstd, mdl.new = mnew, 2 q2 w3 g) v& x& W3 }
cut = c(0.3,0.7), . }- P O1 X8 q- k1 C* F, k
niter = 500,
0 q( [! M4 p4 T0 f( Z- a4 o updown = 'category')& @! P8 z2 O( ]# S- r) o
" I" w9 A3 c+ O. [! t! N
# 结果变量 + 两个只有预测变量的矩阵
6 I, @ p3 x2 p, W) w" r. u- Onribin(event = event, z.std = z.std, z.new = z.new,
8 G6 d( o( X' u7 G3 Z X, ?- |: k3 g cut = c(0.3,0.7), 0 z5 Q+ E! b: D
niter = 500,
8 e3 t2 @7 o4 L/ s( F9 R9 ` updown = 'category')5 K( c8 k" Z& m& ~2 b: Y$ ~
! {8 K% b/ g: n2 U## 结果变量 + 两个模型得到的预测概率
% N+ X7 `2 J$ unribin(event = event, p.std = p.std, p.new = p.new,
% X+ E3 [' V. w. Y cut = c(0.3,0.7), 6 x- D0 ]+ X$ Q6 A; \
niter = 500, * F- A1 z, P' S& R* X
updown = 'category')
0 z+ P. n$ w1 L; E/ X1 a k+ }0 d3 T3 X$ c9 P }9 N2 R0 l
1% L. z3 o6 e+ l: u s4 h
其中,cut是判断风险高低的阈值,我们使用了0.3,0.7,代表0-0.3是低风险,0.3-0.7是中风险,0.7-1是高风险,这个阈值是自己设置的,大家根据经验或者文献设置即可。0 O( k" g% h/ n2 V
' l- J, c* m- B9 U: I: e
niter是使用bootstrap法进行重抽样的次数,默认是1000,大家可以自己设置。
5 ?4 {; z0 @* ]7 L
' v" t6 T5 ^1 cupdown参数,当设置为category时,表示低、中、高风险这种方式;当设置为diff时,此时cut的取值只能设置1个,比如设置0.2,即表示当新模型预测的风险和旧模型相差20%时,认为是重新分类。) J6 P: s8 G( s& ~7 T9 G \
+ i6 z* q4 m0 Q5 N上面的代码运行后结果是这样的:$ _% }3 d/ ]( S# Z" v! @6 _
6 Z& w$ Q0 D9 {$ I8 HUP and DOWN calculation:& f$ ]) g% B3 A7 y9 T
#of total, case, and control subjects at t0: 232 88 144% ?% S( y+ b% V+ Y q T* `
! X$ V* h" L# z$ Z Reclassification Table for all subjects:% S* D, [3 H6 f4 d1 U2 a7 ]
New
: C5 j9 Y. ~& mStandard < 0.3 < 0.7 >= 0.78 n2 P( t0 Q8 F0 g
< 0.3 135 4 0
5 z4 j9 n" Z; F < 0.7 1 31 42 V9 k8 `% S6 S4 S. R, ?; C7 F
>= 0.7 0 2 552 q: M; r& N- c3 p
7 r1 E6 p1 d4 _) ~3 Q# `% W
Reclassification Table for case:# z' u4 c, I* u3 Q: d9 Y
New9 L: Z, a! [0 j- l# q
Standard < 0.3 < 0.7 >= 0.7
! U ]# ?2 j( s, S- Q < 0.3 14 0 0/ t1 f" S# U9 A; Y m
< 0.7 0 18 3
$ L5 r% @0 P, @ S3 n >= 0.7 0 1 52
" V. T3 r% I+ \6 F; W1 y
0 y! @& B# K D% v# L Reclassification Table for control:( z' c; ?. C8 \' U- X h( N) M/ {
New4 C8 d7 {2 e- O/ P
Standard < 0.3 < 0.7 >= 0.7
; ? O; r$ k0 N3 G5 l5 T: X+ J < 0.3 121 4 0) b) a8 N' x5 g. n7 w
< 0.7 1 13 18 J" p0 m! y2 a! W; l2 X
>= 0.7 0 1 3
- q2 r$ ?: T w
% z: k, v" h, x- F+ q) z6 oNRI estimation:
" g; i/ n/ q1 a' M4 I2 u! YPoint estimates:1 b0 M4 x, Z- E& M$ M( j
Estimate
9 ?' Z3 R3 F$ V, L. u1 A6 y8 {$ PNRI 0.001893939" A2 Q# V% |8 G0 s
NRI+ 0.022727273
, h s3 M4 z8 `: W3 rNRI- -0.020833333" D# N8 S6 y# j7 B
Pr(Up|Case) 0.0340909091 m" S' K6 ~- k0 Z& B( _
Pr(Down|Case) 0.0113636361 z3 `5 p# G5 u6 M& K. X' Z
Pr(Down|Ctrl) 0.013888889( t+ ~4 y4 ^- u4 A7 U
Pr(Up|Ctrl) 0.034722222
& G5 |, L. y1 x. h$ M, u8 ^; H: c+ d. [
Now in bootstrap..$ C- e4 I1 }/ O; ~
4 c1 V L- ?$ a
Point & Interval estimates:' K5 Z& ]/ V9 m: u! F/ A o
Estimate Std.Error Lower Upper
! _5 D% a0 ^9 K. ]- J- eNRI 0.001893939 0.027816095 -0.053995513 0.055354449 r/ g: u6 s: E, }/ S4 [
NRI+ 0.022727273 0.021564394 -0.019801980 0.065789474
" f0 d- s( p1 V# rNRI- -0.020833333 0.017312438 -0.058823529 0.007518797/ {# V K: R* x4 x6 F8 K
Pr(Up|Case) 0.034090909 0.019007629 0.000000000 0.072164948
" i! X( h0 M* e& T; uPr(Down|Case) 0.011363636 0.010924271 0.000000000 0.039603960
+ {9 w& E8 R5 \2 J2 K: ?$ T! [8 [- HPr(Down|Ctrl) 0.013888889 0.009334685 0.000000000 0.0352112681 e8 x, |/ j) @2 G! H3 U
Pr(Up|Ctrl) 0.034722222 0.014716046 0.006993007 0.066176471
, G# D: t6 t3 V W/ ]6 Y% G2 F, r3 s* x6 w$ @8 r! P
1
, A# D; }. e4 E7 n首先是3个混淆矩阵,第一个是全体的,第2个是case(结局为1)组的,第3个是control(结局为2)组的,有了这3个矩阵,我们可以自己计算净重分类指数。
4 {4 p( l- e1 K/ @5 X0 W, p4 u5 l6 y) a/ ~1 ~6 v3 z
看case组:" J+ K) A2 Y3 C$ }7 E6 I
8 c' h! T. ]. Z; V8 }' h' f7 W
净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273
) e+ O( B$ S( ~0 u! ?* q
4 V' `8 ~9 p- i( M; e再看control组:7 i1 D' [" |4 P, {$ l, x: d: m
- e4 H1 x- c8 L8 [: C净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.020833333
4 k' b% s! m$ K0 n" B, L
{# Z% Y K7 I' K相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.000315657) b! [$ H8 L- h, Q1 ^
, a% T6 X- U8 b: T) ?" E0 d1 u# }. j6 Y再往下是不做bootstrap时得到的估计值,其中NRI就是绝对净重分类指数,NRI+是case组的净重分类指数,NRI-是control组的净重分类指数(和我们计算的一样哦),最后是做了500次bootstrap后得到的估计值,并且有标准误和可信区间。! ?! v9 R# U8 R2 W( H& y# i! `: }' \7 R
3 v* O/ k. ~' A最后还会得到一张图:9 B; d* R1 u; F
2 ~$ \4 ~& M! x4 f
这张图中的虚线对应的坐标,就是我们在cut中设置的阈值,这张图对应的是上面结果中的第一个混淆矩阵,反应的是总体的情况,case是结果为1的组,也就是发生结局的组,control是结果为0的组,也就是未发生结局的组。
2 Z& G9 V; w p# d) ^/ D( w/ M7 ?8 `
) F5 v" A. t9 WP值没有直接给出,但是可以自己计算。
) ^# M9 h$ z7 O% B5 N
% C! ?+ F/ b ^" |% E# 计算P值 o: V. P* k4 q, B# S
z <- abs(0.001893939/0.027816095)" v# p6 d% C5 k$ M' `, u; ^; U+ z/ u
p <- (1 - pnorm(z))*2' E2 p" S& @" C
p( b. e0 k2 b/ J5 `) y
1
; @1 _, K6 N) O## [1] 0.9457157" C1 y2 ^: `) K5 q, v; o
1
6 @0 s) c& [% Z x3 r+ mPredictABEL包8 I6 A' B$ c/ B& v4 k" U3 F
#install.packages("PredictABEL") #安装R包
# T9 M. R, s4 v5 `1 a* i& wlibrary(PredictABEL)
F' O8 R, I9 o" i9 l2 [3 u; n. i5 V2 p& @( }# P6 I7 ]
# 取出模型预测概率,这个包只能用预测概率计算/ z1 P9 R6 _: }! H! B2 M+ Q
p.std = mstd$fitted.values" I! j. ^/ c# ?
p.new = mnew$fitted.values ) I4 y4 R$ x" \6 p: t( ?) d/ [1 r
1
$ w" e: [( [% [( l; h3 Y然后就是计算NRI:* }' W. `- v* ~+ |% U4 c+ m
" X6 Q5 W6 L+ f. E) s/ P$ c- |: Wdat$event <- event
. Y, t+ ]" w3 e4 Z2 D1 Y5 U7 K0 u( D" d; q3 n
reclassification(data = dat,% [9 ~# A# m& O% l- h. y
cOutcome = 21, # 结果变量在哪一列
+ H, q+ ]5 L* ` predrisk1 = p.std,: A4 q/ b; z) Q0 o' j! T. ?
predrisk2 = p.new,. H' v3 U) w9 l/ y/ A( [$ P0 c g
cutoff = c(0,0.3,0.7,1)0 U% d0 \' m8 ]+ X6 {7 I
)7 [' `+ v4 J) I) G" H( _
1, ?5 Q- b7 ?4 J% i$ R5 \
## _________________________________________
Q& z0 \. \3 x9 t##
5 i2 J$ n6 G% Z## Reclassification table , M2 i7 B [$ Y% p% C/ q* n
## _________________________________________
: w. i T* T& L##
1 p, K1 H- h. D; h t4 ~4 c f3 V2 `5 w## Outcome: absent " Z; V" S/ O, y2 `; T: `% x
##
) \( F7 D' k7 e# U## Updated Model5 x t' y& O' V. J4 O n# n
## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified; T* S. D3 s& f1 ^4 o0 q1 Z `: m
## [0,0.3) 121 4 0 3
6 R2 v, w# S4 s5 w" j## [0.3,0.7) 1 13 1 13. l q C6 R. y+ |- ^6 {4 P
## [0.7,1] 0 1 3 25- N; s9 {: H1 l# a% B1 `4 {3 A
##
3 o) g6 }7 ?8 q Z##
# B3 z. ?" Q* i: Z$ b# N## Outcome: present
: v0 v# [+ w& y0 o: {##
& H/ @) p; n) X' ?## Updated Model
8 G$ \) O G- R, v9 D## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
. \0 n6 `( w8 c0 o0 n## [0,0.3) 14 0 0 0
9 E2 P5 U7 D/ D1 q## [0.3,0.7) 0 18 3 14+ P) b- R: B5 s& w! D
## [0.7,1] 0 1 52 2
R+ Z/ F# A# j## ' H: x, F" r" k5 g0 m
##
! i# @0 g) _5 ?; {## Combined Data R: G* o" ? V4 {! `
##
" O) V5 Y& g) y, [. E## Updated Model
5 a: d9 T* @. @: N- w3 A4 q## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified3 `3 p x9 | |* \( Y' T6 @! M
## [0,0.3) 135 4 0 32 ]5 u, U" U# H$ ]& e8 A+ H
## [0.3,0.7) 1 31 4 14
# J* h6 K) \: p% X## [0.7,1] 0 2 55 47 {% g2 ]9 R# K o( L' N3 v; V
## _________________________________________7 |; e3 K \+ k/ v p4 N9 h% E
## $ F6 _' b- Y* H$ u9 F
## NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806
: c2 P& s2 _0 e1 P4 T" ^## NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048
4 i/ g5 r4 A m- g7 N5 W## IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.28396/ ]- f! ~' T1 U* V' m: N/ g/ t, I0 C
5 q% a) ]7 R& R( ?9 g ?# e
1
4 W+ W( D1 P' B# X. p2 |$ j5 u( v结果得到的是相加净重分类指数,还给出了IDI和P值。两个包算是各有优劣吧,大家可以自由选择。1 |2 `4 O9 s0 g3 ^2 R' U
( V1 o3 U) ]* g+ P, J; A2 e
生存分析的NRI
# d% X( `- P% I* i( k还是使用survival包中的pbc数据集用于演示,这次要构建cox回归模型,因此我们要使用time这一列了。$ P: M7 s' {. S ^8 d; i6 D
/ L% ]9 T1 S& }' snricens包
* A3 {3 ]0 c& M2 X$ f0 _0 Zlibrary(nricens)
9 ` c2 q! `6 @+ x; b- Dlibrary(survival)
9 y; }- P' ~# C; j/ u
4 d) \/ B( V* Y) n3 |- _0 e, Ndat <- pbc[1:312,]
, ?( L3 f$ R- J% c- I; Xdat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡' }7 o5 o) w$ D/ x5 v4 F$ W$ _& j6 c
1
) A, ?+ F2 s' P+ D+ [, L然后准备所需参数:% M* E: I; B) {$ \/ R! d2 I9 ?
# h5 [& z9 \( {) G. p r# E8 @# 两个只由预测变量组成的矩阵9 ^ p7 b6 s7 ?0 y2 p9 G
z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
d/ V! a. h% {: A# h. V) z% Oz.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))
0 o7 _* C& G# j+ y. c1 j3 `! l/ b( {) u
: L! Y8 I% z% R' k6 w8 Z/ H( L9 b# 建立2个cox模型
* L* Y0 c3 D. H7 G: N, Umstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE)
9 H& d4 q- }8 x, o7 z: |6 Mmnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE)
8 F4 G! _! l( ~; S8 y
' _+ ?: b0 _. c5 g# 计算在2000天的模型预测概率,这一步不要也行,看你使用哪些参数
2 ?1 q4 p- i) a; q4 K+ V% [p.std <- get.risk.coxph(mstd, t0=2000)
# J- `. n" l" V* ^p.new <- get.risk.coxph(mnew, t0=2000)
2 h, w9 \; f- u# k* m1' M/ [- n/ ~' e
计算NRI:( c4 F# j) X) ?, ]% k
3 A0 X5 z8 T; X& |/ c1 [2 N' T8 Znricens(mdl.std= mstd, mdl.new = mnew,
% w$ n1 `0 C" ^3 w) C t0 = 2000,
' o* {: G. r+ P) r' ^" e; D0 }# Q2 i cut = c(0.3, 0.7),
2 d; K: y" \; \2 e" R8 ? niter = 1000, 2 W. c7 Y5 P% V# X% W* O, X
updown = 'category')- N/ T0 {% h' ]% _7 o9 J( z
) j* P6 t6 K2 A0 k: m, ~UP and DOWN calculation:/ d% Q5 K9 m8 V; X2 F
#of total, case, and control subjects at t0: 312 88 144
8 n5 a# e- u9 k" h/ f: t$ H2 g: N" X4 F% T, r1 q
Reclassification Table for all subjects:" P3 |+ ?7 z7 y8 \; g" j% x! ?
New
' K7 A+ `4 O* u0 I6 L% l* N% h6 UStandard < 0.3 < 0.7 >= 0.7
8 n) _) L# _0 S( i < 0.3 202 7 0: }4 ~0 ^0 V/ V' g9 x
< 0.7 13 53 6
3 {$ J6 q0 i' J1 d# k. h >= 0.7 0 0 31: a9 F( z1 L5 @' o
; v, F+ r9 \5 A) ]' Z0 }
Reclassification Table for case:
" m" y9 s' a. z, T7 n3 U: L New( ~( E$ ]) f d8 H5 [
Standard < 0.3 < 0.7 >= 0.7* g' W- T" |! r6 r! g, e5 L, Z
< 0.3 19 3 0
+ x6 O+ j; a' T9 b5 d < 0.7 3 32 4
+ T& [7 ]0 m6 F+ {) A' r" t >= 0.7 0 0 27
, F/ t$ g6 W# ^$ I. z
( A* G8 L3 j( U* j! x: h. E Reclassification Table for control:# O% Z ]2 m$ J. V' S( N D
New* N+ j$ d- M+ q( E7 r
Standard < 0.3 < 0.7 >= 0.7
; n- W+ }+ Q8 j H < 0.3 126 3 0. c% a( L: l- P5 T2 h$ Y8 U
< 0.7 5 7 2$ H: g/ @1 {0 {
>= 0.7 0 0 17 {# V: i. `$ F
% i9 N8 a+ P6 g4 r( g; N% gNRI estimation by KM estimator:, l( D8 ?& I) Q+ @7 q( m
2 t3 H# c' H9 j, s: q4 `: @Point estimates:
% q4 y& h+ }# q% @& a! G c Estimate
4 F! a/ R0 ~/ l; e) ?: wNRI 0.05377635
6 r/ G8 d( R+ I9 ANRI+ 0.03748660
5 r( D& X# f' U0 [$ V8 INRI- 0.01628974
( h! O. a# A5 i$ W5 K* FPr(Up|Case) 0.07708938
# d+ M s5 E7 K; rPr(Down|Case) 0.03960278
d; e8 R: Y0 p4 T( P: iPr(Down|Ctrl) 0.04256352
( U4 t f: Q- _% r# cPr(Up|Ctrl) 0.02627378* r0 j4 R/ |: H& ]( c8 }9 s
; \8 H$ o5 `2 nNow in bootstrap..- Q" Y3 L* [; t! A9 N, V, R
% b5 m- j, M' B7 ?: OPoint & Interval estimates:
! ^/ F3 ?% _8 \8 y% d4 E0 P* { Estimate Lower Upper
2 H0 Y' }! U# o8 f, ?+ S9 c3 E$ FNRI 0.05377635 -0.082230381 0.16058172
0 Q7 [. [; R( t6 H! G7 kNRI+ 0.03748660 -0.084245197 0.13231776
+ v* N2 i2 O1 r2 A3 O+ YNRI- 0.01628974 -0.030861213 0.06753616& r5 F% t6 p9 Z6 `- T' I! F
Pr(Up|Case) 0.07708938 0.000000000 0.19102291
2 P z. W; _3 a$ s( t: ~! sPr(Down|Case) 0.03960278 0.000000000 0.15236016) m9 U# M* [3 w- ` k; w
Pr(Down|Ctrl) 0.04256352 0.004671535 0.098631708 t) n- h( v' D! E9 ?* J; r
Pr(Up|Ctrl) 0.02627378 0.006400463 0.05998424
* b1 N: a+ o5 A4 G- h9 f
& s* ^' k( _# I9 A, @3 b1# f# T, z. f0 {* L2 W
3 T* X9 w2 v+ H
Snipaste_2022-05-20_21-49-38+ e$ \0 V) D# E) z+ K8 s/ g
结果的解读和logistic的一模一样。
3 U* m4 T) n3 A0 k2 g9 d- f2 P- Q- s! K$ T) I# W+ C
survNRI包
1 ^9 f/ D, g7 t' W1 R3 Y+ t# 安装R包
4 w2 H% T- r8 d9 y. A9 j0 l sdevtools::install_github("mdbrown/survNRI")
& w* T$ G- u4 v9 X) D* ^4 e1- f! A6 d; |5 u; Y
加载R包并使用,还是用上面的pbc数据集。/ O' ]- x# y d/ T2 X
) I# n/ {$ L; h' nlibrary(survNRI)* d5 f: r% M$ d1 R& o% r
1+ v0 L# d; C* a" O- H7 g6 N
## Loading required package: MASS4 x! Y1 T& ~& ?# t! u
1) x4 g# B9 U8 @; a" t
library(survival)
$ A1 c) ^7 Y+ y; ]
( x4 x' V8 Z h( v# 使用部分数据& T! e; |+ N2 A8 W! X+ [0 V
dat <- pbc[1:312,]! Y' V4 q. ~. g7 `
dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡/ w7 n# @3 r j0 i0 }2 P
3 s: z! h/ m; b& lres <- survNRI(time = "time", event = "status", 2 ?! r5 x+ V$ r0 _% n; C
model1 = c("age", "bili", "albumin"), # 模型1的自变量) w5 A' l5 R: @+ F
model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量
, V4 K8 D- b7 h; b: M" k2 C data = dat, ; o/ ?2 E$ \. `8 K' u* B2 G! J8 N
predict.time = 2000, # 预测的时间点$ Z0 @& G/ y* l. T+ ^4 ?/ r" j
method = "all", 9 U! W& N, h# S: L" E4 X* l
bootMethod = "normal", 9 L0 N9 ^& i4 ]9 w
bootstraps = 500, & k5 o' \0 t8 X; j( B) A
alpha = .05)
* e: s4 b5 ]$ m2 [3 Y, I" H) D! N2 v2 l0 A
1
1 D: i& B6 l" T; [查看结果,$estimates给出了不同组的NRI以及总的NRI,包括了使用不同方法(KM/IPW/SmoothIPW/SEM/Combined)得到的结果;$CI给出了可信区间。! a& Y( ?. e0 \6 i! _' ?! r( @
: ]# w( ^, y) `* z
res" D3 |, X9 K9 z) F2 z; g9 ]
1
' _6 h) q: j6 {7 o: ?## $estimates
' N% p+ T* \6 h7 _/ C* ^## NRI.event NRI.nonevent NRI
0 F7 N* p- v7 w& `2 o. [## KM 0.20445422 0.3187408 0.52319510 U9 L3 N5 p1 E. [) l
## IPW 0.22424434 0.3273544 0.5515987- P& G3 l0 X! u1 w: n$ w3 a% G
## SmoothIPW 0.19645006 0.3144263 0.5108763
8 y& S; s" R7 A## SEM 0.07478611 0.2632127 0.3379988
% ]( k0 H. @" B! x6 b0 s3 t5 h5 ^## Combined 0.19633867 0.3143794 0.5107181
7 c5 A2 S5 E5 q8 Q, t! R##
3 f# u3 x8 ~7 ]; G* B- T## $CI
+ f: _. U' v% {9 F& F j2 }) P## $CI$NRI.event
* b* j) h, Z4 W- f4 r* A0 ]## KM IPW SmoothIPW SEM Combined
$ B9 _: {& [2 c## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723, ], F- r" m( e& O2 r9 a
## upperbound 0.44806768 0.47033936 0.44014214 0.2658309 0.44004967 `! s, ^- R- n: V0 ]+ H
## . A9 Q* M; v% z3 d' [
## $CI$NRI.nonevent" ^% h" i4 g+ O4 I
## KM IPW SmoothIPW SEM Combined3 X' s7 v, h+ H3 K0 U
## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.1286426
3 z7 n/ K9 c4 L, m; u7 o## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.69645490 _0 }3 X* B* V6 n: Z$ _# ]
## - K. J! d. M7 D6 o' P
## $CI$NRI; x& D1 O; B# b0 F& ?% g2 U# ^
## KM IPW SmoothIPW SEM Combined
; S8 i @' H f## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.054434097 R" t( C) h) i+ b
## upperbound 0.89306122 0.92464359 0.87970125 0.64253510 0.87953153
- i- n6 i+ P) S) V- h##
" ^" ]* v O( Q$ G9 a% \& h9 F##
- e7 D" m- Y4 V- W" ^## $bootMethod
2 [$ ^ v+ e4 B2 @4 M; U## [1] "normal"
# R5 |; p* M7 ~5 f##
6 f% A: H, u$ W( D## $predict.time
2 L$ r, o7 H+ B## [1] 2000* |5 [9 R' K' A0 n. h
## / i8 ]$ O: U9 a
## $alpha
$ J* v0 Y7 B) R8 q. h5 w4 _' d## [1] 0.059 o. q+ Q) `+ _# D5 C: I1 C4 S
## 0 [3 O4 Q: F3 F3 o
## attr(,"class")
3 G p) s/ s' c& L) T## [1] "survNRI"' ?) s3 B! e: S: O* X' I
# {, P" w9 A2 `! I; ?. ~8 I2 l5 l5 U2 ^$ T1! k+ R1 D' b: Z8 _3 X7 t; f9 { `
OK,这就是NRI的计算,除此之外,随机森林、决策树、lasso回归、SVM等,这些模型,都是可以计算的NRI的,后面会继续介绍。大家如果有问题欢迎在评论区留言。
2 O2 l3 B1 _! X; C9 W" A6 Y' |9 U, y5 U! d- S5 m/ Y
本文首发于公众号:医学和生信笔记
5 Q( ?) D+ ~2 }# C
% o* F- ?7 ?' o“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。# N9 z: B% f4 l3 V6 [; ^3 l% _$ v
本文由 mdnice 多平台发布
$ R% |$ |, i7 S————————————————- f" u0 c* ?- l+ v- K/ S
版权声明:本文为CSDN博主「阿越就是我」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
3 ^7 ]2 ?1 x: ^' p2 J原文链接:https://blog.csdn.net/Ayue0616/article/details/1267680061 ~6 ~. [( T8 s! W, t* D! v u
9 y' z: a/ }4 @. w; a# U! U/ t/ O) ~ \& R- f
|
zan
|