R語言回歸中的Hosmer-Lemeshow擬合優(yōu)度檢驗

拓端數(shù)據(jù) 2020-03-24

展開全文

原文鏈接：http:///?p=6166

在依賴模型得出結(jié)論或預(yù)測未來結(jié)果之前，我們應(yīng)盡可能檢查我們假設(shè)的模型是否正確指定。也就是說，數(shù)據(jù)不會與模型所做的假設(shè)沖突。對于二元結(jié)果，邏輯回歸是最流行的建模方法。在這篇文章中，我們將看一下 Hosmer-Lemeshow邏輯回歸的擬合優(yōu)度檢驗。

Hosmer-Lemeshow擬合優(yōu)度檢驗

Hosmer-Lemeshow擬合優(yōu)度檢驗是基于根據(jù)預(yù)測的概率或風(fēng)險將樣本分開。具體而言，基于估計的參數(shù)值，對于樣本中的每個觀察，基于每個觀察的協(xié)變量值計算概率。

然后根據(jù)樣本的預(yù)測概率將樣本中的觀察分成g組（我們回過頭來選擇g）。假設(shè)（通常如此）g = 10。然后第一組由具有最低10％預(yù)測概率的觀察組成。第二組由預(yù)測概率次之小的樣本的10％等組成。

在實踐中，只要我們的一些模型協(xié)變量是連續(xù)的，每個觀測將具有不同的預(yù)測概率，因此預(yù)測的概率將在我們形成的每個組中變化。為了計算我們預(yù)期的觀察數(shù)量，Hosmer-Lemeshow測試取組中預(yù)測概率的平均值，并將其乘以組中的觀察數(shù)。測試也執(zhí)行相同的計算，然后計算Pearson擬合優(yōu)度統(tǒng)計量

選擇組的數(shù)量

就我所見，關(guān)于如何選擇組數(shù)g的指導(dǎo)很少。Hosmer和Lemeshow的模擬結(jié)論是基于使用的，建議如果我們在模型中有10個協(xié)變量。

直觀地說，使用較小的g值可以減少檢測錯誤規(guī)范的機(jī)會。

R

首先，我們將使用一個協(xié)變量x模擬邏輯回歸模型中的一些數(shù)據(jù)，然后擬合正確的邏輯回歸模型。

n < - 100
x < - rnorm（n）
xb < - x
pr < - exp（xb）/（1 + exp（xb））
y < - 1 *（runif（n）
mod < - glm（y~x，family = binomial）

接下來，我們將結(jié)果y和模型擬合概率傳遞給hoslem.test函數(shù)，選擇g = 10組：

Hosmer and Lemeshow goodness of fit (GOF) test data: mod$y, fitted(mod) X-squared = 7.4866, df = 8, p-value = 0.4851

這給出p = 0.49，表明沒有合適的不良證據(jù)。我們還可以從我們的hl對象中獲得一個觀察到的與預(yù)期的表：

cbind(hl$observed,hl$expected)
y0 y1 yhat0 yhat1
[0.0868,0.219] 8 2 8.259898 1.740102
(0.219,0.287] 7 3 7.485661 2.514339
(0.287,0.329] 7 3 6.968185 3.031815
(0.329,0.421] 8 2 6.194245 3.805755
(0.421,0.469] 5 5 5.510363 4.489637
(0.469,0.528] 4 6 4.983951 5.016049
(0.528,0.589] 5 5 4.521086 5.478914
(0.589,0.644] 2 8 3.833244 6.166756
(0.644,0.713] 6 4 3.285271 6.714729
(0.713,0.913] 1 9 1.958095 8.041905

為了幫助我們理解計算，現(xiàn)在讓我們自己手動執(zhí)行測試。首先，我們計算模型預(yù)測概率，然后根據(jù)預(yù)測概率的十分位數(shù)對觀測值進(jìn)行分類：

<- mod$fitted\npihatcat <- cut(pihat, brks=c(0,quantile(pi 1,0.9,0.1)),1), els=FALSE)","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet">pihat <- mod$fitted pihatcat <- cut(pihat, brks=c(0,quantile(pi 1,0.9,0.1)),1), els=FALSE)

