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

在依赖模型得出结论或预测未来结果之前,我们应尽可能检查我们假设的模型是否正确指定。

由Kaizong Ye,Coin Ge撰写

也就是说,数据不会与模型所做的假设冲突。对于二元结果,逻辑回归是最流行的建模方法。在这篇文章中,我们将看一下 Hosmer-Lemeshow逻辑回归的拟合优度检验

Hosmer-Lemeshow拟合优度检验


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

×

在数据分析中,对于定类变量和低测度的定序变量,通常不能使用均值、T检验和方差分析等方法来处理。对于不符合正态分布的定类数据或低测度定序数据,其检验方法是利用交叉表技术分行分列计算交叉点的频数,利用卡方距离实施卡方检验,基于频数和数据分布形态分析不同类别的数据是否存在显著性差异,对于定类数据的对比检验,也叫独立性检验。

低测度数据
对于定类变量,其数值大小和顺序并不代表什么意义,对于定类变量和低测度的定序变量,均值和方差都不能描述变量特征,故不能通过分析其平均值、方差等参数开展数据分析。在做统计分析时,对于这类变量通常需要借助中位数、频数、百分比以及不同分布情况,实现数据描述。对于低测度数据,比较典型的研究是关于结构成分的研究,实际上是一种借助频数来分析数据分布形态,并进而发现数据分布差异性的检验。

拟合及拟合优度
由于低测度数据的特点,直接进行基于均值的检验显然是不行的,于是人们借助数学模型,提出了拟合的概念。所谓拟合,就是分析现有观测变量的分布形态,检查其分布能够与某一期望分布(或标准分布)很好地吻合起来。在数学上,拟合的过程就是寻找能很好地温和当前数据序列的数学模型的过程。为了评价拟合的程度,人们提出了判定拟合有效性的机制,这就是拟合优度。拟合优度也借助检验概率的概念来评价数据拟合的质量。

目前,对于低测度数据序列的处理最常见的分析方法是卡方检验。特别是基于交叉表的卡方检验在数据分析中具有重要的地位,它们都建立在拟合概念的基础上。另外,二项分布、游程检验等单样本检验也可以看做是数据拟合的重要应用。与此同时,对定距或定序变量的分布形态判定,也是数据拟合的应用之一,在分布形态判定过程中所获得的检验概率就是该序列与标准分布形态的拟合优度。


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

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


课程

R语言数据分析挖掘必知必会

从数据获取和清理开始,有目的的进行探索性分析与可视化。让数据从生涩的资料,摇身成为有温度的故事。

立即参加

选择组的数量


就我所见,关于如何选择组数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:


R语言Wald检验 vs 似然比检验

阅读文章


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])
}

最后,我们可以通过表格的10×2单元格中的(观察到的预期)^ 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%的情况正确地找到了不合适的证据。具体来说,我们现在将生成跟随具有协变量的逻辑模型,但我们将继续使用线性协变量拟合模型,以便我们的拟合模型被错误地指定。

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

[1] 0.648

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


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498

关注有关新文章的微信公众号


永远不要错过任何见解。当新文章发表时,我们会通过微信公众号向您推送。

技术干货

最新洞察

This will close in 0 seconds