兩種比較模型預測效果的評價指標及其R實現*

李彤陽，林卓琛，葛琪，陳雯，張晉昕

510080 廣州，中山大學 公共衛(wèi)生學院

在醫(yī)學研究特別是臨床醫(yī)學研究中，Logistic回歸是常見的分析方法[1-2]。以兩分類logistic回歸為例，模型的因變量為患病與否，自變量為一系列與該病相關的指標。通過建立Logistic回歸模型我們可以通過每個觀察對象的疾病相關指標，計算出其患病預測概率。目前，Logistic回歸在臨床腫瘤患病風險研究中應用廣泛：如基因多態(tài)性與食管癌患病風險關系研究[3]，2型糖尿病并發(fā)惡性腫瘤風險研究[4]，以及卵巢腫瘤良惡性的鑒別診斷[5]，也有癌癥術后不良反應發(fā)生率及其相關危險因素的研究等[6-7]。

在實際研究中，如果涉及到兩個Logistic模型預測效果的比較，例如兩種診斷試驗的比較，或是在常規(guī)影響因素中添加某一新的生物標志物或感興趣的危險因素，我們常比較兩個模型的靈敏度、特異度、準確度、受試者工作特征(receiver operating characteristic，ROC)曲線[8]以及曲線下面積(area under curve，AUC)[9-10]。通常，Logistic回歸會給出模型的分類表，據此可以看到兩個模型各自的靈敏度、特異度、準確度等指標，若是想進一步整體比較兩模型孰優(yōu)孰劣，則需比較AUC。 AUC綜合了ROC曲線上所有點作為診斷截點的情況，直觀地反映試驗的整體預測效果，在臨床上尤其是診斷試驗的比較中應用廣泛。但在實際應用中，我們常常更關注某一診斷截點下兩個模型對個體的分類能力，這時AUC就不那么好用了；此外，AUC的敏感性有時并不強，其增量并不明顯，專業(yè)意義也不易理解[11-12]。如果我們關心新模型相比于舊模型，對于個體的預測概率改善了多少，或是判斷新加入的生物標志物是否影響到臨床決策，那么只采用常規(guī)的方法就不夠了。近年來，其他比較新舊模型預測效果的方法相繼提出，其中運用最多的是Pencina提出的凈重新分類指數(net reclassification improvement，NRI)和整體鑒別指數(integrated discrimination improvement，IDI)[11,13-15]。NRI關注的是在診斷截點處，通過考察使用新模型后個體預測概率的變化情況，或個體被重新分類的情況，得出新模型比舊模型使得不同組別個體更有利于分到正確組的概率，并對兩組預測概率的改善求和而IDI則是使用兩個模型預測概率均值的差來整體評價新模型改善的程度，故NRI和IDI更貼合應用需求[11,13-14]。

1 NRI和IDI的意義與計算

1.1 常用的模型評價指標

在Logistic回歸模型中，由自變量計算預測概率，當預測概率大于某截點值時，可認為該觀察對象屬于1組，否則屬于0組。而根據課題研究背景，常以1代表患病(y=1)，0代表不患病(y=0)，此時我們就可以認為預測概率高于該截點時，該觀察對象被預測為患病組(反之被預測為不患病組)。

通過診斷試驗四格表可以計算出靈敏度、特異度、約登指數等一系列指標。靈敏度為患病組中個體被正確預測為患者的概率，特異度為不患病組中個體被正確預測為健康者的概率，約登指數為靈敏度與特異度之和減去1。

另外，根據連續(xù)型變量(如預測概率)截點大小的選取，得到的靈敏度和特異度也會相應地改變，這種變化通過連線的形式呈現出來就是ROC曲線，計算ROC曲線下面積則可以得到AUC值，它是對連續(xù)型變量無固定截點時預測能力的概括性刻畫。

1.2 NRI指數

