临床预测模型之综合判别改善指数IDI计算

---

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。

---

IDI，综合判别改善指数，也适用于评价不同模型优劣的，比起NRI，IDI能够从整体角度对模型进行评价，和NRI一起使用效果更佳！

logistic模型的IDI

二分类变量的IDI计算使用PredictABEL包。

使用survival包中的pbc数据集用于演示，这是一份关于原发性硬化性胆管炎的数据，其实是一份用于生存分析的数据，是有时间变量的，但是这里我们用于演示logistic回归，只要不使用time这一列就可以了。

library(survival)

# 只使用部分数据
dat = pbc[1:312,] 
dat = dat[ dat$time > 2000 | (dat$time < 2000 & dat$status == 2), ]

str(dat) # 数据长这样

## 'data.frame':	232 obs. of  20 variables:
##  $ id      : int  1 2 3 4 6 8 9 10 11 12 ...
##  $ time    : int  400 4500 1012 1925 2503 2466 2400 51 3762 304 ...
##  $ status  : int  2 0 2 2 2 2 2 2 2 2 ...
##  $ trt     : int  1 1 1 1 2 2 1 2 2 2 ...
##  $ age     : num  58.8 56.4 70.1 54.7 66.3 ...
##  $ sex     : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
##  $ ascites : int  1 0 0 0 0 0 0 1 0 0 ...
##  $ hepato  : int  1 1 0 1 1 0 0 0 1 0 ...
##  $ spiders : int  1 1 0 1 0 0 1 1 1 1 ...
##  $ edema   : num  1 0 0.5 0.5 0 0 0 1 0 0 ...
##  $ bili    : num  14.5 1.1 1.4 1.8 0.8 0.3 3.2 12.6 1.4 3.6 ...
##  $ chol    : int  261 302 176 244 248 280 562 200 259 236 ...
##  $ albumin : num  2.6 4.14 3.48 2.54 3.98 4 3.08 2.74 4.16 3.52 ...
##  $ copper  : int  156 54 210 64 50 52 79 140 46 94 ...
##  $ alk.phos: num  1718 7395 516 6122 944 ...
##  $ ast     : num  137.9 113.5 96.1 60.6 93 ...
##  $ trig    : int  172 88 55 92 63 189 88 143 79 95 ...
##  $ platelet: int  190 221 151 183 NA 373 251 302 258 71 ...
##  $ protime : num  12.2 10.6 12 10.3 11 11 11 11.5 12 13.6 ...
##  $ stage   : int  4 3 4 4 3 3 2 4 4 4 ...

dim(dat) # 232 20

## [1] 232  20

然后就是准备计算IDI所需要的各个参数。

# 定义结局事件，0是存活，1是死亡
event = ifelse(dat$time < 2000 & dat$status == 2, 1, 0)

# 建立2个模型
mstd = glm(event ~ age + bili + albumin, family = binomial(), data = dat, x=TRUE)
mnew = glm(event ~ age + bili + albumin + protime, family = binomial(), data = dat, x=TRUE)

# 取出模型预测概率
p.std = mstd$fitted.values
p.new = mnew$fitted.values

接下来就是使用PredictABEL计算IDI：

#install.packages("PredictABEL") #安装R包
library(PredictABEL)  

dat$event <- event

reclassification(data = dat,
                 cOutcome = 21, # 结果变量在哪一列
                 predrisk1 = p.std,
                 predrisk2 = p.new,
                 cutoff = c(0,0.3,0.7,1)
                 )

