rms | 如何绘制模型带置信区间的预测曲线
前段时间,有读者咨询如何绘制带置信区间的预测曲线。本篇介绍的rms
工具包可以很简便地解决这个问题。
该包是Regression Modeling Strategies的附加学习资源,小编从网上找到了这本书的目录:https://link.springer.com/content/pdf/bfm%3A978-3-319-19425-7%2F1.pdf。
1 lm()
/ols()
在rms
工具包中,使用Predict()
函数预测的模型结果可以直接使用plot()
或ggplot()
函数进行绘制,并自带置信区间,但是建模函数必须使用rms
工具包中的相关函数。比如对线性模型而言,基础包stats
中的函数是lm()
,而rms
工具包中是ols()
函数。
ols()
函数的语法结构如下:
ols(formula, data=environment(formula),weights, subset, na.action=na.delete, method="qr", model=FALSE,x=FALSE, y=FALSE, se.fit=FALSE, linear.predictors=TRUE,penalty=0, penalty.matrix, tol=1e-7, sigma,var.penalty=c('simple','sandwich'), ...)
可以看出,
ols()
函数和lm()
函数的用法基本上是一致的。
library(rms)
model.1 <- lm(mpg ~ wt, data = mtcars)
model.2 <- ols(mpg ~ wt, data = mtcars)summary(model.1)
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10
print(model.2)
## Linear Regression Model
##
## ols(formula = mpg ~ wt, data = mtcars)
##
## Model Likelihood Discrimination
## Ratio Test Indexes
## Obs 32 LR chi2 44.73 R2 0.753
## sigma3.0459 d.f. 1 R2 adj 0.745
## d.f. 30 Pr(> chi2) 0.0000 g 5.822
##
## Residuals
##
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
##
## Coef S.E. t Pr(>|t|)
## Intercept 37.2851 1.8776 19.86 <0.0001
## wt -5.3445 0.5591 -9.56 <0.0001
两个函数的模型参数估算结果也是一样的。
2 predict()
/Predict()
在stats
包中,与模型预测相关的有两个函数:fitted()
和predict()
。前者是根据模型结果对原数据中的样本的响应变量(因变量)进行预测:
head(fitted(model.1))
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## 23.28261 21.91977 24.88595 20.10265
## Hornet Sportabout Valiant
## 18.90014 18.79325
predict()
函数则可以对新数据的样本进行预测:
## S3 method for class 'lm'
predict(object, newdata, se.fit = FALSE, scale = NULL, df = Inf,interval = c("none", "confidence", "prediction"),level = 0.95, type = c("response", "terms"),terms = NULL, na.action = na.pass,pred.var = res.var/weights, weights = 1, ...)
object:模型对象;
newdata:新数据;必须为数据框结构;若省略,则默认为原数据;
interval:区间类型,有两种类型:置信区间(confidence)和预测区间(prediction);
level:置信水平;默认为95%。
关于置信区间和预测区间的区别,小编不是很清楚,有兴趣的可以参见这个链接的解释:https://blog.sciencenet.cn/blog-54276-288414.html。
predict(model.1, data.frame(wt = seq(1,6,0.5)),interval = "prediction")
## fit lwr upr
## 1 31.940655 25.135231 38.74608
## 2 29.268419 22.654123 35.88271
## 3 26.596183 20.128114 33.06425
## 4 23.923947 17.554110 30.29378
## 5 21.251711 14.929874 27.57355
## 6 18.579476 12.254262 24.90469
## 7 15.907240 9.527355 22.28712
## 8 13.235004 6.750452 19.71956
## 9 10.562768 3.925916 17.19962
## 10 7.890533 1.056933 14.72413
## 11 5.218297 -1.852790 12.28938predict(model.1, data.frame(wt = seq(1,6,0.5)),interval = "confidence")
## fit lwr upr
## 1 31.940655 29.18042 34.700892
## 2 29.268419 27.02030 31.516535
## 3 26.596183 24.82389 28.368481
## 4 23.923947 22.55284 25.295059
## 5 21.251711 20.12444 22.378987
## 6 18.579476 17.43342 19.725534
## 7 15.907240 14.49018 17.324295
## 8 13.235004 11.40347 15.066543
## 9 10.562768 8.24913 12.876406
## 10 7.890533 5.06154 10.719526
## 11 5.218297 1.85595 8.580644
而在rms
工具包中,模型预测对应的是Predict()
函数(首字母大写):
Predict(object, ..., fun=NULL, funint=TRUE,type = c("predictions", "model.frame", "x"),np = 200, conf.int = 0.95,conf.type = c("mean", "individual","simultaneous"),usebootcoef=TRUE, boot.type=c("percentile", "bca", "basic"),posterior.summary=c('mean', 'median', 'mode'),adj.zero = FALSE, ref.zero = FALSE,kint=NULL, ycut=NULL, time = NULL, loglog = FALSE,digits=4, name,factors=NULL, offset=NULL)
conf.type:区间类型;置信区间(mean)、预测区间(individual)。
Predict()
函数的“新数据”使用...
表示,直接写出参与预测的变量即可:
predict.2 <- Predict(model.2, wt = seq(1,6,0.5))
predict.2
## wt yhat lower upper
## 1 1.0 31.940655 29.18042 34.700892
## 2 1.5 29.268419 27.02030 31.516535
## 3 2.0 26.596183 24.82389 28.368481
## 4 2.5 23.923947 22.55284 25.295059
## 5 3.0 21.251711 20.12444 22.378987
## 6 3.5 18.579476 17.43342 19.725534
## 7 4.0 15.907240 14.49018 17.324295
## 8 4.5 13.235004 11.40347 15.066543
## 9 5.0 10.562768 8.24913 12.876406
## 10 5.5 7.890533 5.06154 10.719526
## 11 6.0 5.218297 1.85595 8.580644
##
## Response variable (y): mpg
##
## Limits are 0.95 confidence limits
可以看出,
predict.2
与predict()
函数预测的置信区间一致。
也可以不指定参与预测变量的数值,但需要提前定义变量的范围和等分的个数。指定范围的方法如下:
ddist <- datadist(mtcars)
options(datadist = "ddist")
等分数目由np
参数确定,默认值为200:
Predict(model.2, wt, np = 10)
## wt yhat lower upper
## 1 1.736000 28.007124 25.989733 30.02451
## 2 2.131194 25.895018 24.237593 27.55244
## 3 2.526389 23.782913 22.429583 25.13624
## 4 2.921583 21.670807 20.520507 22.82111
## 5 3.316778 19.558702 18.453202 20.66420
## 6 3.711972 17.446596 16.210345 18.68285
## 7 4.107167 15.334491 13.837242 16.83174
## 8 4.502361 13.222385 11.388690 15.05608
## 9 4.897556 11.110280 8.898861 13.32170
## 10 5.292750 8.998174 6.385598 11.61075
##
## Response variable (y): mpg
##
## Limits are 0.95 confidence limits
3 plot()
/ggplot()
Predict()
函数预测的结果可以直接使用plot()
函数和ggplot()
函数进行绘制。在加载rms
工具包时会自动加载ggplot2
工具包。
需要注意的是,rms
包中的plot()
函数采用的是lattice绘图系统,而非基础绘图系统(graphics
),因此不能叠加graphics
包中的函数。
plot(predict.2)
以下代码会报错:
plot(predict.2)
points(mtcars$wt, mtcars$mpg)Error in plot.xy(xy.coords(x, y), type = type, ...) : plot.new has not been called yet
在ggplot2绘图系统中,不需要额外使用几何图形函数:
ggplot(predict.2)
还可以使用图层叠加功能:
ggplot(predict.2) +geom_point(data = mtcars, aes(wt, mpg))
4 其他情况的简要说明
如果要向上述模型中添加二次项,在lm()
函数中可以使用I(wt^2)
或poly(wt,2)
的形式,但若想使用Predict()
函数对模型进行预测,就必须使用rms
工具包中的函数进行建模,该包添加二次项的函数为pol()
。
以下代码会报错:
ols(mpg ~ wt + I(wt^2), data = mtcars)Error in if (!length(fname) || !any(fname == zname)) { : 需要TRUE/FALSE值的地方不可以用缺少值
以下为正确的写法:
model.3 <- ols(mpg ~ pol(wt,2), data = mtcars)
predict.3 <- Predict(model.3, wt = seq(1,6,0.5))
predict.3ggplot(predict.3) +geom_point(data = mtcars, aes(wt, mpg))
如果模型中有多个自变量,可以使用全部变量进行预测,也可以只使用个别变量进行预测。
使用全部变量进行预测:
ddist <- datadist(mtcars)
options(datadist = "ddist")
model.4 <- ols(mpg ~ wt + drat, data = mtcars)
predict.41 <- Predict(model.4, np = 10)
plot(predict.41)
使用单个变量进行预测:
predict.42 <- Predict(model.4, wt = seq(1,6,0.5))
plot(predict.42)
未参与预测的变量实际上取的是中位数:
median(mtcars$drat)
## [1] 3.695
也可以指定未参与预测变量的取值:
predict.43 <- Predict(model.4, wt = seq(1,6,0.5),drat = 4,np = 10)
plot(predict.43)
rms
包中的模型函数基本上是自成体系的,如logistic回归用lrm()
函数,以及一系列的样条曲线函数,如rcs()
函数(linear tail-restricted cubic spline function)。
ddist <- datadist(mtcars)
options(datadist = "ddist")
model.5 <- ols(mpg ~ rcs(wt) + drat, data = mtcars)
predict.5 <- Predict(model.5, wt)
plot(predict.5)
rms | 如何绘制模型带置信区间的预测曲线相关推荐
- matlab画置信区间图,matlab绘制带置信区间的双y轴图形 | 学步园
matlab的双y轴网上有很多方法,但是带置信区间的双y轴就很少了,并且由于网上给的例子一般都是使用红蓝两色,对于只想使用黑色或者灰色的俺们来说太鲜艳啦~ 上图为使用matlab绘制的双y轴带置信区间 ...
- AI实战:搭建带注意力机制的 seq2seq 模型来做数值预测
AI实战:搭建带注意力机制的 seq2seq 模型来做数值预测 seq2seq 框架图 环境依赖 Linux python3.6 tensorflow.keras 源码搭建模型及说明 依赖库 impo ...
- R语言拟合ARIMA模型并使用拟合模型进行预测推理、使用autoplot函数可视化ARIMA模型预测结果、可视化包含置信区间的预测结果
R语言拟合ARIMA模型并使用拟合模型进行预测推理.使用autoplot函数可视化ARIMA模型预测结果.可视化包含置信区间的预测结果 目录
- Python数据分析案例-分别使用时间序列ARIMA、SARIMAX模型与Auto ARIMA预测国内汽车月销量
1. 前言 模型: ARIMA模型(英语:Autoregressive Integrated Moving Average model),差分整合移动平均自回归模型,又称整合移动平均自回归模型(移动也 ...
- Python通过ARIMA模型进行时间序列分析预测
ARIMA模型预测 时间序列分析预测就是在已有的和时间有关的数据序列的基础上构建其数据模型并预测其未来的数据,例如航空公司的一年内每日乘客数量.某个地区的人流量,这些数据往往具有周期性的规律.如下图所 ...
- Matlab用向量误差修正VECM模型蒙特卡洛Monte Carlo预测债券利率时间序列和MMSE 预测
最近我们被客户要求撰写关于VECM的研究报告,包括一些图形和统计输出. 此示例说明如何从 VEC( q ) 模型生成 Monte Carlo 预测.该示例将生成的预测与最小均方误差 (MMSE) 预测 ...
- python大气模型算法_关于预测空气质量机器学习哪个算法简单?
尽管线性模型是最简单的机器学习技术之一,但它们仍然是进行预测的强大工具. 代码格式乱可以看原文链接:http://tecdat.cn/?p=11387tecdat.cn 这尤其是由于线性模型特别容易 ...
- Matlab用向量误差修正VECM模型蒙特卡洛Monte Carlo预测债券利率时间序列和MMSE 预测...
原文链接:http://tecdat.cn/?p=27246 此示例说明如何从 VEC( q ) 模型生成 Monte Carlo 预测.该示例将生成的预测与最小均方误差 (MMSE) 预测和来自V ...
- python短期预测图_Python中利用长短期记忆模型LSTM进行时间序列预测分析
原文链接:http://tecdat.cn/?p=6663 此示例中,神经网络用于使用2011年4月至2013年2月期间的数据预测都柏林市议会公民办公室的能源消耗. 每日数据是通过总计每天提供的15分 ...
最新文章
- 互联网企业安全高级指南3.7.1 攻防驱动修改
- ASP.NET虚拟主机安全漏洞解决方案
- 独家 | 人工智能和大数据是如何联系在一起的?
- 指南--安装 NVU HTML 编辑器
- 协议簇:TCP 解析: 连接断开
- python给函数添加属性_如何在python中自动向类添加属性?
- 【SpringCloud】服务调用OpenFeign
- void* 与 shared_ptr的相互转换
- Illustrator中文版教程,如何在 Illustrator 中创建几何图案?
- python做数据分析对数学要求_Python数据分析之Pandas
- 【单目标优化求解】基于matlab粒子群混沌混合蝴蝶优化算法求解最优目标问题(HPSOBOA)【含Matlab源码 1538期】
- DoTween函数汇总
- 通过图片url 获取图片file对象
- 有没有手机版_iQOO Neo 855版性价比神机:不到两千,充电一局玩十局
- 压缩包密码,办公文档密码破解实例讲解!
- WEB常见的HTTP错误代码404 500等
- Linux基础-磁盘阵列RAID
- 怒肝半月!Python 学习路线+资源大汇总
- 用VS2010打开VS2013、VS2015、VS2017等高版本项目
- 大家好,欢迎您加入这个集体!
热门文章
- Spring Cloud与微服务学习总结(2)——Spring Cloud相较于Dubbo等RPC服务框架的优势
- Java基础学习总结(84)——Java面向对象六大原则和设计模式
- Css学习总结(1)——20个很有用的CSS技巧
- aliyun gradle 代理_android studio gradle国内代理设置
- docker mysql5.7.19_Docker19.03.13下安装Mysql57
- linux安装KVM
- Vue的数据依赖实现原理简析
- iOS网络开发(5)请求的缓存
- CentOS6离线bash漏洞—再修复方法
- Redis 2.8.9源码 - Redis中的字符串实现 sds