本內(nèi)容為【科研私家菜】R語言機(jī)器學(xué)習(xí)與臨床預(yù)測(cè)模型系列課程

R小鹽準(zhǔn)備介紹R語言機(jī)器學(xué)習(xí)與預(yù)測(cè)模型的學(xué)習(xí)筆記

你想要的R語言學(xué)習(xí)資料都在這里， 快來收藏關(guān)注【科研私家菜】

01 什么是綜合判別改善指數(shù)？

綜合判別改善指數(shù)(Integrated Discrimination Improvement, IDI)。這個(gè)指標(biāo)也用于判斷預(yù)測(cè)模型改善情況，與上一講的凈重新分類指數(shù)(NRI)有類似也有不同。
IDI是由Pencina等人于2008年提出的一個(gè)非常新的判別指標(biāo)。由于它考慮了不同切點(diǎn)的情況，可以用來反映模型的整體改善狀況，在一定程度上補(bǔ)齊了NRI的短板。同時(shí)，雖然AUC也考慮到了不同切點(diǎn)，但是AUC的改善情況在臨床中不易解釋，IDI也因此彌補(bǔ)了AUC的缺陷，可以形象地展示研究對(duì)象被準(zhǔn)確重新判別的比例。
因此，當(dāng)我們?cè)谶M(jìn)行2個(gè)疾病模型比較，或者2個(gè)指標(biāo)診斷效能比較時(shí)，除了傳統(tǒng)的ROC曲線及其AUC，也可以同時(shí)給出NRI和IDI，更加全面多層次的展示模型的改善情況。

02 IDI計(jì)算公式

IDI的計(jì)算其實(shí)也比較簡單，它反映的是兩個(gè)模型預(yù)測(cè)概率差距上的變化，因此是基于疾病模型對(duì)每個(gè)個(gè)體的預(yù)測(cè)概率計(jì)算所得。它的計(jì)算方法為：

其中Pnew,events、Pold,events表示在患者組中，新模型和舊模型對(duì)于每個(gè)個(gè)體預(yù)測(cè)疾病發(fā)生概率的平均值，兩者相減表示預(yù)測(cè)概率提高的變化量，對(duì)于患者來說，預(yù)測(cè)患病的概率越高，模型越準(zhǔn)確，因此差值越大則提示新模型越好。

在模型比較時(shí)，將兩部分相減即可得到IDI，總體來說IDI越大，則提示新模型預(yù)測(cè)能力越好。與NRI類似，若IDI＞0，則為正改善，說明新模型比舊模型的預(yù)測(cè)能力有所改善，若IDI＜0，則為負(fù)改善，新模型預(yù)測(cè)能力下降，若IDI=0，則認(rèn)為新模型沒有改善。

我們可以通過計(jì)算Z統(tǒng)計(jì)量，來判斷IDI與0相比是否具有統(tǒng)計(jì)學(xué)顯著性，統(tǒng)計(jì)量Z近似服從正態(tài)分布，公式如下：

02 二分類變量R語言實(shí)現(xiàn)

示例數(shù)據(jù)來自于survival包里自帶的一份梅奧診所的數(shù)據(jù)，記錄了418位患者的臨床指標(biāo)與原發(fā)性膽汁性肝硬化（PBC）的關(guān)系。其中前312個(gè)患者來自于RCT研究，其他患者來自于隊(duì)列研究。我們用前312例患者的數(shù)據(jù)來預(yù)測(cè)2000天時(shí)間點(diǎn)上是否發(fā)生死亡。此處需要說明的是原始數(shù)據(jù)是一個(gè)生存數(shù)據(jù)，我們重新定義一個(gè)二分類的結(jié)局，死亡or 存活，不考慮時(shí)間因素。先載入這份數(shù)據(jù)，如圖4.所示。這個(gè)表中的結(jié)局變量是status，0 = 截尾（刪失），1 = 接受肝移植，2 = 死亡。膽我們的研究目的“死亡與否”是個(gè)二分類變量，所以要做變量變換。再看time一欄，有的不夠2000天，這些樣本要么是沒到2000天就死亡了，要么是刪失了。我們要?jiǎng)h掉2000天內(nèi)刪失的數(shù)據(jù)。

##hereconsider pbc dataset in survival package as an example

library(survival)

dat=pbc[1:312,]

dat$sex=ifelse(dat$sex==''f'',1,0)

##subjectscensored before 2000 days are excluded

dat=dat[dat$time>2000|(dat$time<2000&dat$status==2),]

##predcitingthe event of ''death'' before 2000 days

event=ifelse(dat$time<2000&dat$status==2,1,0)

##standardprediction model: age, bilirubin, and albumin

z.std=as.matrix(subset(dat,select=c(age,bili,albumin)))

##newprediction model: age, bilirubin, albumin, and protime

z.new=as.matrix(subset(dat,select=c(age,bili,albumin,protime)))

##glmfit (logistic model)

mstd=glm(event~.,binomial(logit),data.frame(event,z.std),x=TRUE)

mnew=glm(event~.,binomial(logit),data.frame(event,z.new),x=TRUE)

##UsingPredictABEL package

library(PredictABEL)

pstd<-mstd$fitted.values

pnew<-mnew$fitted.values

##用cbind函數(shù)把前面定義的event變量加入數(shù)據(jù)集，并定義為dat_new

dat_new=cbind(dat,event)

##計(jì)算NRI,同時(shí)報(bào)告了IDI,IDI計(jì)算與cutoff點(diǎn)設(shè)置無關(guān)。

##cOutcome指定結(jié)局變量的列序號(hào)

##predrisk1,predrisk2為新舊logistic回歸模型

reclassification(data=dat_new,cOutcome=21,

                 predrisk1=pstd,predrisk2=pnew,

                 cutoff=c(0,0.2,0.4,1))
# 計(jì)算的IDI為0.44%，說明新模型較舊模型預(yù)測(cè)能力僅改善0.44%。

03 生存資料模型

生存資料的NRI與分類結(jié)局的NRI差別在于前者需要構(gòu)建Cox回歸模型，所以我們首先構(gòu)建新舊Cox回歸模型，計(jì)算這兩個(gè)模型的NRI.

##here consider pbc dataset in survival package as an example

library(survival)

dat=pbc[1:312,]

dat$time=as.numeric(dat$time)

##定義生存結(jié)局

dat$status=ifelse(dat$status==2,1, 0)

##定義時(shí)間點(diǎn)

t0=365*5

##基礎(chǔ)回歸模型

indata0=as.matrix(subset(dat,select=c(time,status,age,bili,albumin)))

##增加1個(gè)預(yù)測(cè)變量新模型

indata1=as.matrix(subset(dat,select=c(time,status,age,bili,albumin,protime)))

##舊模型中預(yù)測(cè)變量矩陣

covs0<-as.matrix(indata0[,c(-1,-2)])

關(guān)注R小鹽，關(guān)注科研私家菜(VX_GZH: SciPrivate)，有問題請(qǐng)聯(lián)系R小鹽。讓我們一起來學(xué)習(xí) R語言機(jī)器學(xué)習(xí)與臨床預(yù)測(cè)模型