1. 背景知識  

Logistic回歸可以用于建立二分類結(jié)局變量的臨床預(yù)測模型，在醫(yī)學(xué)研究中可用于預(yù)測如有效/無效、發(fā)生/未發(fā)生、復(fù)發(fā)/未復(fù)發(fā)等二值化的臨床事件。預(yù)測模型有優(yōu)劣之分，好的模型不僅可以較準(zhǔn)確的預(yù)測終點(diǎn)事件發(fā)生概率（校準(zhǔn)度好），也可以很好地區(qū)分?jǐn)?shù)據(jù)集中發(fā)生終點(diǎn)事件概率不同的對象（區(qū)分度好），還可以從中發(fā)現(xiàn)某種因素對終點(diǎn)事件可能的影響（獨(dú)立風(fēng)險因素或保護(hù)因素）。因此如何判斷并驗(yàn)證模型的優(yōu)劣就顯得尤為重要。有關(guān)模型優(yōu)劣的評估指標(biāo)我們前文已經(jīng)述及，本文主要介紹Logistic回歸模型的外驗(yàn)證問題。

模型驗(yàn)證可采取Out time validation，比方說使用2005-2010年的樣本建模后，在2010-2015年的樣本中對模型進(jìn)行驗(yàn)證，這樣可以評價模型隨著時間推移，預(yù)測是否仍然準(zhǔn)確；可采取Across modeling techniques，即對于某個數(shù)據(jù)集，既可以采用邏輯回歸建模，也可以使用判別分析、決策樹、支持向量機(jī)等方式進(jìn)行建模，從中選擇在testing data中表現(xiàn)最好的模型。無論使用何種手段，使用不同于建模時所用的數(shù)據(jù)集對模型進(jìn)行外部驗(yàn)證都是非常重要的一環(huán)。

建模時，通常會先樣本數(shù)據(jù)集中抽取部分樣本用于建模，該部分樣本稱為training set。建模完成后對模型進(jìn)行驗(yàn)證時，先使用數(shù)據(jù)集中的保留樣本進(jìn)行模型的初步評價（internal validation），這部分樣本稱為testing set。在單個數(shù)據(jù)集中表現(xiàn)良好的模型在其他數(shù)據(jù)集中的表現(xiàn)不一定令人滿意，因此還需要在全新的數(shù)據(jù)集中對模型進(jìn)行外部驗(yàn)證，這個數(shù)據(jù)集稱為validation set。

2. 校準(zhǔn)度評價

下面我們使用ResourceSelection包中的hoslem.test()來進(jìn)行Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn),通常用來評價Logistic模型的校準(zhǔn)度。首先加載所需r包：

#install.packages('ResourceSelection')

library(ResourceSelection)

一、建立Logistic回歸模型 我們模擬一個數(shù)據(jù)集（training set）用于建模，這里將該數(shù)據(jù)集全部樣本拿來建模，也可以抽取部分樣本建模，剩余樣本即testing set用于對模型進(jìn)行內(nèi)部驗(yàn)證：set.seed(123)

n <- 100

x <- rnorm(n)

xb <- x

pr <- exp(xb)/(1+exp(xb))#根據(jù)邏輯回歸生成一個0-1的概率

y <- 1*(runif(n) < pr)#根據(jù)發(fā)生概率pr確定單個病人事件是否發(fā)生，pr越大，y=1的概率越大，該樣本發(fā)生終點(diǎn)事件的概率越大

intern.data <- data.frame(x=x,y=y)

mod <- glm(y~x,intern.data,family=binomial)#生成模型mod

