風(fēng)險(xiǎn)預(yù)測(cè)模型評(píng)價(jià)之NRI的R語言計(jì)算

風(fēng)險(xiǎn)預(yù)測(cè)模型評(píng)價(jià)之NRI的R語言計(jì)算R語言代碼的盛（暴）宴（擊），除了NRI的運(yùn)算，還有蠻多預(yù)處理的操作，在哪都能用得著。大家做好戰(zhàn)斗準(zhǔn)備。R里有2個(gè)包專事計(jì)算NRI，分別為nricens和PredictABEL。從最后結(jié)果來說，nricens計(jì)算出來的是絕對(duì)NRI，PredictABEL則為相對(duì)NRI。但我們已經(jīng)知道計(jì)算原理了呀，而且它們都能生成新舊模型分類的對(duì)照表，所以其實(shí)只用其中一個(gè)包就都可以計(jì)算了。不過它們還是有一些小小差異，我們就以logistic回歸模型為例，分別看一下這兩個(gè)包，供大家參考選擇。Cox模型參數(shù)較多也較復(fù)雜，但我相信你看完這篇的講解就能看懂幫助文檔中的cox案例，算是留個(gè)小作業(yè)給你吧~擬合模型先用一份示例數(shù)據(jù)做個(gè)模型。這是survival包里帶有的一份梅奧診所的數(shù)據(jù)，記錄了418位患者的臨床指標(biāo)，觀察這些因素與原發(fā)性膽汁性肝硬化（PBC）的關(guān)系。其中前312位參加了RCT，其他的只參加了觀察隊(duì)列。我們用前312份樣本，做個(gè)預(yù)測(cè)2000天時(shí)間點(diǎn)上死亡與否的模型。先載入這份數(shù)據(jù)看一下。library(survival)

###logistic回歸

egData<-pbc[1:312,]做一個(gè)logistic回歸，我們需要一個(gè)結(jié)局事件作為因變量，它必需是個(gè)分類變量；其次需要若干自變量，它們可以是分類也可以是連續(xù)。這個(gè)表中的結(jié)局是status，0=截尾（刪失），1=接受移植，2=死亡。研究目的“死亡與否”是個(gè)二分類變量，所以要做些變換。再看time一欄，有的不夠2000天，這些樣本要么是沒到2000天就死亡了，要么是刪失了。我們要?jiǎng)h掉2000天內(nèi)刪失的數(shù)據(jù)。egData=egData[egData$time>2000|

(egData$time<2000&egData$status==2),]“[]”表示篩選條件，|表示“或”，&為“和”。所以條件句就是egData中的time一列大于2000的保留，或小于2000但同時(shí)狀態(tài)為死亡的也保留。最后一個(gè)“,”別忘了，其在條件句的前面表示對(duì)列進(jìn)行選擇，在其后表示選擇行。選好后做一個(gè)event向量，把status的三個(gè)狀態(tài)變成死亡=1,未死亡=0。event=ifelse(egData$time<2000&egData$status==2,1,0)ifelse(test,yes,no)大法好啊，前面一個(gè)test是邏輯判斷句，其值為真時(shí)返回yes的值，為假時(shí)返回no的值。所以本句中test就是當(dāng)time<2000，且status為2（死亡）時(shí)，記為1，否則為0。然后把event合并入原來的表格。egData=cbind(egData,event)cbind()是以列合并，另有rbind()以行合并。這樣event就成了最后一列，為結(jié)局事件。然后選擇模型的自變量（predictors）。太多了，選取其中幾個(gè)示例。就以年齡、膽紅素、白蛋白為舊模型（standard），三者加上一個(gè)凝血酶原時(shí)間為新模型（new）。一般做logistic回歸是用glm(因變量~自變量1+自變量2+……+自變量n,

family=binomial('logit'),

data=數(shù)據(jù)表)，但如果自變量較多的話，前面那個(gè)運(yùn)算式就會(huì)很長很長，萬一這些自變量還是基因名或編號(hào)，就很想死了。所以順便講一個(gè)簡化的辦法，即把那一串先寫成formula。fml.std<-as.formula(paste('event~',

paste(colnames(egData)[c(5,11,13)],

collapse='+')))這里有好幾層函數(shù)，paste()會(huì)把括號(hào)中的元素粘貼起來，collapse是其中的間隔。colnames()是獲取表格的列名，[]中的數(shù)值向量為所選擇的列序號(hào)。這樣如果是一個(gè)超大表格，你選中第10~70列還可以寫成“10:70”。好了，同樣寫出新模型的formula：fml.new<-as.formula(paste('event~',

paste(colnames(egData)[c(5,11,13,19)],

collapse='+')))

可以查看一下，新模型的formula寫成這個(gè)效果：然后像上邊說的那樣用glm()擬合兩個(gè)模型。mstd=glm(fml.std,family=binomial('logit'),

data=egData,x=TRUE)

mnew=glm(fml.new,family=binomial('logit'),

data=egData,x=TRUE)

