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 < -  100x < -  rnorm(n)xb < -  xpr < -  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%的时间内正确地找到了不合适的证据。具体来说,我们现在将生成跟随具有协变量的逻辑模型,但我们将继续使用线性协变量拟合模型,以便我们的拟合模型被错误地指定。

 

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

[1] 0.648

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

大数据部落——中国专业的第三方数据服务提供商,提供定制化的一站式数据挖掘和统计分析咨询服务

统计分析和数据挖掘咨询服务 :y0.cn/teradat(咨询服务请联系官网客服

QQ:3025393450

 

【服务场景】

科研项目;

公司项目外包 ;线上线下一对一培训 ;数据采集;学术研究;报告撰写;市场调查。

【大数据部落】提供定制化的一站式数据挖掘和统计分析咨询服务

分享最新的大数据资讯,每天学习一点数据分析,让我们一起做有态度的数据人

微信客服号:lico_9e

QQ交流群:186388004 

欢迎选修我们的R语言数据分析挖掘必知必会课程!

欢迎关注微信公众号,了解更多数据干货资讯!

原文地址:https://www.cnblogs.com/tecdat/p/11448145.html

时间: 2024-08-24 23:12:11

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

在R语言环境中无法载入rJava包的解决办法

问题描述: 安装包xlsx包后,运行library("xlsx")后弹出错误窗口: RGui (64-bit): Rgui.exe - 系统错误 无法启动此程序,因为计算机中丢失 jvm.dll.尝试重新安装该程序以解决此问题. 在R语言环境中的错误是: 载入需要的程辑包:rJava Error : loadNamespace()里算'rJava'时.onLoad失败了,详细内容: 调用: inDL(x, as.logical(local), as.logical(now), ...)

R语言学习中的小bug:R中矩阵相乘错误于A %*% B: 需要数值/复数矩阵/矢量参数

遇到了小bug: R中矩阵相乘错误于A %*% B: 需要数值/复数矩阵/矢量参数 看到网上别人的做法,发现了用class(A)和class(B)之后才发现,是因为读入的时候数据的类型不对,A.B的类型并不是matrix,才导致了这个问题. 用as.matrix来变型一下,就OK了. R语言学习中的小bug:R中矩阵相乘错误于A %*% B: 需要数值/复数矩阵/矢量参数,布布扣,bubuko.com

R语言编程中的常见错误

R语言编程中的常见错误有一些错误是R的初学者和经验丰富的R程序员都可能常犯的.如果程序出错了,请检查以下几方面.? 使用了错误的大小写.help().Help()和HELP()是三个不同的函数(只有第一个是正确的).? 忘记使用必要的引号.install.packages("gclus")能够正常执行,然而Install.packages(gclus)将会报错.? 在函数调用时忘记使用括号.例如,要使用help()而非help.即使函数无需参数,仍需加上().? 在Windows上,路

R语言putty中直接使用X11(Xming)绘图

1.下载Xming地址 http://pan.baidu.com/s/1o6ilisU,安装,推荐默认安装在C盘,推荐快捷方式放在与putty快捷方式同一个文件夹: 2.打开putty,在SSH的X11选项中勾选Enable X11 forwarding,保存putty设置: 3.运行Xming,最小化在系统托盘,使用putty进入一个服务器: 4.打开R语言 > x=c(-5:5)> y=x*x> plot(x,y,type='b') 5.绘图显示曲线

R语言并行计算中的内存控制

R语言使用向量化计算,因此非常容易在集群上进行并行计算.parallel 包提供了非常方便的函数用来进行并行计算,但有一个问题是并行时对于内存中的对象会拷贝多份,因此会比较占内存,这里提供一个比较简易的方法在内存中共享对象从而达到降低内存占用的目的. cl<-makeCluster(10, type="FORK") result_list <- parLapply(cl, list, function) stopCluster(cl) 非常简单,在创建集群的时候添加type

R语言 回归的多面性

回归是一个令人困惑的词,因为它有许多特殊变种(见表8-1).对于回归模型的拟合,R提 供的强大而丰富的功能和选项也同样令人困惑.例如,2005年Vito Ricci创建的列表表明,R中做 回归分析的函数已超过了205个(http://cran.r-project.org/doc/contrib/Ricci-refcardregression.pdf). 表8-1 回归分析的各种变体 回归类型 用 途 简单线性 用一个量化的解释变量预测一个量化的响应变量 多项式 用一个量化的解释变量预测一个量化的

R语言-回归

定义: 回归是统计学的核心,它其实是一个广义的概念,通常指那些用一个或多个预测变量来预测响应变量.既:从一堆数据中获取最优模型参数 1.线性回归 1.1简单线性回归 案例:女性预测身高和体重的关系 结论:身高和体重成正比关系 1 fit <- lm(weight ~ height,data = women) 2 summary(fit) 3 plot(women$height,women$weight,xlab = 'Height inches',ylab = 'Weight pounds')

R | R语言表达式中常用的符号

符号 用途 ~ 分隔符号,左边为响应变量,右边为解释变量,eg:要通过x.z和w预测y,代码为y~x+z+w + 分隔预测变量 : 表示预测变量的交互项 eg:要通过x.z及x与z的交互项预测y,代码为y~x+z+x:z * 表示所有可能交互项的简洁方式,代码y~x*z*w可展开为y~x+z+w+x:z+x:w+z:w+x:z:w ^ 表示交互项达到某个次数,代码y~(x+z+w)^2可展开为y~x+z+w+x:z+x:w+z:w . 表示包含除因变量外的所有变量,eg:若一个数据框包含变量x.

R语言数据挖掘中的,“回归分析”是如何操作的?

回归分析是对多个自变量(又称为预测变量)建立一个函数来预测因变量(又称为响应变量的值). 例如,银行根据房屋贷款申请人的年龄.收入.开支.职业.负担人口,以及整体信用限额等因素,来评估申请人的房贷风险. 线性回归 线性回归是利用预测变量的一个线性组合函数,来预测响应变量的统计分析方法,该线性回归模型的形式如下: y = c0 + c1x1 + c2x2 + -+ ckxk; x1, x2,- xk为预测变量,y为对预测的响应变量. 下面将在澳大利亚消费者价格指数(CPI)的数据上使用函数lm做线