凈重新分類指數NRI[11,13-15]，評價的是在兩模型采用最優(yōu)診斷截點進行預測時，與舊模型相比，使用新模型使得個體的預測結果得到改善的概率，包括事件發(fā)生組改善概率與不發(fā)生組改善概率兩部分(在R中分別以NRI+和NRI-表示)。R中對于NRI指數的計算有兩種不同的方式，一是將預測風險概率視為連續(xù)變量，對應于R中的diff方法；一是將風險概率劃分為若干區(qū)間，對應于R中的category方法。

1.2.1 diff法NRI指數 diff法NRI指數以某一概率變化值(概率變化的大小，而非概率截點)作為判斷新舊模型預測概率變化的依據，只要個體的預測概率變化超過了某一變化值，例如0.02，則認為該個體在新模型下預測的概率和舊模型的不同。相應的的計算公式如下：

NRI=P(up|event)-P(down|event)+
P(down|nonevent)-P(up|nonevent)

(1)

為了便于理解這個公式，我們以預測結局為患病、不患病兩種情況的研究為例，其中up|event和down|event分別表示患病組中個體采用新模型后預測概率值上升和下降的情形(概率變化的幅度超過了設定值才認為是上升或下降)，down|nonevent和up|nonevent表示不患病組中，個體采用新模型后預測的概率值下降和上升的情形(上升和下降的定義同前)，并分別對應于各情形發(fā)生的概率。在實際計算中，可以使用以下公式計算NRI指數：

(2)

其中#up_event和 #down_event 分別表示患病組中新模型預測的患病概率上升和下降的人數(上升和下降的定義同前)，相應的，#up_nonevent和#down_nonevent分別表示不患病組中新模型預測的患病概率上升和下降的人數(上升和下降的定義同前)。#event表示患病組人數，#nonevent表示不患病組人數。

由上述公式可知，NRI指數越大，則與舊模型相比，在新模型中個體無論是本來屬于患病組還是不患病組，均有更大機會被分至實際所屬的組別，預測結果的正確率得以提高，使用新模型后患病組計算出來的預測概率(P(y=1))更高了，而不患病組計算出來的預測概率(P(y=1))更低了。

1.2.2 Category法NRI指數 Category法NRI指數是將預測風險概率以若干個概率截點劃分為感興趣的概率區(qū)間，并基于一個重新分類表，考察個體被新模型預測后重新分類到不同概率區(qū)間的情況，它不關注個體概率變化了多少，而是直接分析個體發(fā)生某事件的預測概率所屬的區(qū)間類別得到優(yōu)化的概率。仍以預測結局為患病和不患病的研究為例，當僅設置一個概率截點時，由于實際上是將概率區(qū)間二分類了，其重新分類表相當于一個患病組和不患病組分開統(tǒng)計的配對四格表(見表1)，可按公式(3)計算NRI指數：

(3)

表1 二分類category法NRI指數新舊模型預測結果的重新分類表

Table 1. Reclassification Table for Category-Based NRI with One Cut-off

Old modelNew modelPredicted probability

*Cut: Probability cut-off.

結合公式(3)和表1，#event1和#nonevent0分別表示患病組和不患病組中新模型預測正確而舊模型預測錯誤的人數(b1和c2)，#event0和#nonevent1分別表示患病組和不患病組中新模型預測錯誤而舊模型預測正確的人數(c1和b2)。#event表示患病組人數，#nonevent表示不患病組人數?？梢姡捎赼1，d1，a2，d2是兩種模型預測的概率區(qū)間相同的情況，在計算NRI指數時并未用到。

同時不難發(fā)現，此時category法的NRI指數等于新模型靈敏度與特異度之和減去舊模型靈敏度與特異度之和，在數值上也相當于新模型的約登指數減去舊模型的約登指數。

