在依靠一个模型得出结论或预测未来结果之前，我们应该尽可能地检查我们所假设的模型是否正确指定。也就是说，数据与模型的假设没有冲突。对于二元结果，逻辑回归是最流行的建模方法。在这篇文章中，我们将看看流行的、但有时受到批评的Hosmer-Lemeshow逻辑回归的拟合度检验。
Logistic回归模型
我们假定我们有二分类结局变量及协变量
. logistic回归模型假定
这表明
未知参数
通常是通过最大似然法估计的。在R中通过glm (generalized linear model)函数实现，
Hosmer-Lemeshow拟合度测试
Hosmer-Lemeshow拟合度检验的基础是根据预测的概率或风险将样本划分开来。具体来说，根据样本中每个观测值的估计参数值，根据每个观测值的协变量值，计算出Y=1的概率。
然后根据预测的概率将样本中的观察值分成g组（我们稍后再来讨论g的选择）。假设（正如通常的做法）g=10。那么第一组由预测概率最低的10%的观测值组成。第二组由预测概率次小的10%的样本组成，等等。
暂时假设第一组中的所有观测值的预测概率都是0.1。那么，如果我们的模型是正确指定的，我们将期望这些观测值中Y=1的比例为10%。当然，即使模型被正确指定，观察到的比例也会在一定程度上偏离10%，但不会偏离太多。如果这组观察值中Y=1的比例反而是90%，这就提示我们的模型没有准确预测概率（风险），即表明我们的模型没有很好地拟合数据。
在实践中，只要我们的一些模型协变量是连续的，每个观察值都会有不同的预测概率，因此在我们形成的每个组中，预测概率都会有所不同。为了计算我们会有多少个Y=1的观察值，Hosmer-Lemeshow检验取组内预测概率的平均值，然后乘以组内观察值的数量。该检验还对Y =0进行了同样的计算。 的情况下进行同样的计算，然后计算出皮尔逊拟合度统计数字
表示第l组中观察到的Y=0观察值的数量，
表示第l组中观察到的Y=1观察值的数量。
和
表示预期的0数目。
在1980年的一篇论文中，Hosmer-Lemeshow通过模拟表明，（只要p+1<g）当模型被正确指定时，他们的检验统计量近似地遵循（g-2）自由度上的卡方分布。这意味着，给定我们的拟合模型，p值可以用计算出的检验统计量作为相应的卡方分布的右尾概率来计算。如果p值较小，则表明拟合效果不佳。
应该强调的是，p值大并不意味着模型拟合得好，因为缺乏反对无效假设的证据并不等同于支持替代假设的证据。特别是，如果我们的样本量较小，测试的p值高可能只是测试检测错误规范的能力较低的结果，而不是表明拟合度好。
选择组数
Hosmer和Lemeshow从模拟中得出的结论是基于使用（g>p+1），这表明如果我们在模型中有10个协变量，我们应该选择g>11，尽管这在教科书或软件包中似乎没有提到。
直观地说，使用较小的g值应该能减少检测错误指定的机会。然而，如果我们把g选得太大，每组的数字可能会非常小，以至于很难确定观察到的和预期的差异是由于偶然性还是表明模型的错误指定。
许多人（如Paul Allison）强调的另一个问题是，对于一个给定的数据集，如果改变g，有时会得到一个完全不同的p值，这样，在选择一个g时，我们可能会得出我们的模型不适合的结论，但在选择另一个时，我们会得出没有证据表明不适合。这的确是测试的一个令人不安的方面。
Hosmer-Lemeshow在R中的实现
R的glm函数不能进行Hosmer-Lemeshow检验，但许多其他R库都有函数来进行检验。下面我说明了使用ResourceSelection库中的hoslem.test函数来做这件事.
首先，我们将从一个带有一个协变量x的逻辑回归模型中模拟一些数据，然后拟合出正确的逻辑回归模型。这意味着我们的模型是正确指定的，而且我们应该有希望检测不到拟合不良的证据。
library(ResourceSelection)
set.seed(43657)
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组。
hl <- hoslem.test(mod$y, fitted(mod), g=10)
hl
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, breaks=c(0,quantile(pihat, probs=seq(0.1,0.9,0.1)),1), labels=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])
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])
}最后，我们可以通过表格的10x2单元中的（观察-预期）^2/预期之和来计算Hosmer-Lemeshow检验统计量。
hosmerlemeshow <- sum((obsevents-expevents)^2 / expevents)
hosmerlemeshow
[1] 7.486643与hoslem.test函数的检验统计值一致。
改变组数
接下来，让我们看看当我们选择g=5、g=6，直至g=15时，测试的p值是如何变化的。我们可以用一个简单的for循环来做这件事。
for (i in 5:15) {
print(hoslem.test(mod$y, fitted(mod), g=i)$p.value)
}which gives
[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值。我们将这样做1,000次，并将测试的p值存储在一个数组中。
pvalues <- array(0, 1000)
for (i in 1:1000) {
n <- 100
x <- rnorm(n)
xb <- x
pr <- exp(xb)/(1+exp(xb))
y <- 1*(runif(n) < pr)
mod <- glm(y~x, family=binomial)
pvalues[i] <- hoslem.test(mod$y, fitted(mod), g=10)$p.value
}完成后，我们可以计算出P值小于0.05的比例。由于这里的模型是正确指定的，我们希望这个所谓的第一类错误率不大于5%。
mean(pvalues<=0.05) [1] 0.04
因此，从1000次模拟中，Hosmer-Lemeshow测试给出了一个显著的P值，表明在4%的情况下拟合度很差。因此，该测试在我们预期的5%的范围内错误地提示了拟合不良--它似乎是在工作。
现在让我们改变模拟，使我们拟合的模型被错误地指定，并且应该对数据的拟合程度很差。希望我们能发现，Hosmer-Lemeshow检验能在超过5%的情况下正确地找到拟合不良的证据。但我们将继续以线性
为协变量来拟合模型，这样我们的拟合模型就会被错误地指定。要做到这一点，我们只需将生成线性预测器的那条线改为等于
:
for (i in 1:1000) {
n <- 100
x <- rnorm(n)
xb <- x^2
pr <- exp(xb)/(1+exp(xb))
y <- 1*(runif(n) < pr)
mod <- glm(y~x, family=binomial)
pvalues[i] <- hoslem.test(mod$y, fitted(mod), g=10)$p.value
}计算P值小于0.05的比例，我们发现
mean(pvalues<=0.05)
[1] 0.648因此，Hosmer-Lemeshow检验在65%的情况下为我们提供了拟合不良的重要证据。它没有更经常地检测出不良拟合度，至少部分是功率的结果--随着样本量的增加，检测不良拟合度（或实际上任何东西！）的功率将更大。
当然，这只是一个非常简单的模拟研究，但很高兴看到，至少对于我们所使用的设置来说，该测试的表现是我们所希望的。
最后，一个评论。像Hosmer-Lemeshow这样的 "全局 "拟合度检验的一个局限性是，如果我们得到一个显著的P值，表明拟合度很差，那么该检验就不能说明模型在哪些方面拟合度很差。
补充材料解读：
Hosmer-Lemeshow统计量的计算分6个步骤进行，[2]以170名志愿者的咖啡因数据为例。
1. 计算所有n个受试者的p(成功率)
使用逻辑回归的系数计算每个受试者的p(成功率)。具有相同解释变量值的受试者将具有相同的估计成功概率。下表显示了p（成功），即逻辑模型所预测的志愿者获得A级的预期比例。
2. 将p(成功率)从最大值到最小值排序
步骤1的表格按p(成功率)排序，即预期比例。如果每个志愿者服用的剂量不同，表中就会有170个不同的值。因为只有21个独特的剂量值，所以只有21个独特的p（成功）值。
3. 3. 将有序的数值分为Q百分位数组
将p(success)的有序值划分为Q组。组的数量，Q，通常为10。由于p(success)的值是并列的，每组的受试者人数可能不一样。Hosmer-Lemeshow检验的不同软件实现采用不同的方法来处理具有相同p(成功)的受试者，因此创建Q组的切点可能不同。此外，使用不同的Q值会产生不同的切割点。步骤4中的表格显示了咖啡因数据的Q=10区间。
4. 建立一个观察到的和预期的计数表
每个区间的成功和失败的观察数是通过计算该区间的受试者而得到的。一个区间内的预期成功数是该区间内受试者的成功概率之和。
下表显示了Bilder和Loughin的R函数HLTest()选择的p(success)区间的切点，观察到的和预期的数量为A和非A。
5. 从表中计算出Hosmer-Lemeshow统计量
Hosmer-Lemeshow统计量是用介绍中给出的公式计算的，对于咖啡因的例子来说，它是17.103。
{\displaystyle H=\sum _{q=1}^{10}\left({\frac {(Observed.A-Expected.A)^{2}}{Expected.A}}+{\frac {(Observed.not.A-Expected.not.A)^{2}}{Expected.not.A}}\right)}
6. 计算P值
将计算出的Hosmer-Lemeshow统计量与自由度为Q-2的卡方分布进行比较，以计算出P值。
在咖啡因的例子中，有Q=10组，有10-2=8个自由度。df=8的卡方统计量为17.103，其p值为p=0.029。p值低于α=0.05，所以拒绝了观察到的和预期的比例在所有剂量中都是相同的这一无效假设。计算的方法是得到8个自由度的右尾卡方分布的累积分布函数，即cdf_chisq_rt（x，8），或1-cdf_chisq_lt（x，8）。
Reference：
The Hosmer-Lemeshow goodness of fit test for logistic regression – The Stats Geek
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
