R语言回归中的Hosmer-Lemeshow以及calibration curve校正曲线
Hosmer-Lemeshow test
详细的Hosmer-Lemeshow test使用方法的介绍请参考:看这里
calibration curve是评价模型的一个重要指标,不信的话请看文章1,文章2,总之很重要的!!!
最近需要画校正曲线,但网上的相关资料实在是太少了,花了点时间研究了一下,跟大家分享一下如何使用Hosmer-Lemeshow test来绘制calibration curve,(p.s. 都是个人理解,可能有错,欢迎讨论):
首先,按照方法介绍中的样例,我们生成一下数据:
set.seed(43657)
n <- 100
x <- rnorm(n)
xb <- x
pr <- exp(xb)/(1+exp(xb))
y <- 1*(runif(n) < pr)
这里的x是变量,y相当于标签(0,1),这里做的是二分类的。然后我们放到glm函数中fit一下:
mod <- glm(y~x, family=binomial)
然后,将结果y和模型拟合概率传递给hoslem.test函数,选择g = 10组。在这之前,要记得先安装并加载package,安装方法为install.paskages(‘ResourceSelection’),加载方法为library(ResourceSelection)。
hl <- hoslem.test(mod$y, fitted(mod), g=10)
hl
这里的mod$y和 fitted(mod)分别表示标签和预测概率,可以换成自己的数据。hl输出为:
Hosmer and Lemeshow goodness of fit (GOF) testdata: mod$y, fitted(mod)
X-squared = 7.4866, df = 8, p-value = 0.4851
我们还可以从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
其中y0, y1表示观察的结果,yhat0,yhat2表示预期的结果。为了理解这些数字都代表些什么,我们手动计算一下,首先,我们计算模型的预测概率,然后根据预测概率的十分位数对观察结果进行分组(这里其实就是将模型预测的概率进行排序,然后将数据划分成10组):
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的数量。
这里提一下,上面我们选择了g=10,那么根据模型拟合出来的概率mod$fitted和十分位数,我们将所有数据分为10组,下面是本例中10组的分组情况,可以看到每个数据在哪个组别:
10 10 9 8 2 9 2 1 6 5 2 1 1 10 5 6 4 4 8 3 6 3 9 2 9 3 6 5 4 7 4 9 3 5 8 4 7 4 3 2 2 8 6 1 3 5 7 2 7 2 8 7 6 1 1 4 8 9 10 10 1 5 1 4 6 9 10 6 7 3 8 2 2 9 5 10 7 6 5 4 7 6 9 1 10 3 3 8 10 7 5 7 10 1 8 8 4 9 3 5
为了计算期望的0和1的数量,我们找到每个组中预测概率的平均值,然后将其乘以组大小,这里为10:
meanprobs <- array(0, dim=c(10,2))
expevents <- array(0, dim=c(10,2))
obsevents <- array(0, dim=c(10,2))
stdprobs <- array(0, dim=c(10,2))
obsprobs <- array(0, dim=c(10,2))for (i in 1:10) {meanprobs[i,1] <- mean(pihat[pihatcat==i])stdprobs[i,1] <- sd(pihat[pihatcat==i])expevents[i,1] <- sum(pihatcat==i)*meanprobs[i,1]obsevents[i,1] <- sum(y[pihatcat==i])obsprobs[i,1] <- sum(y[pihatcat==i])/sum(pihatcat==i)meanprobs[i,2] <- mean(1-pihat[pihatcat==i])stdprobs[i,2] <- sd(1-pihat[pihatcat==i])expevents[i,2] <- sum(pihatcat==i)*meanprobs[i,2]obsevents[i,2] <- sum(1-y[pihatcat==i])obsprobs[i,2] <- 1-(sum(y[pihatcat==i])/sum(pihatcat==i))
}
这里meanprobs 表示每组的平均概率,expevents 表示期望事件,obsevents 表示观察事件,我加了两个参数,方便绘制calibration curve,其中,stdprobs 表示每组的标准差,obsprobs 表示每组的观察概率。
然后,我们可以通过上面计算得到的表格的10x2单元格中的(观察到的预期)^ 2 /预期的总和来计算Hosmer-Lemeshow检验统计量:
hosmerlemeshow <- sum((obsevents-expevents)^2 / expevents)
hosmerlemeshow
可以看到输出为:
[1] 7.486643
和hoslem.test函数计算出来的一样对吧。
calibration curve
现在,我们有了观察的概率表obsprobs, 和预测的平均概率meanprobs ,以及预测概率的方差stdprobs,可以绘制calibration curves了。
先绘制一个空图,
yy <- meanprobs[,c(1)]
xx <- obsprobs[,c(1)]
plot(xx, yy, type = "n", xlab = "Predicted probability", ylab = "Oberserved probability")
然后,添加数据,
points(xx, yy, type = "p", pch = 19, col = "seagreen", lty = 1)
通过改变type的值可以选择绘制散点图还是折线图,或者是别的类型的图。然后,我们在图上添加误差棒
plot_error <- function(xx, yy, sd, len = 1, col = "black") {len <- len * 0.05arrows(x0 = xx, y0 = yy, x1 = xx, y1 = yy - sd, col = col, angle = 90, length = len)arrows(x0 = xx, y0 = yy, x1 = xx, y1 = yy + sd, col = col, angle = 90, length = len)
}
plot_error(xx, yy, sd = stdprobs[,c(1)], col = "seagreen")
最后,加上对角线和legend,就可以绘制出一张校正曲线图啦。
abline(0,1,lty=3,lwd=1,col=c(rgb(0,0,0,maxColorValue=255)))labs <- c("RANDOM TEST")
legend("top", legend = labs, cex = 0.8, lty = 1, lwd = 2, pch = 19, col = c("seagreen"), inset = 0.01, horiz = TRUE, box.col = "white")
P.S 这里使用的是简单的glm模型来做模型得到的概率,如果有现成的概率和标签,也是可以直接绘制calibration curve的~
R语言回归中的Hosmer-Lemeshow以及calibration curve校正曲线相关推荐
- R语言回归中的Hosmer-Lemeshow拟合优度检验
在依赖模型得出结论或预测未来结果之前,我们应尽可能检查我们假设的模型是否正确指定.也就是说,数据不会与模型所做的假设冲突. 我们围绕回归技术进行一些咨询,帮助客户解决独特的业务问题.对于二元结果,逻辑 ...
- R语言回归模型中的Pr(>|t|)如何解读?Pr(>|t|)如何计算?
R语言回归模型中的Pr(>|t|)如何解读?Pr(>|t|)如何计算? 目录 |t|)如何解读?Pr(>|t|)如何计算?">R语言回归模型中的Pr(>|t|) ...
- R语言回归模型构建、回归模型基本假设(正态性、线性、独立性、方差齐性)、回归模型诊断、car包诊断回归模型、特殊观察样本分析、数据变换、模型比较、特征筛选、交叉验证、预测变量相对重要度
R语言回归模型构建.回归模型基本假设(正态性.线性.独立性.方差齐性).回归模型诊断.car包诊断回归模型.特殊观察样本分析.数据变换.模型比较.特征筛选.交叉验证.预测变量相对重要度 目录
- R语言回归模型协方差分析(Analysis of Covariance)
R语言回归模型协方差分析(Analysis of Covariance) 目录 R语言回归模型协方差分析(Analysis of Covariance) 输入数据 ANCOVA分析
- 在R语言环境中无法载入rJava包的解决办法
问题描述: 安装包xlsx包后,运行library("xlsx")后弹出错误窗口: RGui (64-bit): Rgui.exe - 系统错误 无法启动此程序,因为计算机中丢失 ...
- R语言回归模型残差标准误差计算实战(Residual Standard Error):计算残差标准误、残差标准误解读
R语言回归模型残差标准误差计算实战(Residual Standard Error):计算残差标准误.残差标准误解读 目录
- R语言程序设计中的for循环实战
R语言程序设计中的for循环实战 目录 R语言程序设计中的for循环实战 R语言控制结构(Control Structures) R语言程序设计中的fo
- R语言ggplot2可视化:拟合二次曲线(quadratic curve)并使用ggplot2进行可视化、可视化两个响应变量和一个预测变量的二次曲线
R语言ggplot2可视化:拟合二次曲线(quadratic curve)并使用ggplot2进行可视化.可视化两个响应变量和一个预测变量的二次曲线 目录
- R语言回归表达式中常用的符号
内容源于课程,侵删.方便大家学习R语言. 符 号 用 途 ~ 分隔符号,左边为响应变量,右边为解释变量. 例如,要通过 x, z, w 预测 y,其R代码为y ~ x ...
最新文章
- 技术宅硬核跨年,DIY墨水屏日历:自动刷新位置、天气,随机播放2000多条「毒鸡汤」| 开源...
- MyEclipse6.0注册码算法代码,MyEclipse7.0注册码算法代码
- chrome安装或者重装 【windows】
- Ubuntu下安装tilix终端仿真器
- python求123逆序数_应用Python来计算排列中的逆序数个数
- Hive 内置函数权威指南,操作大全
- stylus 设置全局样式_在vue-cli 3中, 给stylus、sass样式传入共享的全局变量
- 在nodejs项目装一个库的多个版本
- 广播的接收与U盘广播
- 什么是UKey?Ukey在密评中的应用 双因素身份认证 安当加密
- Python3 post请求上传文件
- Cadence元器件封装库
- iOS开发者账号申请
- 【云周刊】第209期:Perseus(擎天):统一深度学习分布式通信框架 [弹性人工智能]... 1
- php 极简 新闻,开源极简新闻信息分享工具
- Altium Designer 2020 学习笔记(一)-----原理图及原理图库部分(配动态图操作演示)
- iis7.0介绍和完整版如何下载
- SUSE配置zypper
- FreeRDP的编译和使用。
- margin、padding、border
热门文章
- 微信墙设计原理php,php开发|微信墙|源码|微信留言板|微信表白墙|吐槽墙|心愿墙|源码...
- 错失4000亿美金!微软做错了什么?
- 删除字符串中各单词间多余空格
- 拓扑绝缘体 量子计算机,特定条件下石墨烯可变身拓扑绝缘体
- 新帽式滑环应用安防高速球机
- eclipse取消撤销_eclipse取消撤销_Eclipse项目被删除如何撤消它?或恢复
- 课程设计-ASP.NET作文网站
- 差分算法和遗传算法的区别
- 泄露数十万用户信息,硬件钱包Ledger遭遇信任危机
- Qt: GIF图片播放器(QMovie类)