對模型進(jìn)行Hosmer-Lemeshow test，該檢驗(yàn)將數(shù)據(jù)劃分為一定的組數(shù)g，此處參數(shù)g的含義我們在前述章節(jié)闡述校準(zhǔn)度概念時已作出解釋。假設(shè)我們預(yù)測100個人的結(jié)局發(fā)生概率，不是指我們真用模型預(yù)測出來有病/無病，模型只給我們有病的概率，根據(jù)概率大于某個截斷值(比如說0.5)來判斷有病/無病。100個人，我們最終通過模型得到了100個概率，也就是100個0-1之間的數(shù)，我們將這100個數(shù)，按照從小到大排列，再依次將這100個人分成10組，每組10個人，實(shí)際的概率其實(shí)就是這10個人中發(fā)生疾病的比例，預(yù)測的概率就是每組預(yù)測得到的10個比例的平均值，然后比較這兩個數(shù)，一個作為橫坐標(biāo)，一個作為縱坐標(biāo)，就得到了一致性曲線圖(Calibration Plot)，也可計算這個圖的95%區(qū)間范圍。在Logistic回歸模型中，校準(zhǔn)度也可以通過擬合優(yōu)度檢驗(yàn)(Hosmer-Lemeshow goodness-of-fit test)來度量。

hl <- hoslem.test(mod$y, fitted(mod), g=10)

#第一個參數(shù)為目標(biāo)事件是否真實(shí)發(fā)生，

#第二個參數(shù)為根據(jù)變量x預(yù)測的事件發(fā)生概率，

#第三個參數(shù)為分組參數(shù)g

#在Hosmer-Lemeshow檢驗(yàn)中，計算擬合優(yōu)度統(tǒng)計量，若模型擬合良好，該統(tǒng)計量則應(yīng)服

#從自由度為g-2的卡方分布。若假設(shè)檢驗(yàn)的P值<檢驗(yàn)水準(zhǔn)α，則提示模型的擬合不佳。

#輸出Hosmer-Lemeshow test的結(jié)果

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，尚不能認(rèn)為該模型擬合不良。

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檢驗(yàn)列聯(lián)表，

#y0代表未發(fā)生事件的次數(shù)，

#y1代表發(fā)生事件的次數(shù)，

#yhat0代表模型預(yù)測事件未發(fā)生的概率，

#yhat1代表模型預(yù)測事件發(fā)生概率

二、用同樣的方法模擬生成的外部數(shù)據(jù)集，對建好的模型進(jìn)行外部驗(yàn)證：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為根據(jù)logistic回歸模擬的事件發(fā)生概率

y.e <- 1*(runif(n.e) < pr.e)#y為事件實(shí)際是否發(fā)生，0為沒有發(fā)生，1為發(fā)生

exter.data <- data.frame(x=x.e,y=y.e)

#模擬好的外部數(shù)據(jù)集

#用內(nèi)部數(shù)據(jù)生成的模型預(yù)測外部數(shù)據(jù)的事件發(fā)生概率

pr.e <- predict(mod,exter.data,type = c('response'))

#使用建好模型mod預(yù)測新數(shù)據(jù)集得到的預(yù)測概率，

#第一個參數(shù)為模型對象glm，

#第二個參數(shù)為要驗(yàn)證的數(shù)據(jù)集

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，尚不能認(rèn)為模型擬合不良，提示模型在新數(shù)據(jù)集中表現(xiàn)也較好。若P<0.05，則說明模型擬合不良。

三、Hosmer-Lemeshow檢驗(yàn)統(tǒng)計量的計算原理如下：計算模型預(yù)測的概率，根據(jù)預(yù)測概率從大到小將數(shù)據(jù)分為10組，即上述函數(shù)中的參數(shù)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檢驗(yàn)統(tǒng)計量

meanprobs <- array(0, dim=c(10,2))

#建立一個10行2列的矩陣用于存放事件發(fā)生和未發(fā)生的平均概率值

expevents <- array(0, dim=c(10,2))

#建立一個10行2列的矩陣用于存放根據(jù)概率值計算的事件發(fā)生數(shù)和未發(fā)生數(shù)

obsevents <- array(0, dim=c(10,2))

#建立一個10行2列的矩陣用于存放事件的實(shí)際發(fā)生數(shù)和未發(fā)生數(shù)

