R语言回归中的Hosmer-Lemeshow拟合优度检验

原文链接：http://tecdat.cn/?p=6166

在依赖模型得出结论或预测未来结果之前，我们应尽可能检查我们假设的模型是否正确指定。也就是说，数据不会与模型所做的假设冲突。对于二元结果，逻辑回归是最流行的建模方法。在这篇文章中，我们将看一下 Hosmer-Lemeshow逻辑回归的拟合优度检验。

Hosmer-Lemeshow拟合优度检验

Hosmer-Lemeshow拟合优度检验是基于根据预测的概率或风险将样本分开。具体而言，基于估计的参数值，对于样本中的每个观察，基于每个观察的协变量值计算概率。

然后根据样本的预测概率将样本中的观察分成g组（我们回过头来选择g）。假设（通常如此）g = 10。然后第一组由具有最低10％预测概率的观察组成。第二组由预测概率次之小的样本的10％等组成。

在实践中，只要我们的一些模型协变量是连续的，每个观测将具有不同的预测概率，因此预测的概率将在我们形成的每个组中变化。为了计算我们预期的观察数量，Hosmer-Lemeshow测试取组中预测概率的平均值，并将其乘以组中的观察数。测试也执行相同的计算，然后计算Pearson拟合优度统计量

选择组的数量

就我所见，关于如何选择组数g的指导很少。Hosmer和Lemeshow的模拟结论是基于使用的，建议如果我们在模型中有10个协变量 。

直观地说，使用较小的g值可以减少检测错误规范的机会。

R

首先，我们将使用一个协变量x模拟逻辑回归模型中的一些数据，然后拟合正确的逻辑回归模型。

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

接下来，我们将结果y和模型拟合概率传递给hoslem.test函数，选择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，表明没有合适的不良证据。我们还可以从我们的hl对象中获得一个观察到的与预期的表：

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

为了帮助我们理解计算，现在让我们自己手动执行测试。首先，我们计算模型预测概率，然后根据预测概率的十分位数对观测值进行分类：

pihat <- mod$fitted
pihatcat <- cut(pihat, brks=c(0,quantile(pi 1,0.9,0.1)),1),  els=FALSE)

接下来，我们循环通过组1到10，计算观察到的0和1的数量，并计算预期的0和1的数量。为了计算后者，我们找到每组中预测概率的均值，并将其乘以组大小，这里是10：

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单元格中的（观察到的预期）^ 2 /预期的总和来计算Hosmer-Lemeshow检验统计量：

[1] 7.486643

与hoslem.test函数的测试统计值一致。

改变组的数量 接下来，让我们看看测试的p值如何变化，因为我们选择g = 5，g = 6，直到g = 15。我们可以通过一个简单的for循环来完成：

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值有所改变，但它们都显然不重要，所以他们给出了类似的结论，没有证据表明不合适。因此，对于此数据集，选择不同的g值似乎不会影响实质性结论。

通过模拟检查Hosmer-Lemeshow测试

要完成，让我们进行一些模拟，以检查Hosmer-Lemeshow测试在重复样本中的表现。首先，我们将从先前使用的相同模型重复采样，拟合相同（正确）模型，并使用g = 10计算Hosmer-Lemeshow p值。我们将这样做1000次，并将测试p值存储在一个数组中：

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值，表明不合适。所以测试错误地表明在我们预期的5％限制内不合适 - 它似乎工作正常。

现在让我们改变模拟，以便我们适合的模型被错误地指定，并且应该很难适应数据。希望我们会发现Hosmer-Lemeshow测试在5％的时间内正确地找到了不合适的证据。具体来说，我们现在将生成跟随具有协变量的逻辑模型，但我们将继续使用线性协变量拟合模型，以便我们的拟合模型被错误地指定。

[1] 0.648

我们发现，计算p值小于0.05的比例

因此，Hosmer-Lemeshow测试为我们提供了65％的不合适的重要证据。