##  _________________________________________
##  
##      Reclassification table    
##  _________________________________________
## 
##  Outcome: absent 
##   
##              Updated Model
## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
##     [0,0.3)       121         4       0               3
##     [0.3,0.7)       1        13       1              13
##     [0.7,1]         0         1       3              25
## 
##  
##  Outcome: present 
##   
##              Updated Model
## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
##     [0,0.3)        14         0       0               0
##     [0.3,0.7)       0        18       3              14
##     [0.7,1]         0         1      52               2
## 
##  
##  Combined Data 
##   
##              Updated Model
## Initial Model [0,0.3) [0.3,0.7) [0.7,1]  % reclassified
##     [0,0.3)       135         4       0               3
##     [0.3,0.7)       1        31       4              14
##     [0.7,1]         0         2      55               4
##  _________________________________________
## 
##  NRI(Categorical) [95% CI]: 0.0019 [ -0.0551 - 0.0589 ] ; p-value: 0.94806 
##  NRI(Continuous) [95% CI]: 0.0391 [ -0.2238 - 0.3021 ] ; p-value: 0.77048 
##  IDI [95% CI]: 0.0044 [ -0.0037 - 0.0126 ] ; p-value: 0.28396

IDI在最后一行，同时给出了95%的可信区间和P值；还给出了NRI和P值。

生存资料的IDI

生存资料的IDI计算使用survIDINRI包计算。

# 安装R包
install.packages("survIDINRI")

加载R包并使用，还是用上面的pbc数据集。

library(survIDINRI)
## Loading required package: survC1

# 使用部分数据
dat <- pbc[1:312,]
dat$status <- ifelse(dat$status==2, 1, 0) # 0表示活着，1表示死亡

str(dat)

## 'data.frame':	312 obs. of  20 variables:
##  $ id      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ time    : int  400 4500 1012 1925 1504 2503 1832 2466 2400 51 ...
##  $ status  : num  1 0 1 1 0 1 0 1 1 1 ...
##  $ trt     : int  1 1 1 1 2 2 2 2 1 2 ...
##  $ age     : num  58.8 56.4 70.1 54.7 38.1 ...
##  $ sex     : Factor w/ 2 levels "m","f": 2 2 1 2 2 2 2 2 2 2 ...
##  $ ascites : int  1 0 0 0 0 0 0 0 0 1 ...
##  $ hepato  : int  1 1 0 1 1 1 1 0 0 0 ...
##  $ spiders : int  1 1 0 1 1 0 0 0 1 1 ...
##  $ edema   : num  1 0 0.5 0.5 0 0 0 0 0 1 ...
##  $ bili    : num  14.5 1.1 1.4 1.8 3.4 0.8 1 0.3 3.2 12.6 ...
##  $ chol    : int  261 302 176 244 279 248 322 280 562 200 ...
##  $ albumin : num  2.6 4.14 3.48 2.54 3.53 3.98 4.09 4 3.08 2.74 ...
##  $ copper  : int  156 54 210 64 143 50 52 52 79 140 ...
##  $ alk.phos: num  1718 7395 516 6122 671 ...
##  $ ast     : num  137.9 113.5 96.1 60.6 113.2 ...
##  $ trig    : int  172 88 55 92 72 63 213 189 88 143 ...
##  $ platelet: int  190 221 151 183 136 NA 204 373 251 302 ...
##  $ protime : num  12.2 10.6 12 10.3 10.9 11 9.7 11 11 11.5 ...
##  $ stage   : int  4 3 4 4 3 3 3 3 2 4 ...

构建参数需要的值：

# 两个只由预测变量组成的矩阵
z.std = as.matrix(subset(dat, select = c(age, bili, albumin)))
z.new = as.matrix(subset(dat, select = c(age, bili, albumin, protime)))

然后使用IDI.INF()函数计算IDI：

res <- IDI.INF(indata = dat[,c(2,3)],
               covs0 = z.std,
               covs1 = z.new,
               t0 = 2000, # 时间点
               npert = 500 # 重抽样次数
               )

IDI.INF.OUT(res) # 提取结果

##     Est.  Lower Upper p-value
## M1 0.020 -0.003 0.058   0.080
## M2 0.202 -0.042 0.384   0.064
## M3 0.011  0.000 0.036   0.052

m1：IDI的值，可信区间，P值

m2：NRI的值，可信区间，P值

m3：Median improvement in risk score，可信区间，p值。

以上就是IDI的计算方法。

除此之外，随机森林、决策树、lasso回归等也是可以计算IDI的，后面会继续介绍。

---

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。

---