這樣一長串運(yùn)算式用剛才命名好的fml.std和fml.new代替就好了。x=TRUE是將來用nricens包計(jì)算時(shí)要求用到的，表示輸出結(jié)果中是否包含所用到的數(shù)據(jù)表，平時(shí)可以不寫。模型就這樣做完了~先不急著計(jì)算NRI，先看看它的總體情況。summary(mstd)運(yùn)行這句就得到該模型的描述特征。殘差、相關(guān)系數(shù)、各個(gè)自變量的統(tǒng)計(jì)顯著性等，注意倒數(shù)第二行的AIC，就是上一期提到的赤池信息準(zhǔn)則，表示模型校準(zhǔn)度，很少有人匯報(bào)呢?？梢杂猛瑯拥姆椒纯葱履Ｐ汀＿@里就-不展開了，進(jìn)入下一環(huán)節(jié)。NRI的計(jì)算?先看nricens包的方法。library(nricens)NRI<-nribin(mdl.std=mstd,mdl.new=mnew,

updown='category',cut=c(0.3,0.6),

niter=10000,alpha=0.05)

填上新舊兩個(gè)模型。Cut是判斷風(fēng)險(xiǎn)高低的臨界值，現(xiàn)在我們寫了2個(gè)，也就是0~29%為低風(fēng)險(xiǎn)，30%~59%為中風(fēng)險(xiǎn)，60%~100%為高風(fēng)險(xiǎn)?，F(xiàn)實(shí)中可以查閱相關(guān)文獻(xiàn)進(jìn)行設(shè)置，預(yù)測(cè)風(fēng)險(xiǎn)達(dá)到多少需要怎樣干預(yù)之類的。Updown為定義一個(gè)樣本的風(fēng)險(xiǎn)是否變動(dòng)的方式，category是指分類值，即我就熟悉的低、中、高風(fēng)險(xiǎn)，另有一種diff，為連續(xù)值。選diff時(shí)，cut就設(shè)1個(gè)值，比如0.02，即認(rèn)為當(dāng)預(yù)測(cè)的風(fēng)險(xiǎn)在新舊模型中相差2%時(shí)，即被認(rèn)為是重新分類了。這種方法用的比較少。后面幾個(gè)參數(shù)就比較有意思了，niter為重復(fù)取樣的次數(shù)，即boostrap方法，不做的話將其設(shè)為0就好了；做的話建議至少1000次，這也是默認(rèn)值，但我（讀書少）見過的研究都10000次。然后將統(tǒng)計(jì)顯著性alpha設(shè)為0.05。這樣就可以看到輸出的結(jié)果：如果不做bootstrap，就是這個(gè)結(jié)果。有重新分類情況的詳表，最后是NRI和各種變動(dòng)的概率。第一個(gè)NRI如前所述，是絕對(duì)NRI，大家可以根據(jù)之前的知識(shí)和上邊的詳表自己計(jì)算驗(yàn)證一下，此時(shí)可手動(dòng)計(jì)算出相對(duì)NRI。其他指標(biāo)隨便看看。如果做了bootstrap，就會(huì)多出一個(gè)表：因?yàn)樽隽?0000次重復(fù)取樣，相當(dāng)于有10000個(gè)NRI，于是就有了標(biāo)準(zhǔn)誤和置信區(qū)間，剛才我們?cè)O(shè)alpha=0.05，所以后面的Lower和Upper就是95%置信區(qū)間的下界和上界。同時(shí)，做不做bootstrap都會(huì)得到一張圖，表示各數(shù)據(jù)點(diǎn)在新舊模型中的分布。默認(rèn)的Case和Control標(biāo)簽我覺得不太嚴(yán)謹(jǐn)，Case代表結(jié)局事件中編號(hào)為“1”的組，也就是發(fā)生了結(jié)局（死亡），Control為“0”，未發(fā)生。其實(shí)是positive和negative比較貼切吧。反正就這個(gè)意思。這張圖也和重新分類表的意思差不多，看看就好。?再看PredictABEL包的做法library(PredictABEL)

pstd<-mstd$fitted.values

pnew<-mnew$fitted.values

先把兩個(gè)模型中的預(yù)測(cè)風(fēng)險(xiǎn)值提出來，也就是模型中的fitted.value。這個(gè)包只能從預(yù)測(cè)風(fēng)險(xiǎn)計(jì)算，剛才的nricens包可以用模型，也可以用預(yù)測(cè)風(fēng)險(xiǎn)（把mdl.std和mdl.new參數(shù)換成p.std和p.new）。reclassification(data=egData,cOutcome=21,

predrisk1=pstd,predrisk2=pnew,

cutoff=c(0,0.30,0.60,1))cOutcome是結(jié)局事件的列序號(hào)，剛才我們不是把event放到最后了么，即第21列，填上。兩個(gè)預(yù)測(cè)風(fēng)險(xiǎn)值也相應(yīng)填上。這里的cutoff跟剛才的不一樣，還要填上前面的0和后面的1，成為完整的0~100%的區(qū)間。然后得到一個(gè)重新分類表：跟上邊nricens做的差不多了。不過這個(gè)包沒有bootstrap的選項(xiàng)。接著看下面的結(jié)果，這里的分類NRI是咱們上回說的相加NRI，同樣可以根據(jù)上一期的知識(shí)手動(dòng)計(jì)算一下。記得咱們并沒有設(shè)置bootstrap吧？可這里也有個(gè)95%置信區(qū)間，只是內(nèi)部調(diào)用了一個(gè)更為簡陋的只能計(jì)算連續(xù)NRI的improveProb()函數(shù)做的，而且連續(xù)變化的臨界值也不太透明，遂不管了。最后還有個(gè)IDI是指，發(fā)生和未發(fā)生結(jié)局事件樣本的平均預(yù)測(cè)