當我們設置多個概率截點時，模型預測的風險概率就被劃分成不同的概率區(qū)間(規(guī)定新舊模型的預測概率區(qū)間采用相同的劃分方式)，此時我們認為患病組和不患病組個體患病風險概率的分布不同，理論上患病組的人更多地集中在風險概率高的區(qū)間，不患病組的人更多地集中在風險概率低的區(qū)間，因為模型計算出來的預測概率越大，越有利于該個體預測為患病，同理，模型預測的概率越小，越有利于個體被預測為不患病。此時NRI指數計算公式的形式與公式(2)相同，含義略有差別：

公式中#up_event和#down_event分別表示患者采用新模型預測后所屬的概率區(qū)間比舊模型預測的概率區(qū)間等級升高和降低的人數， #up_nonevent 和#down_nonevent分別表示不患病者采用新模型預測后所屬的概率區(qū)間比舊模型預測的概率區(qū)間等級升高和降低的人數，相應的重新分類表(以設定2個概率截點為例)見表2，其中x1，x5，x9；y1，y5，y9表示預測結果相同的情形。實際上，公式(2)用在此處，僅僅是升高和降低的含義發(fā)生了改變，在diff法中升高和降低的前提是個體預測概率的差值超過了設定的概率變化值，在category法中升高和降低的判斷是依據個體所屬概率區(qū)間等級的改變。另外，category法中只設定一個概率截點的情形也適用于公式(2)，相應含義與設定多個概率截點時相同。

表2 三分類category法NRI指數新舊模型預測結果的重新分類表

Table 2. Reclassification Table for Category-Based NRI with Two Cut-offs

Old modelNew modelPredicted probabili-ty

*Cut: Probability cut-off.

1.2.3 NRI指數的假設檢驗 對于Category法NRI指數，可以使用公式(4)計算其對應的z值，實施總體NRI指數是否等于0的假設檢驗[13]，其中#up_event、#down_event、#up_nonevent 和 #down_nonevent取category法使用公式(2)時對應的含義。

(4)

而對于diff法的NRI指數，可使用Bootstrap(自助)法算出其對應的置信區(qū)間和P值。

1.3 IDI指數

1.3.1 IDI指數的定義與計算 整體鑒別指數IDI[11,13-14]，或稱綜合判別改善指數，它考慮了不同診斷截點的情況，并且從平均概率的角度反映了模型的整體改善狀況。計算公式如下：

(5)

由公式(5)可以發(fā)現，基于預測概率越高越可能為患病組的設定， IDI指數越大，說明新模型可使兩組個體被判定為本應所屬類別的指向性更明確，鑒別能力得以提高。

2 實例分析

本文采用《SPSS 8.0 for Windows統(tǒng)計軟件使用指南》[16]中Logistic回歸分析章節(jié)50例急性淋巴細胞白血病病人資料(文中例14.3)。數據共包含5個變量，其含義分別為：X1：入院治療時取外周血中白細胞數(千個/mm3)；X2：浸潤淋巴結分級(根據浸潤程度分0，1，2 三級)；X3：治療緩解后病人有無鞏固治療(有=1，無=0)；t：隨訪病人生存時間，按月計算；研究結局Y：病人生存時間是否小于1年(一年或以上=0，一年以下=1)。本次研究目的是建立合適的Logistic回歸模型，對淋巴細胞白血病一年內是否死亡進行預測?，F分別建立兩個模型，其中模型1納入自變量入院治療時取外周血中白細胞數和浸潤淋巴結分級，模型2在模型1自變量的基礎上再加上治療緩解后病人有無鞏固治療，模型建立后使用NRI與IDI對兩模型預測效果進行評價。

在R中可使用nricens包和PredictABEL包進行指標的計算，相應的R代碼如下：

#讀入需要使用的包

#install.packages("PredictABEL")

# install.packages("nricens")

library(PredictABEL)

library(nricens)

#建立計算最佳截點函數

opt.cut=function(perf, pred){

cut.ind=mapply(FUN=function(x, y, p){

d=(x-0)^2+(y-1)^2

ind=which(d==min(d))

c(sensitivity=y[[ind]],specificity=1-x[[ind]],cutoff=p[[ind]])

}, perf@x.values, perf@y.values, pred@cutoffs)

}

