1. 背景知識
Logistic回歸可以用於建立二分類結局變量的臨床預測模型，在醫學研究中可用於預測如有效/無效、發生/未發生、復發/未復發等二值化的臨床事件。預測模型有優劣之分，好的模型不僅可以較準確的預測終點事件發生概率（校準度好），也可以很好地區分數據集中發生終點事件概率不同的對象（區分度好），還可以從中發現某種因素對終點事件可能的影響（獨立風險因素或保護因素）。因此如何判斷並驗證模型的優劣就顯得尤為重要。有關模型優劣的評估指標我們前文已經述及，本文主要介紹Logistic回歸模型的外驗證問題。
模型驗證可採取Out time validation，比方說使用2005-2010年的樣本建模後，在2010-2015年的樣本中對模型進行驗證，這樣可以評價模型隨著時間推移，預測是否仍然準確；可採取Across modeling techniques，即對於某個數據集，既可以採用邏輯回歸建模，也可以使用判別分析、決策樹、支持向量機等方式進行建模，從中選擇在testing data中表現最好的模型。無論使用何種手段，使用不同於建模時所用的數據集對模型進行外部驗證都是非常重要的一環。
建模時，通常會先樣本數據集中抽取部分樣本用於建模，該部分樣本稱為training set。建模完成後對模型進行驗證時，先使用數據集中的保留樣本進行模型的初步評價（internal validation），這部分樣本稱為testing set。在單個數據集中表現良好的模型在其他數據集中的表現不一定令人滿意，因此還需要在全新的數據集中對模型進行外部驗證，這個數據集稱為validation set。
2. 校準度評價
下面我們使用ResourceSelection包中的hoslem.test()來進行Hosmer-Lemeshow擬合優度檢驗,通常用來評價Logistic模型的校準度。首先加載所需r包：
#install.packages("ResourceSelection")
library(ResourceSelection)
一、建立Logistic回歸模型 我們模擬一個數據集（training set）用於建模，這裡將該數據集全部樣本拿來建模，也可以抽取部分樣本建模，剩餘樣本即testing set用於對模型進行內部驗證：
set.seed(123)
n <- 100
x <- rnorm(n)
xb <- x
pr <- exp(xb)/(1+exp(xb))#根據邏輯回歸生成一個0-1的概率
y <- 1*(runif(n) < pr)#根據發生概率pr確定單個病人事件是否發生，pr越大，y=1的概率越大，該樣本發生終點事件的概率越大
intern.data <- data.frame(x=x,y=y)
mod <- glm(y~x,intern.data,family=binomial)#生成模型mod
對模型進行Hosmer-Lemeshow test，該檢驗將數據劃分為一定的組數g，此處參數g的含義我們在前述章節闡述校準度概念時已作出解釋。假設我們預測100個人的結局發生概率，不是指我們真用模型預測出來有病/無病，模型只給我們有病的概率，根據概率大於某個截斷值(比如說0.5)來判斷有病/無病。100個人，我們最終通過模型得到了100個概率，也就是100個0-1之間的數，我們將這100個數，按照從小到大排列，再依次將這100個人分成10組，每組10個人，實際的概率其實就是這10個人中發生疾病的比例，預測的概率就是每組預測得到的10個比例的平均值，然後比較這兩個數，一個作為橫坐標，一個作為縱坐標，就得到了一致性曲線圖(Calibration Plot)，也可計算這個圖的95%區間範圍。在Logistic回歸模型中，校準度也可以通過擬合優度檢驗(Hosmer-Lemeshow goodness-of-fit test)來度量。
hl <- hoslem.test(mod$y, fitted(mod), g=10)
#第一個參數為目標事件是否真實發生，
#第二個參數為根據變量x預測的事件發生概率，
#第三個參數為分組參數g
#在Hosmer-Lemeshow檢驗中，計算擬合優度統計量，若模型擬合良好，該統計量則應服
#從自由度為g-2的卡方分布。若假設檢驗的P值<檢驗水準α，則提示模型的擬合不佳。
#輸出Hosmer-Lemeshow test的結果
hl
## Hosmer and Lemeshow goodness of fit (GOF) test
## data: mod$y, fitted(mod)
## X-squared = 6.4551, df = 8, p-value = 0.5964
P值為0.5964，尚不能認為該模型擬合不良。
cbind(hl$observed,hl$expected)
## y0 y1 yhat0 yhat1
## [0.111,0.298] 7 3 7.692708 2.307292
## (0.298,0.396] 8 2 6.491825 3.508175
## (0.396,0.454] 5 5 5.764301 4.235699
## (0.454,0.494] 6 4 5.243437 4.756563
## (0.494,0.564] 7 3 4.739571 5.260429
## (0.564,0.624] 4 6 4.077834 5.922166
## (0.624,0.669] 2 8 3.532070 6.467930
## (0.669,0.744] 2 8 2.910314 7.089686
## (0.744,0.809] 1 9 2.213029 7.786971
## (0.809,0.914] 2 8 1.334912 8.665088
#生成Hosmer-Lemeshow檢驗列聯表，
#y0代表未發生事件的次數，
#y1代表發生事件的次數，
#yhat0代表模型預測事件未發生的概率，
#yhat1代表模型預測事件發生概率
二、用同樣的方法模擬生成的外部數據集，對建好的模型進行外部驗證：
set.seed(123)
n.e <- 150#.e代表external
x.e <- rnorm(n.e)#自變量x
xb.e <- x.e
pr.e <- exp(xb.e)/(1+exp(xb.e))#pr.e為根據logistic回歸模擬的事件發生概率
y.e <- 1*(runif(n.e) < pr.e)#y為事件實際是否發生，0為沒有發生，1為發生
exter.data <- data.frame(x=x.e,y=y.e)
#模擬好的外部數據集
#用內部數據生成的模型預測外部數據的事件發生概率
pr.e <- predict(mod,exter.data,type = c("response"))
#使用建好模型mod預測新數據集得到的預測概率，
#第一個參數為模型對象glm，
#第二個參數為要驗證的數據集
hl.e <- hoslem.test(y.e,pr.e, g=10)
hl.e
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: y.e, pr.e
## X-squared = 10.313, df = 8, p-value = 0.2437
P=0.2437>0.05，尚不能認為模型擬合不良，提示模型在新數據集中表現也較好。若P<0.05，則說明模型擬合不良。
三、Hosmer-Lemeshow檢驗統計量的計算原理如下：計算模型預測的概率，根據預測概率從大到小將數據分為10組，即上述函數中的參數g。
pihat <- mod$fitted
pihatcat <- cut(pihat, breaks=c(0,quantile(pihat, probs=seq(0.1,0.9,0.1)),1), labels=FALSE)
然後計算Hosmer-Lemeshow檢驗統計量
meanprobs <- array(0, dim=c(10,2))
#建立一個10行2列的矩陣用於存放事件發生和未發生的平均概率值
expevents <- array(0, dim=c(10,2))
#建立一個10行2列的矩陣用於存放根據概率值計算的事件發生數和未發生數
obsevents <- array(0, dim=c(10,2))
#建立一個10行2列的矩陣用於存放事件的實際發生數和未發生數
for (i in 1:10) {
meanprobs[i,1] <- mean(pihat[pihatcat==i])#計算每組事件發生的平均概率
expevents[i,1] <- sum(pihatcat==i)*meanprobs[i,1]#預測事件數=組內樣本個數*平均概率
obsevents[i,1] <- sum(y[pihatcat==i])#實際事件數
#同樣的方法計算預測事件未發生數和實際事件未發生數
meanprobs[i,2] <- mean(1-pihat[pihatcat==i])
expevents[i,2] <- sum(pihatcat==i)*meanprobs[i,2]
obsevents[i,2] <- sum(1-y[pihatcat==i])
}
計算 Hosmer-Lemeshow 檢驗統計量。已經證明了若模型擬合良好，該統計量則應服從自由度為g-2的卡方分布。無效假設為：Hosmer-Lemeshow 檢驗統計量服從自由度為g-2的卡方分布。備擇假設為：Hosmer-Lemeshow。檢驗統計量不服從自由度為g-2的卡方分布。檢驗水準α=0.05
hosmerlemeshow <- sum((obsevents-expevents)^2 / expevents)
hosmerlemeshow
## [1] 6.455077
計算得出的統計量同上述函數計算的結果一致
四、使用校正曲線圖直觀評價模型
#install.packages("PredictABEL")
library(PredictABEL)
#使用plotCalibration繪製校正曲線
#參數說明：data為要驗證的數據集，
#cOutcome指定結局變量為哪一列，
#predRisk為模型預測的發生概率可以通過predict()函數計算得出，
#groups為組數，
#rangeaxis為坐標軸的範圍。
#現在在外部數據中進行校正圖繪製
plotCalibration(data=exter.data,
cOutcome=2,
predRisk=pr.e,
groups=10,
rangeaxis=c(0,1))
## $Table_HLtest
## total meanpred meanobs predicted observed
## [0.111,0.259) 15 0.199 0.333 2.99 5
## [0.259,0.372) 15 0.306 0.333 4.59 5
## [0.372,0.426) 15 0.393 0.333 5.90 5
## [0.426,0.477) 15 0.449 0.267 6.73 4
## [0.477,0.536) 15 0.499 0.400 7.49 6
## [0.536,0.593) 15 0.562 0.467 8.44 7
## [0.593,0.655) 15 0.624 0.533 9.36 8
## [0.655,0.725) 15 0.685 0.733 10.28 11
## [0.725,0.808) 15 0.764 0.600 11.46 9
## [0.808,0.914] 15 0.866 0.733 12.99 11
##
## $Chi_square
## [1] 10.329
##
## $df
## [1] 8
##
## $p_value
## [1] 0.2427
校正圖的橫軸為預測的風險率，縱軸為實際風險率，每個點代表一個組，其基本思想同Hosmer-Lemeshow檢驗一致
僅僅評價模型的校準度還不夠。假設某項手術後病人的短期死亡率為0.1%，現在有這樣一個糟糕的模型，無論病人的健康狀況如何、是否吸菸、是否有糖尿病，該模型預測的所有病人死亡率均為1%，最後術後1000例病人中有1例短期死亡。儘管模型的預測與事實相符（均為0.1%的短期死亡率）即擬合優度較高，但這個模型並沒有將死亡風險高的病人同死亡風險低的病人區分開來，即模型的區分度不夠，因此我們需要評價模型的區分度（即模型將短期死亡率高的患者同短期死亡率低的患者區分開的能力）進一步對模型進行外部驗證。。
3. 區分度評價
一、使用上文擬合好的Logistic回歸模型繪製ROC曲線
#提取模型預測的概率值
pr <- predict(mod,type=c("response"))
#install.packages("pROC")
library(pROC)
#使用roc函數，方程左側為實際事件發生與否，右側為模型預測的事件發生率
roccurve <- roc(y ~ pr)
#繪製ROC曲線
plot.roc(roccurve,xlim = c(1,0),ylim=c(0,1))
#獲取ROC曲線的AUC值
auc(roccurve)
## Area under the curve: 0.7455
在intern.data中，模型的AUC為0.7455
二、接下來驗證模型在外部數據集中的區分度
pr.e <- predict(mod,newdata = exter.data, type=c("response"))
roccurve <- roc(y.e ~ pr.e)
#繪製ROC曲線
plot.roc(roccurve)
#獲取ROC曲線的AUC值
auc(roccurve)
## Area under the curve: 0.6723
可見，模型在外部數據集中進行驗證，AUC=0.6723，模型在外部數據集驗證中區分度較好。
4. 總結
以上我們總結了對於Logistic回歸模型進行外部驗證的方法，包括校準度評估和區分度評估。
5. 參考文獻
[1]. https://cran.r-project.org/web/packages/ResourceSelection/index.html
[2]. https://cran.r-project.org/web/packages/PredictABEL/index.html
[3]. Xavier Robin, Natacha Turck, Alexandre Hainard, Natalia Tiberti, Frédérique Lisacek, Jean-Charles Sanchez and Markus Müller (2011). 「pROC: an open-source package for R and S+ to analyze and compare ROC curves」. BMC Bioinformatics, 12, p. 77. DOI: 10.1186/1471-2105-12-77
作者：王子威；審校：周支瑞
作者簡介：王子威，海軍軍醫大學附屬長海醫院泌尿外科研究生在讀，R語言愛好者。
*掃描下圖二維碼或點擊「閱讀原文」查看周老師全部課程