接下來，我們循環(huán)通過組1到10，計算觀察到的0和1的數(shù)量，并計算預(yù)期的0和1的數(shù)量。為了計算后者，我們找到每組中預(yù)測概率的均值，并將其乘以組大小，這里是10：

<- array(0, dim=c(10,2))\nexpevents <- array(0, dim=c(10,2))\nobsevents <- array(0, dim=c(10,2))\n\nfor (i in 1:10) {\n\tmeanprobs[i,1] <- mean(pihat[pihatcat==i])\n \n\tobsevents[i,2] <- sum(1-y[pihatcat==i])\n}","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet">meanprobs <- array(0, dim=c(10,2)) expevents <- array(0, dim=c(10,2)) obsevents <- array(0, dim=c(10,2)) for (i in 1:10) { meanprobs[i,1] <- mean(pihat[pihatcat==i]) obsevents[i,2] <- sum(1-y[pihatcat==i]) }

最后，我們可以通過表格的10x2單元格中的（觀察到的預(yù)期）^ 2 /預(yù)期的總和來計算Hosmer-Lemeshow檢驗統(tǒng)計量：


[1] 7.486643

與hoslem.test函數(shù)的測試統(tǒng)計值一致。

改變組的數(shù)量
接下來，讓我們看看測試的p值如何變化，因為我們選擇g = 5，g = 6，直到g = 15。我們可以通過一個簡單的for循環(huán)來完成：

for（i in 5:15）{ print（hoslem.test（mod $ y，fits（mod），g = i）$ p.value） }

[1] 0.4683388
[1] 0.9216374
[1] 0.996425
[1] 0.9018581
[1] 0.933084
[1] 0.4851488
[1] 0.9374381
[1] 0.9717069
[1] 0.5115724
[1] 0.4085544
[1] 0.8686347

雖然p值有所改變，但它們都顯然不重要，所以他們給出了類似的結(jié)論，沒有證據(jù)表明不合適。因此，對于此數(shù)據(jù)集，選擇不同的g值似乎不會影響實質(zhì)性結(jié)論。

通過模擬檢查Hosmer-Lemeshow測試

要完成，讓我們進(jìn)行一些模擬，以檢查Hosmer-Lemeshow測試在重復(fù)樣本中的表現(xiàn)。首先，我們將從先前使用的相同模型重復(fù)采樣，擬合相同（正確）模型，并使用g = 10計算Hosmer-Lemeshow p值。我們將這樣做1000次，并將測試p值存儲在一個數(shù)組中：

< - array（0,1000）\n\nfor（i in 1：1000）{\n\tn < - 100\n\tx < - rnorm（n）\n \tpr < - exp（xb）/（1 + exp（xb））\n \tmod < - glm（y~x，family = binomial）\n }","classes":{"has":1}}" data-cke-widget-upcasted="1" data-cke-widget-keep-attr="0" data-widget="codeSnippet">pvalues < - array（0,1000） for（i in 1：1000）{ n < - 100 x < - rnorm（n） pr < - exp（xb）/（1 + exp（xb）） mod < - glm（y~x，family = binomial） }

完成后，我們可以計算出p值小于0.05的比例。由于此處正確指定了模型，因此我們希望這種所謂的類型1錯誤率不大于5％：


[1] 0.04

因此，在1,000次模擬中，Hosmer-Lemeshow測試在4％的情況下給出了顯著的p值，表明不合適。所以測試錯誤地表明在我們預(yù)期的5％限制內(nèi)不合適 - 它似乎工作正常。

現(xiàn)在讓我們改變模擬，以便我們適合的模型被錯誤地指定，并且應(yīng)該很難適應(yīng)數(shù)據(jù)。希望我們會發(fā)現(xiàn)Hosmer-Lemeshow測試在5％的時間內(nèi)正確地找到了不合適的證據(jù)。具體來說，我們現(xiàn)在將生成跟隨具有協(xié)變量的邏輯模型，但我們將繼續(xù)使用線性協(xié)變量擬合模型，以便我們的擬合模型被錯誤地指定。

我們發(fā)現(xiàn)，計算p值小于0.05的比例