#讀取數據

setwd("D:data")

dat<-read.csv("data.csv",na.strings="#NULL!")

#擬合模型

mod1<-glm(y～X1+X2,family=binomial(link=′logit′),data=dat)

mod2<-glm(y～X1+X2+X3,family=binomial(link='logit'),data=dat)

pre1<-predict(mod1,newdata=dat,type='response')

pre2<-predict(mod2,newdata=dat,type='response')

pred1<-prediction(pre1,daty);pref1<-performance(pred1,"tpr","fpr")

pred2<-prediction(pre2,daty);pref2<-performance(pred2,"tpr","fpr")

cutpoint1<-opt.cut(pref1,pred1)[3];cutpoint2<-opt.cut(pref2,pred2)[3]

pre_case1<-pre1>cutpoint1;pre_case2=pre2>cutpoint1

table(daty,pre_case1);table(daty,pre_case2)

#使用nricens包計算NRI指數

updown="category",cut=cutpoint2,niter=10000,alpha=0.05)

updown="diff",cut=0,niter=10000,alpha=0.05)

#使用PredictABEL包計算NRI指數和IDI指數

reclassification(data=dat,cOutcome=5,predrisk1=pre1,

predrisk2=pre2,cutoff=c(0,cutpoint2,1))

結果如表3和表4：

表3 基于是否考慮“有無鞏固治療”所得兩模型預測概率的重新分類表

Table 3. Predictive Probabilities of Two Models With or Without the Consolidation Therapy Variable

Model 1Model 2Predictive probability <0.825≥0.825Death within 1year (n=30) <0.825821 ≥0.82501Survival >1 year (n=20) <0.825182 ≥0.82500

本例中選取的概率截點0.825，是通過擬合模型預測概率的ROC曲線，并取約登指數最大值對應的概率確定的。如果已經確定好概率截點，可直接給出相應的賦值。表3相當于二分類category法NRI指數的重新分類表。讀者也可據表3按公式(3)計算NRI(category)，結果與表4中軟件計算的值相同。此外，結果報告時還應給出一年內死亡與未死亡組各自的概率的改善[15]，這與R輸出組與NRI+和NRI-對應(此處以略)。

由表4可知，兩個R軟件包計算的NRI指數結果一致， nricens包沒有提供IDI的計算結果；實際上兩個軟件包都會給出相應的置信區(qū)間，此處受篇幅限制未列出。由表4可知，模型2(包含“有無鞏固治療”的模型)與模型1(不包含“有無鞏固治療”的模型)相比，得到的NRI指數和IDI指數均為正數，即模型2的分類能力優(yōu)于模型1，也即，包含“有無鞏固治療”的模型比不包含的模型分類能力好。具體好在何處？以category法NRI指數為例，值為0.600，說明使用模型2相對于模型1，使得一年內死亡和未死亡的患者更有利于得到正確分類的概率之和為0.600(注意，此處不可理解或表示為百分比或概率[15]),P<0.001，結果具有統(tǒng)計學意義；類似地，IDI指數是0.298，說明采用模型2相比于模型1，在一年內死亡組(Y=1)中預測患者一年內死亡的概率均值的提高量和在一年內未死亡組(Y=0)預測患者一年內死亡的概率均值的降低量之和為0.298，說明模型2預測能力比模型1要好，P<0.001，可見結果具有統(tǒng)計學意義。

表4 基于是否考慮“有無鞏固治療”所得兩個統(tǒng)計模型的比較

Table 4. Two Statistical Models With or Without the Consolidation Therapy Variable

Software packageNRI (category)PNRI (diff)PIDIPnricens0.6002.00×10-81.1671.95×10-07———*———*PredictABEL0.6001.11×10-81.1671.39×10-070.2982.93×10-6

NRI: Net reclassification improvemen; IDI: Integrated discrimination improvement.

* No corresponding results provided.

3 討 論

