- 在线时间
- 1630 小时
- 最后登录
- 2024-1-29
- 注册时间
- 2017-5-16
- 听众数
- 81
- 收听数
- 1
- 能力
- 120 分
- 体力
- 544437 点
- 威望
- 12 点
- 阅读权限
- 255
- 积分
- 168709
- 相册
- 1
- 日志
- 0
- 记录
- 0
- 帖子
- 5324
- 主题
- 5250
- 精华
- 18
- 分享
- 0
- 好友
- 163
TA的每日心情![](source/plugin/dsu_paulsign/img/emot/kx.gif) | 开心 2021-8-11 17:59 |
---|
签到天数: 17 天 [LV.4]偶尔看看III 网络挑战赛参赛者 网络挑战赛参赛者 - 自我介绍
- 本人女,毕业于内蒙古科技大学,担任文职专业,毕业专业英语。
![](plugin.php?id=eis_qrcode2:make_qrcode&tid=491532) 群组: 2018美赛大象算法课程 群组: 2018美赛护航培训课程 群组: 2019年 数学中国站长建 群组: 2019年数据分析师课程 群组: 2018年大象老师国赛优 |
- o. ?: d; w! G' y0 j6 g& [: J3 ^( e净重新分类指数NRI的计算
5 p7 p5 H$ Q! D8 h“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
+ A: ~# n. _2 kNRI,net reclassification index,净重新分类指数,是用来比较模型准确度的,这个概念有点难理解,但是非常重要,在临床研究中非常常见,是评价模型的一大利器!& h5 a% J8 N! _8 q. Y
5 d" F7 P) V" l$ E在R语言中有很多包可以计算NRI,但是能同时计算logistic回归和cox回归的只有nricens包,PredictABEL可以计算logistic模型的净重分类指数,survNRI可以计算cox模型的净重分类指数。. i; j: @/ ~ k6 m$ M* x+ [ r5 `+ J
( d: [+ I" Y$ C6 ?logistic的NRI. A5 |/ L+ n, U/ Y
nricens包0 S2 A: g7 s( D% K0 {
PredictABEL包
$ ~) ?1 l9 K% N2 e1 Q( E( V9 H生存分析的NRI5 i G7 w$ \5 h E: H
nricens包+ ~2 m' c- D6 F9 r
survNRI包! T% v: e* U# }0 w4 T
logistic的NRI9 }- c% c/ e8 N& U( R* \! ^5 m
nricens包# U. [+ J e0 E
#install.packages("nricens") # 安装R包 z* R# V8 a! S% C
library(nricens)
3 ^2 x+ c+ x* N4 R9 N7 l1) Z* Q; [9 @( J* I! r9 x' m: b4 r
## Loading required package: survival+ v5 ~' ^* m! A: `2 y
1
+ V9 r/ n7 G' ~; u5 ]+ N使用survival包中的pbc数据集用于演示,这是一份关于原发性硬化性胆管炎的数据,其实是一份用于生存分析的数据,是有时间变量的,但是这里我们用于演示logistic回归,只要不使用time这一列就可以了。) A0 C8 u, b$ l8 k$ ^( p/ ~
, q" O+ \4 `' llibrary(survival); V( v7 k) P1 i& d( Z# w% W. H
, }" f3 s8 p: D6 o; Q
# 只使用部分数据
6 u$ X% z) ~0 `/ }1 [dat = pbc[1:312,] & h( U! W( v G8 Z+ n5 f- M# [
dat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]
5 X$ I5 J- ^3 `$ I& p
; Q( r+ x' [6 F l) K8 ?str(dat) # 数据长这样0 i" \, K0 T5 B! H2 n
1
% E5 _4 V2 I( R9 A. k# e## 'data.frame': 232 obs. of 20 variables:
2 R" A1 d9 }' f## $ id : int 1 2 3 4 6 8 9 10 11 12 ...' X! ^" ^) b" i
## $ time : int 400 4500 1012 1925 2503 2466 2400 51 3762 304 ...
3 S0 H* b" V) \## $ status : int 2 0 2 2 2 2 2 2 2 2 ...+ a' X) a+ N. _
## $ trt : int 1 1 1 1 2 2 1 2 2 2 ...0 X+ p7 b0 Y4 I5 B
## $ age : num 58.8 56.4 70.1 54.7 66.3 ...' Y2 x) F& u% w
## $ sex : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ... ~/ e: W! ~" L, m: P0 h
## $ ascites : int 1 0 0 0 0 0 0 1 0 0 ...0 i5 G0 H4 M7 B }0 C) G2 }
## $ hepato : int 1 1 0 1 1 0 0 0 1 0 ...2 d9 u! w! X" d5 F4 V1 V! {+ d* K4 M- c
## $ spiders : int 1 1 0 1 0 0 1 1 1 1 ...
0 b' G" J+ \6 J: Y0 n- l8 y## $ edema : num 1 0 0.5 0.5 0 0 0 1 0 0 ...
{) ]& G3 {' x( K U## $ bili : num 14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...
- i* H4 k/ R) R( d8 n- K% g" K0 v## $ chol : int 261 302 176 244 248 280 562 200 259 236 ...
4 a8 B7 q+ Z6 S3 Q# v `6 x1 R## $ albumin : num 2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...- {' T$ F' ~: Z, p% V- u
## $ copper : int 156 54 210 64 50 52 79 140 46 94 ...
" k2 u' u0 _4 U& t, k2 V# D. @## $ alk.phos: num 1718 7395 516 6122 944 ...$ D# A5 o, u' |8 z9 S% S! C
## $ ast : num 137.9 113.5 96.1 60.6 93 ...4 e3 ?: [7 u- a+ {) T9 ^: q
## $ trig : int 172 88 55 92 63 189 88 143 79 95 ...* U- ~* u$ h. v8 S
## $ platelet: int 190 221 151 183 NA 373 251 302 258 71 .... Z% k, r, f6 ^& U" h+ z! r H8 x
## $ protime : num 12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...
$ t0 f9 G( Z1 K: E/ Q/ H## $ stage : int 4 3 4 4 3 3 2 4 4 4 ...
9 K4 E) L% A, W: E! B/ r9 z
& Y2 S# M! k$ ]5 s0 B1
4 W8 }1 z) s W. Jdim(dat) # 232 20
) Q1 \5 H# r" P* o1
' h: x& O8 l9 M" B1 h## [1] 232 20 T( ^4 M# ^- x K5 f; @
1
9 e/ |+ W0 Q$ S, K6 H然后就是准备计算NRI所需要的各个参数。9 G) n: x8 L( b8 J# Q1 C0 c9 @9 X7 `
6 \8 C6 K* H+ M# 定义结局事件,0是存活,1是死亡5 v) h# R! ^; a. C
event = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)
) Q: |) a" o+ m0 M: R# d/ l/ i* q. ~# U2 R8 S: s' r+ R1 A' _
# 两个只由预测变量组成的矩阵 t" ~- _) m; E) f9 D; ~2 e' Q
z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
" h& }4 X" {6 C- I# s4 iz.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))& q. v# D0 f3 e5 D+ D2 _0 V
1 b- ]( B" w+ k7 \! t& X' E# 建立2个模型3 s$ Y y# e( q0 P6 \
mstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)2 N2 y ]6 @# X9 f
mnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)
+ s& j% E/ m8 ?0 g- c% I$ N
6 @: e' \( I$ K$ @4 p3 t# 取出模型预测概率' [. C4 ?- J( V, J: F3 G" f
p.std = mstd$fitted.values5 B, X* X G: {) I1 s
p.new = mnew$fitted.values
" F4 @ s; H6 t& X3 b$ M7 Q& ~
6 P0 K j$ `1 D h; R1 z" ]1
, B8 U4 `5 ]3 @8 Y1 ^9 y& g3 M6 T然后就是计算NRI,对于二分类变量,使用nribin()函数,这个函数提供了3种参数使用组合,任选一种都可以计算出来(结果一样),以下3组参数任选1组即可。 mdl.std, mdl.new 或者 event, z.std, z.new 或者 event, p.std, p.new。
B5 S, r% f+ W: b; L+ g. Y! p, z1 w" I0 ]
# 这3种方法算出来都是一样的结果
- x; W7 Y T+ ?! z, n9 \0 t
1 ]$ [9 y7 _. @5 f _0 D# 两个模型+ ` d3 ~! u- }9 l/ o# T/ O6 p
nribin(mdl.std = mstd, mdl.new = mnew, $ l( Y2 V! {- {: t2 j$ U& u
cut = c(0.3,0.7), 1 m1 l; L9 N" x( f8 O
niter = 500, / @+ A' G7 {+ L6 |
updown = 'category'); t% q: C- E i9 J% ?
3 u: X+ @/ O! O* ]- Z5 |2 N# 结果变量 + 两个只有预测变量的矩阵" q: a# E" M. f" M
nribin(event = event, z.std = z.std, z.new = z.new, 6 @5 g8 N: u6 e* @, B' `7 q9 ?
cut = c(0.3,0.7), 9 U5 n7 b, s3 q2 h$ t6 f6 z3 a, W
niter = 500, ) ^6 _3 D& o. c( U$ X& Y
updown = 'category')# B0 i1 F3 r3 c# t; S
+ @+ Y! T! y: t2 e3 a; n
## 结果变量 + 两个模型得到的预测概率
! ~( d& `; O0 a# |nribin(event = event, p.std = p.std, p.new = p.new, / ~8 h( I* \; A$ m
cut = c(0.3,0.7), ) X5 X/ t7 ^( `% F
niter = 500,
. k& ^: H" I# E+ I: ^ updown = 'category')
/ {4 ?0 o/ w" c. A8 _ t( t2 H4 N9 X# o; k, E3 i% {
13 l* ^" I {. D$ B) K; n# D
其中,cut是判断风险高低的阈值,我们使用了0.3,0.7,代表0-0.3是低风险,0.3-0.7是中风险,0.7-1是高风险,这个阈值是自己设置的,大家根据经验或者文献设置即可。3 n6 W# y4 u, m2 C/ J( r6 Y" G
: p1 M- h1 [# u) Yniter是使用bootstrap法进行重抽样的次数,默认是1000,大家可以自己设置。9 l/ v; q# h, k
5 {' `! r5 X$ o/ Uupdown参数,当设置为category时,表示低、中、高风险这种方式;当设置为diff时,此时cut的取值只能设置1个,比如设置0.2,即表示当新模型预测的风险和旧模型相差20%时,认为是重新分类。
) f) j2 i0 [. S( _7 i6 B
9 b% \1 r# o/ U+ g; I/ I! c上面的代码运行后结果是这样的:* K: M- M0 L9 T
* r$ p0 \$ M8 y) O& L
UP and DOWN calculation:
- i+ @) ^( i1 O* @! x, t# b/ N" h #of total, case, and control subjects at t0: 232 88 144
: Z! e7 W+ v& X/ O9 j0 L; ~9 q
2 v. p* v) \- {- t! P Reclassification Table for all subjects:* K' A, A4 ` { i ^( O5 {
New3 ^& h' W5 o7 ~ w+ l* `( a
Standard < 0.3 < 0.7 >= 0.7
' u: e+ C& r+ i5 o < 0.3 135 4 0* y: l# w( N/ k# O' u# ]* G
< 0.7 1 31 4
: X' g5 O$ c) V/ [2 J/ i- [, M >= 0.7 0 2 55
! K0 a, e2 J8 J; n- P% V& w, o# C! q3 S! x# m% O
Reclassification Table for case:
7 j7 W8 N* @3 E; G New
7 n+ ?8 Z2 |0 WStandard < 0.3 < 0.7 >= 0.7
) d5 j' x4 n( R0 R! r2 ^ A < 0.3 14 0 0" V7 b' l& |8 N% x- p
< 0.7 0 18 30 H# d; S( e! e2 t( Q# F
>= 0.7 0 1 52" i" W+ Z% y1 A- N0 E1 t
, V( Y" V: r# S- m$ @
Reclassification Table for control:- G; x; u9 C; H4 M/ Z: ?
New
0 e- d2 G6 _; v$ g! ]3 CStandard < 0.3 < 0.7 >= 0.7. s; p0 m* D4 Q& J
< 0.3 121 4 00 h! y3 n; E _0 U3 V& k
< 0.7 1 13 1
; w# ], j- T/ Y6 W3 T6 u, g >= 0.7 0 1 3
) k; w n9 g6 Y; F1 n8 q- z5 G: |7 k$ Y+ _$ W. @ Y" D
NRI estimation:
5 C1 A; y e6 L8 H( IPoint estimates:4 E* V5 P. d$ ]
Estimate
; L9 e/ c% o; i( X% A1 W$ E. ~0 YNRI 0.001893939
* F, S, F; {( n+ O/ hNRI+ 0.022727273/ V8 e, w6 v1 a5 V1 e
NRI- -0.020833333" n9 d& o5 B. U* A2 L5 p5 r
Pr(Up|Case) 0.034090909
: Z) ~- ?2 m, K& \0 z9 o1 p" Q* JPr(Down|Case) 0.011363636; \. W- [- s$ ]5 Y5 w' T
Pr(Down|Ctrl) 0.0138888896 I0 p0 h/ Q/ x: A
Pr(Up|Ctrl) 0.034722222
& v& |$ H* N7 g( v ]6 S! v; X) K- k7 m8 I+ s
Now in bootstrap..% B" @! w1 p) V0 N3 |
0 g7 M( p7 J; H& ^/ @/ r
Point & Interval estimates:6 u4 J' X' C& C" S3 y" M
Estimate Std.Error Lower Upper( a+ ]1 ~: \% U1 D
NRI 0.001893939 0.027816095 -0.053995513 0.0553544495 Q. L) W: R, J5 Y/ x3 J
NRI+ 0.022727273 0.021564394 -0.019801980 0.065789474
- p0 A% M" U2 f' LNRI- -0.020833333 0.017312438 -0.058823529 0.007518797
- ^( v- D: w4 B& U) SPr(Up|Case) 0.034090909 0.019007629 0.000000000 0.072164948! n1 n2 C& {- \+ l
Pr(Down|Case) 0.011363636 0.010924271 0.000000000 0.039603960
. O# l, X5 g/ H% Q" tPr(Down|Ctrl) 0.013888889 0.009334685 0.000000000 0.035211268* k- ^1 x# _5 d* h/ y9 F
Pr(Up|Ctrl) 0.034722222 0.014716046 0.006993007 0.0661764714 ~5 A3 F8 R+ c. r" z
0 {, Q8 X1 N: L% T, I# H/ d1+ J# H0 W8 x/ ^) r: ?2 _) j
首先是3个混淆矩阵,第一个是全体的,第2个是case(结局为1)组的,第3个是control(结局为2)组的,有了这3个矩阵,我们可以自己计算净重分类指数。
! K* G) o7 e7 u) u1 L% V4 W3 l0 {, @* R5 v" x( t' b1 u
看case组:$ [, I: [( x* R$ \1 H1 o
6 [7 y; c; c: k# z1 x7 }
净重分类指数 = ((0+3)-(0+1)) / 88 ≈ 0.022727273: U2 [8 V% [ V
7 W9 V% p. m! R( B* \7 o再看control组:% }$ C' K$ y9 I f! t
# Z. F; L9 J" ]! ]0 U
净重分类指数 = ((1+1)-(4+1)) / 144 ≈ -0.0208333333 ?: O4 O# f1 Q: v* {5 o
3 K& {" w( ?8 B9 b. r: i) n
相加净重分类指数 = case组净重分类指数 + control组净重分类指数 = 2/88 - 3/144 ≈ 0.000315657" |% H0 _5 H B
7 ^: X3 S3 I1 v1 I$ {, H再往下是不做bootstrap时得到的估计值,其中NRI就是绝对净重分类指数,NRI+是case组的净重分类指数,NRI-是control组的净重分类指数(和我们计算的一样哦),最后是做了500次bootstrap后得到的估计值,并且有标准误和可信区间。6 y1 p+ l9 G- H9 ]! H) D/ @; s
7 T: F# F4 V7 b- m% f& J
最后还会得到一张图:
9 ^8 o2 | G9 x, d% M4 e0 V* ]1 |; {$ H2 J
这张图中的虚线对应的坐标,就是我们在cut中设置的阈值,这张图对应的是上面结果中的第一个混淆矩阵,反应的是总体的情况,case是结果为1的组,也就是发生结局的组,control是结果为0的组,也就是未发生结局的组。- P* D* ?$ O! n6 A6 \% N$ X
2 ]( A+ B1 H$ {
P值没有直接给出,但是可以自己计算。# |" `+ Y$ p5 ?6 _0 ]/ Z" ?
/ O! @2 C( K: z, s. e# 计算P值9 k, w( c3 S8 ~, J: k2 ?9 X
z <- abs(0.001893939/0.027816095)
7 F4 o- y5 S- G, g: M3 np <- (1 - pnorm(z))*2
5 E& e1 i4 ?2 d4 l( Bp
4 Z* I9 Y# |- W8 i1
* L$ t- ]1 K) r3 @$ u% g9 V## [1] 0.9457157
& e- v3 t2 e& H; N2 G17 Z, v4 ~5 s: C, E# V9 Q
PredictABEL包0 T5 c4 G, g ~ M8 ?
#install.packages("PredictABEL") #安装R包/ }* _. W* W" G
library(PredictABEL)
0 M/ P7 [% s: \2 O0 J9 a4 W: Q" R. x. P1 @# ]! z
# 取出模型预测概率,这个包只能用预测概率计算
& a: a# M" r: c8 V8 A, w* g: C( op.std = mstd$fitted.values4 `, F0 H7 C I$ u
p.new = mnew$fitted.values
] Z1 w @( [- i; M$ x# I1
, x; C z, E0 I; u2 P- J然后就是计算NRI:( U: v% G& H, T5 ~ a6 Y* W
" V/ y3 m/ `! x" G: s4 }
dat$event <- event& H# v( o+ }9 P" n" i& X* y
$ ?& |8 k% w: d3 j1 A% r/ N
reclassification(data = dat,. V( ^. `+ v/ W) X4 U0 \
cOutcome = 21, # 结果变量在哪一列
0 ]" U/ e1 _ C% |1 o" Q predrisk1 = p.std,3 ~1 s- |0 d) ~$ I6 i/ ?: Q; a
predrisk2 = p.new,* z! e) |) T5 l
cutoff = c(0,0.3,0.7,1)5 Y* G+ q* `' u) t4 k( h3 t
)
5 Y9 b5 p2 t& W! C1
4 `, K, g. @# J$ y: q" o& m## _________________________________________
. h; r/ K( T3 ^& f8 W##
$ d% w/ D k+ `5 i' h## Reclassification table % f: t5 |5 N! J: L
## _________________________________________4 H: y( K3 I' t
##
9 a& M6 J* ~+ u' @9 M## Outcome: absent % T1 K' r' W& c9 H4 y9 Z
##
# P+ e# n7 G K( s' E7 {/ @2 b4 Y## Updated Model
) Y4 S& G+ C8 g! S, C- T5 S## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
, M+ N- H0 Y" {( g* e* g& `% |## [0,0.3) 121 4 0 3
! ?5 X; o* g$ m6 X, f## [0.3,0.7) 1 13 1 13* b2 m/ |- F4 f+ _- E: @
## [0.7,1] 0 1 3 25# F' x) @8 Q1 L% p, n+ C
##
1 v/ Z8 L. |7 m% W- y8 ~! ]8 {## - Z6 F5 o) p# l0 o2 B6 c8 P( C5 x
## Outcome: present
5 \0 v7 q9 u5 I## $ T' \/ j K0 ]$ o& K" Y
## Updated Model
2 `, t4 T) d0 i- y/ D: y2 y0 O' u' W## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
! H" Z2 S% j" G3 t- s" V## [0,0.3) 14 0 0 0
2 K5 l, L! v* W## [0.3,0.7) 0 18 3 14
y H/ l U" g* _$ E9 A## [0.7,1] 0 1 52 2
% ~$ p5 j& C% |##
& {& Q7 `$ V9 K& e# c" }7 V## $ S; {1 ^! o a$ {
## Combined Data
8 _; ?, u- }4 p( \( U2 T##
3 N- F8 |. [( {5 R" E* e0 k" D## Updated Model _7 R, R0 J$ Q+ G: A
## Initial Model [0,0.3) [0.3,0.7) [0.7,1] % reclassified
+ \1 v: q/ b) d, [" ]## [0,0.3) 135 4 0 3' M6 \/ A% A$ |# u2 v3 y) z& o+ ?
## [0.3,0.7) 1 31 4 14
+ `1 ] @4 C9 e/ f7 i## [0.7,1] 0 2 55 4
: l( ^% z3 M! s: A5 z8 |6 X( X## _________________________________________
* @* [8 J* g5 n f## 0 J7 I) b1 T4 @, g) n0 z
## NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806 - R2 U! `- S0 K
## NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048
0 Y% {) z' P; N, }+ e0 |## IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.28396/ ]2 \! U- \! V- ?( P9 \; H/ r; ^0 [
2 e" }, i" V( z E1 N6 J1- k' [+ ?! U: ^+ ?
结果得到的是相加净重分类指数,还给出了IDI和P值。两个包算是各有优劣吧,大家可以自由选择。
: u7 q. h5 N% z$ E, l1 m! N6 z' C. U$ r; D; |' V/ e, i3 V. I
生存分析的NRI
9 ]" _) M( s0 M还是使用survival包中的pbc数据集用于演示,这次要构建cox回归模型,因此我们要使用time这一列了。: ~" c9 y5 z( b7 I
! ]& G/ x; J, Dnricens包
8 z$ K7 v! ^* y1 P, }library(nricens): [/ q& t: X1 j
library(survival)# Q" h1 n' C: m, \
" S# i1 e6 G) S, X& Fdat <- pbc[1:312,], B& [ u6 } e# ~
dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡
I$ o1 h3 ?, W6 x' K' W# R1
. a' {" u x4 M Z, v- }然后准备所需参数:3 }- `& q( Q$ i U9 W. ]. U& e
7 T a# \8 O# s; v0 I% Q# 两个只由预测变量组成的矩阵
7 e# r) i8 z, { K0 ?$ y; hz.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
* h8 l) K3 A M8 R! T2 t2 t& g: xz.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))
1 N, ] R, U2 {% O) a# F+ _- ~0 Y
, k4 _' t/ K, J/ E/ m# 建立2个cox模型. v4 Q9 U$ G* F3 Q& u/ \: ?
mstd <- coxph(Surv(time,status) ~ age + bili + albumin, data = dat, x=TRUE)
# A6 N' k6 f9 S+ Amnew <- coxph(Surv(time,status) ~ age + bili + albumin + protime, data = dat, x=TRUE)
; `2 i1 X6 Q9 N7 ?3 ~
4 _4 f1 @1 \0 {; B8 h3 f0 e# C# 计算在2000天的模型预测概率,这一步不要也行,看你使用哪些参数
9 l" U+ o+ c% X; g. Gp.std <- get.risk.coxph(mstd, t0=2000)& h, N2 o' v7 T9 `- E
p.new <- get.risk.coxph(mnew, t0=2000)
; j( o1 c9 S: ^+ t1 l) C$ f1: z9 g( ]; ]: d" f5 h# ^; J
计算NRI:) C' G f( Z u
$ g1 s+ ~% F% g( hnricens(mdl.std= mstd, mdl.new = mnew, * V) r8 F' g# p1 L
t0 = 2000,
( i" J% u* W% t# {% j+ k& V cut = c(0.3, 0.7),
" e0 q* a& X- T: M niter = 1000,
1 l0 T. E7 f1 ]) u) V! d' }- \ updown = 'category')
9 [9 O) ?+ w- r5 A( ?7 |
6 P! O* `9 ^6 @$ @UP and DOWN calculation:: s6 p/ z7 \+ u5 R/ l: d
#of total, case, and control subjects at t0: 312 88 1441 n/ a6 n0 ~1 j9 y
% f2 q- a N; r& d
Reclassification Table for all subjects:
! x& |; ]3 T/ d) y1 ~5 B$ I% k New
3 {) v6 Q9 N: L) K. k, i1 FStandard < 0.3 < 0.7 >= 0.7
: ?6 G; `4 ~* z! i: e; u < 0.3 202 7 00 u. ^; y# {7 x! ]; Z* c8 _
< 0.7 13 53 6
; f( O! {3 }, D2 \# m6 M >= 0.7 0 0 31
& F2 g' c' y* ~, ?, `2 e% q+ r. H" l+ @: c0 S" M. p
Reclassification Table for case:2 ~8 U, Y: P# \6 `3 `( \6 m
New
9 [& O! m' d/ w; y- R8 eStandard < 0.3 < 0.7 >= 0.7
/ A! e8 X3 n' x' x. A < 0.3 19 3 0
, u, ~' F6 l* ^' B* y < 0.7 3 32 4
3 }: L( W( S5 Q1 S' O >= 0.7 0 0 27: C( w- |$ |) B- U
; \4 B+ b" ~" ? u8 r; {& J& `
Reclassification Table for control:
3 R( B+ j$ ~) D' }1 i( R New4 T0 `8 E' }* \/ T1 {
Standard < 0.3 < 0.7 >= 0.75 z4 Z7 G& `* K( @% A
< 0.3 126 3 0
0 P) u9 n& C9 N0 T < 0.7 5 7 2
) P( b- j3 R4 y7 L >= 0.7 0 0 1! S' o1 F* ]5 L4 Q
2 p, p, M6 E- U. UNRI estimation by KM estimator:5 m5 m& g) ~: z& [
! f; F+ M8 F' {& z
Point estimates:: e$ U3 V I! Z1 L" s& \4 A
Estimate
5 b2 O% Q) _( i3 ^8 A6 K3 m: UNRI 0.05377635" {1 ]5 m7 p8 A% C
NRI+ 0.03748660
. S% T6 s3 x% |7 LNRI- 0.01628974
9 n E" ~& ]( U" r0 _7 Y9 W* qPr(Up|Case) 0.077089388 d" R9 G. \0 F! T& g
Pr(Down|Case) 0.03960278& z$ g5 D/ T4 d( ]8 L1 M
Pr(Down|Ctrl) 0.04256352
( b6 E# N( o# V% O1 G8 \) [Pr(Up|Ctrl) 0.02627378
2 G3 v' @4 J. g6 ~) O' U2 }4 e8 ^% z; C) p3 u+ B0 J. k; B9 g* u2 X- y
Now in bootstrap..* L9 L% h. s4 v, x& B
0 J% ~2 J M5 p# p) F& M. UPoint & Interval estimates:- E0 q- L- S9 |
Estimate Lower Upper
5 l }6 c. C8 m) s; M8 ^% KNRI 0.05377635 -0.082230381 0.16058172
" R3 t" G, o' W# l1 S6 n, x4 BNRI+ 0.03748660 -0.084245197 0.13231776' w6 b# x f7 H2 i, {
NRI- 0.01628974 -0.030861213 0.06753616
2 {. ?+ g2 d. k0 d2 Y8 y* s1 D {Pr(Up|Case) 0.07708938 0.000000000 0.19102291' d: o7 F3 X; q+ ~
Pr(Down|Case) 0.03960278 0.000000000 0.15236016# ~0 i4 ^- o. ?: |8 l
Pr(Down|Ctrl) 0.04256352 0.004671535 0.09863170; u' `9 \* L ^( _
Pr(Up|Ctrl) 0.02627378 0.006400463 0.05998424
: N) R9 Q- t& {; O: T. p6 q* e0 l6 p
, L1 K: F5 ?0 {: x4 I/ A* z1
6 z- f* B( ]9 I8 M
1 i7 S9 C# | k5 }- ~, U: O4 V5 {Snipaste_2022-05-20_21-49-38, S8 Q3 c5 ]2 G- p9 ?
结果的解读和logistic的一模一样。% s8 B7 z4 ^5 F9 b
$ _, z" C) |, ?" D% x
survNRI包: O! W- J9 q, y/ E
# 安装R包
3 C8 ~( b* b# u4 }devtools::install_github("mdbrown/survNRI")
9 |( Z) u O! L: i1
* ^! [" t! V; `7 X I6 V$ g加载R包并使用,还是用上面的pbc数据集。4 a1 k& c8 o% D4 R) w/ h# D8 F
7 j* P& n, N @! O1 F8 p, f* Wlibrary(survNRI)
0 Z' U7 q' D; V1
$ G6 H' {& D5 F' S## Loading required package: MASS
5 m$ [2 W& {+ }1 E9 f1. N) g; v# w9 {& u9 R
library(survival)- P! H$ f8 a q# w, K9 X! C" X
3 @! V: ?" y* E: |; Y: M! B- }# 使用部分数据/ m4 Z! i! N" k" l
dat <- pbc[1:312,]
! F: T; j- K5 ^9 p+ {) Pdat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着,1表示死亡6 h5 k) Y5 `$ I/ N, A8 _6 s- ?
( a" }2 R) ]9 _& _' i, _res <- survNRI(time = "time", event = "status",
4 k* L6 f4 }& s* ? model1 = c("age", "bili", "albumin"), # 模型1的自变量% K/ g/ D" d# m |; @% C7 y4 m
model2 = c("age", "bili", "albumin", "protime"), # 模型2的自变量: Z1 }. L* `8 O; f) o; K; O) e, o2 ~4 s
data = dat, 0 O" g m. s* \! k3 [8 @; A/ Z
predict.time = 2000, # 预测的时间点, s, @' g# w* y/ ^% \2 D
method = "all",
' N a4 Z7 f; O- w4 v* t bootMethod = "normal", 5 e! Z. ~6 A5 C3 Q& R6 b; E/ Q
bootstraps = 500,
' V) e" p. @, e. B% D alpha = .05)
7 m. g% V" r3 u* }( `! d
" x! ]3 ~ O) Q" c- c1
# F1 h/ E. B$ e查看结果,$estimates给出了不同组的NRI以及总的NRI,包括了使用不同方法(KM/IPW/SmoothIPW/SEM/Combined)得到的结果;$CI给出了可信区间。% V; [) M# Q: |* j, D8 n
& m, V. E! e. g+ O1 H
res9 a7 R1 {3 c0 T* }7 W3 H
1
2 d% u7 M; f8 h \, T+ ?) U6 ?' J## $estimates! n2 y; D: H- s* D' k7 I
## NRI.event NRI.nonevent NRI% \' R5 W7 K1 z+ o
## KM 0.20445422 0.3187408 0.52319516 U" f8 T5 ~- ?% a% P; \0 k! h- u
## IPW 0.22424434 0.3273544 0.5515987
% m6 E' @4 n7 `9 u## SmoothIPW 0.19645006 0.3144263 0.5108763
% }- g& n: n/ @; y3 l: ]## SEM 0.07478611 0.2632127 0.3379988+ Y3 N& R2 y' @6 n6 z. Y0 Y
## Combined 0.19633867 0.3143794 0.5107181% G: N! }. \; R6 N. X- V9 x
##
% w. |( T- S9 E4 [, d P! N, ?## $CI# s0 A9 f$ O* L( M
## $CI$NRI.event! ?: l0 h$ B$ C" p7 c
## KM IPW SmoothIPW SEM Combined
4 r4 {. x) O; p# C## lowerbound -0.03915924 -0.02185068 -0.04724202 -0.1162587 -0.0473723% Y8 s1 c( R. }) _- a( E7 B! @
## upperbound 0.44806768 0.47033936 0.44014214 0.2658309 0.4400496
$ e1 V! J( |2 h4 E( X/ b z## ]3 \9 t# o9 m3 O! A F
## $CI$NRI.nonevent
4 E. s9 a2 U! s; W## KM IPW SmoothIPW SEM Combined. V0 h% o$ r/ w% I) A8 c$ \
## lowerbound 0.1317108 0.1396315 0.1286685 0.08638933 0.12864260 W5 K+ c' x n, R5 |
## upperbound 0.7102251 0.7393216 0.6966341 0.51482212 0.6964549
& P1 d) |& N# k; Q1 v## 1 w4 d8 c7 P: A" t4 r4 {
## $CI$NRI
5 s O; b: k* ?$ z6 C## KM IPW SmoothIPW SEM Combined
4 \3 N" |, g7 N" |0 e0 w% a9 f## lowerbound -0.05112533 -0.04569046 -0.05439863 -0.04132364 -0.05443409$ _* |- [1 z* I
## upperbound 0.89306122 0.92464359 0.87970125 0.64253510 0.87953153
! M& W" |+ _* E! i- q* z) N0 U5 C## * U7 P Z& C d. Y- `/ `& ]
##
2 b ^; E0 K* ?7 i. N" i## $bootMethod+ J6 F* W$ { a
## [1] "normal"3 W% U5 f9 q2 l) k) k
##
, ]0 I0 X7 m0 {## $predict.time: {) Y0 Q0 C- }0 `; f6 ]9 {" E
## [1] 2000
- [; u# E5 y0 u+ I+ Q## ( w& l: O, M K' c9 s! D8 s6 E) T3 n
## $alpha
; N5 S+ f, o' @9 C) ^5 z## [1] 0.05
7 k/ e/ {6 z' t; f5 p: `+ k## 9 C E( v! p8 M4 R% c9 B# L6 f$ S% O' ]
## attr(,"class")* z/ z1 r2 }$ ]! d
## [1] "survNRI"; P: t6 z. z* ]" W* I8 ?
; Y, \/ }. n( V- R' x6 {1$ s5 r6 E5 D0 A. n3 V4 p+ l
OK,这就是NRI的计算,除此之外,随机森林、决策树、lasso回归、SVM等,这些模型,都是可以计算的NRI的,后面会继续介绍。大家如果有问题欢迎在评论区留言。
( }8 i3 p% q, k0 q" P
9 c/ Y% p; j3 k5 y) q# i. i本文首发于公众号:医学和生信笔记
2 [% B- W# g3 V. y8 i8 [& @4 x. E9 [
“ 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。# c4 `2 j* O- M$ l: d% P
本文由 mdnice 多平台发布
t7 W- h' w, t1 P1 e, j; O————————————————5 N/ `$ D; x! Z" ?/ B
版权声明:本文为CSDN博主「阿越就是我」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
; m; f& e p! x; g; d0 I原文链接:https://blog.csdn.net/Ayue0616/article/details/126768006
: e& G$ {: w9 v
( i" ~* m1 ~( I/ W
0 D' N9 y0 h$ n D; } |
zan
|