for (i in 1:10) {

meanprobs[i,1] <- mean(pihat[pihatcat==i])#計算每組事件發(fā)生的平均概率

expevents[i,1] <- sum(pihatcat==i)*meanprobs[i,1]#預(yù)測事件數(shù)=組內(nèi)樣本個數(shù)*平均概率

obsevents[i,1] <- sum(y[pihatcat==i])#實(shí)際事件數(shù)

#同樣的方法計算預(yù)測事件未發(fā)生數(shù)和實(shí)際事件未發(fā)生數(shù)

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 檢驗(yàn)統(tǒng)計量。已經(jīng)證明了若模型擬合良好，該統(tǒng)計量則應(yīng)服從自由度為g-2的卡方分布。無效假設(shè)為：Hosmer-Lemeshow 檢驗(yàn)統(tǒng)計量服從自由度為g-2的卡方分布。備擇假設(shè)為：Hosmer-Lemeshow。檢驗(yàn)統(tǒng)計量不服從自由度為g-2的卡方分布。檢驗(yàn)水準(zhǔn)α=0.05hosmerlemeshow <- sum((obsevents-expevents)^2 / expevents)

hosmerlemeshow

## [1] 6.455077

計算得出的統(tǒng)計量同上述函數(shù)計算的結(jié)果一致

四、使用校正曲線圖直觀評價模型#install.packages('PredictABEL')

library(PredictABEL)

#使用plotCalibration繪制校正曲線

#參數(shù)說明：data為要驗(yàn)證的數(shù)據(jù)集，

#cOutcome指定結(jié)局變量為哪一列，

#predRisk為模型預(yù)測的發(fā)生概率可以通過predict()函數(shù)計算得出，

#groups為組數(shù)，

#rangeaxis為坐標(biāo)軸的范圍。

#現(xiàn)在在外部數(shù)據(jù)中進(jìn)行校正圖繪制

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

校正圖的橫軸為預(yù)測的風(fēng)險率，縱軸為實(shí)際風(fēng)險率，每個點(diǎn)代表一個組，其基本思想同Hosmer-Lemeshow檢驗(yàn)一致

僅僅評價模型的校準(zhǔn)度還不夠。假設(shè)某項手術(shù)后病人的短期死亡率為0.1%，現(xiàn)在有這樣一個糟糕的模型，無論病人的健康狀況如何、是否吸煙、是否有糖尿病，該模型預(yù)測的所有病人死亡率均為1%，最后術(shù)后1000例病人中有1例短期死亡。盡管模型的預(yù)測與事實(shí)相符（均為0.1%的短期死亡率）即擬合優(yōu)度較高，但這個模型并沒有將死亡風(fēng)險高的病人同死亡風(fēng)險低的病人區(qū)分開來，即模型的區(qū)分度不夠，因此我們需要評價模型的區(qū)分度（即模型將短期死亡率高的患者同短期死亡率低的患者區(qū)分開的能力）進(jìn)一步對模型進(jìn)行外部驗(yàn)證。 。

3. 區(qū)分度評價

一、使用上文擬合好的Logistic回歸模型繪制ROC曲線

#提取模型預(yù)測的概率值

pr <- predict(mod,type=c('response'))

#install.packages('pROC')

library(pROC)#使用roc函數(shù)，方程左側(cè)為實(shí)際事件發(fā)生與否，右側(cè)為模型預(yù)測的事件發(fā)生率

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

二、接下來驗(yàn)證模型在外部數(shù)據(jù)集中的區(qū)分度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

可見，模型在外部數(shù)據(jù)集中進(jìn)行驗(yàn)證，AUC=0.6723，模型在外部數(shù)據(jù)集驗(yàn)證中區(qū)分度較好。

4. 總結(jié)

以上我們總結(jié)了對于Logistic回歸模型進(jìn)行外部驗(yàn)證的方法，包括校準(zhǔn)度評估和區(qū)分度評估。

5. 參考文獻(xiàn)

[1]. https://cran./web/packages/ResourceSelection/index.html 

[2]. https://cran./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

作者：王子威；審校：周支瑞