生物醫(yī)學數據中經常用到Logistic回歸模型來對疾病發(fā)生進行預測。除常見的靈敏度、特異度、約登指數、AUC等指標外，NRI指數和IDI指數也可用于評價模型預測的效果，且是全面的定量分析指標。R軟件提供的nricens包和PredictABEL包能方便地實現計算。

ROC曲線與AUC值刻畫了模型在選用不同診斷截點時的靈敏度和特異度，AUC的提高通常可作為模型預測能力改善的第一判斷標準，然而實際應用中，當兩個模型的AUC曲線存在交叉時，或AUC的增量極小時，僅僅采用AUC可能就很難說明問題了，此時NRI和IDI指數則是更好的選擇。當設置了一至多個概率截點，關心截點處的分類情況時， category法的NRI指數就更合適，如果目前還無法確定明確的截點，或需整體評價兩模型的預測能力，則IDI指數和AUC更適合[11]。

本文中雖然給出了兩種計算NRI指數的方法，但是實際應用中通常采用category法，即給出概率截點劃分預測概率區(qū)間，并以此作為個體分類的依據。由于category法NRI指數的計算與設定的概率截點有關，因此，預測概率區(qū)間的劃分應該選擇適宜的界值，需要結合臨床專業(yè)的需求來確定。區(qū)間劃分不恰當時，可能得不到統(tǒng)計上有意義的結論，或者所得結論缺少臨床實踐的指導意義[11];同時，由于計算時只考慮了已調變化，未考慮變化的等級數，因此可能存在信息缺失[12];這可能是NRI指數使用存在的缺陷。此外，在模型比較時，需注意概率截點選擇的側重點是否一致，比如有的模型以提升特異度為側重點，有的模型以提升靈敏度為側重點，這就需要在使用NRI和IDI指數的過程時考慮兩個模型的“側重類型”是否相同，側重類型不同的模型間沒有可比性。當前研究大多數研究考慮靈敏度和特異度各占一半的比重，但某些特殊的模型例如用于初篩的模型則偏重靈敏度，因此在模型的評價上需更加謹慎的應用NRI和IDI指數[17]。對于IDI指數來說，在計算總改善情況時，通?；疾〗M和不患病組占有同等的比重，但對于側重于靈敏度的模型，計算時應給予患病組概率改變量更多的權重。對此類模型，可以使用調整過權重后的IDI。有興趣的讀者可做進一步的研究[18]。

理論上，除了將NRI和IDI指數應用在Logistic回歸模型外，也可應用于其他可以輸出預測概率的模型，如決策樹、隨機森林、神經網絡等。但是對于這兩個指標在這些模型中的適用性，仍需進一步的研究與討論。

綜上，NRI和IDI指數是近年來新興的模型比較評價指標，在模型比較中能夠輕松處理AUC等不能解決的復雜情形，具有獨特的優(yōu)勢，但二者在腫瘤流行病學研究中應用還較少。本文結合臨床腫瘤風險預測實例，詳細闡述了兩種指標的計算原理和適用情況，以及使用R語言計算兩種指標估計值及假設檢驗的方法，以期為臨床腫瘤研究工作者提供更為詳盡的參考。當涉及兩模型優(yōu)劣的比較，或考察某一因素的加入對模型預測效果的改善程度時，除了使用約登指數及AUC外，建議研究者不妨同時考慮NRI和IDI指數，以期更加全面地展示模型的改善情況。

作者聲明：本文全部作者對于研究和撰寫的論文出現的不端行為承擔相應責任；并承諾論文中涉及的原始圖片、數據資料等已按照有關規(guī)定保存，可接受核查。

學術不端：本文在初審、返修及出版前均通過中國知網(CNKI)科技期刊學術不端文獻檢測系統(tǒng)的學術不端檢測。

同行評議：經同行專家雙盲外審，達到刊發(fā)要求。

利益沖突：所有作者均聲明不存在利益沖突。

文章版權：本文出版前已與全體作者簽署了論文授權書等協議。