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

###logistic回歸

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

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

family=binomial('logit'),

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

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

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

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

collapse='+')))

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

data=egData,x=TRUE)

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

data=egData,x=TRUE)

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

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

niter=10000,alpha=0.05)

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

pstd<-mstd$fitted.values

pnew<-mnew$fitted.values

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