R语言广义可加模型在空气环境污染方面的应用(1)
粉丝私信我希望复制一篇文章的图片,图片来源于文章:Wu C, Yan Y, Chen X, Gong J, Guo Y, Zhao Y, Yang N, Dai J, Zhang F, Xiang H. Short-term exposure to ambient air pollution and type 2 diabetes mortality: A population-based time series study. Environ Pollut. 2021 Nov 15;289:117886. doi: 10.1016/j.envpol.2021.117886. Epub 2021 Jul 31. PMID: 34371265.
文章有3个图片,是一个关于空气污染和糖尿病发病率的图片,图形如下:
文章内容比较多,我打算通过两节内容来演示,今天我们来演示一下文章1-2的图片生成,我这里没有环境污染和糖尿病的数据,使用的既往的美国芝加哥1987年至 2000年大气污染与死亡数据(公众号回复:芝加哥2,可以获得数据)做实例分析,我们先导入需要的R包和数据看看
library(tsModel)
library(ggplot2)
library(nlme)
library(mgcv)
bc<-read.csv("E:/r/test/chicago.csv",sep=',',header=TRUE)
我们先来看看数据的构成,death:死亡人数 (per day),pm10:大气污染物pm10的中位数值,pm25median,o3median:二氧化硫的中位数值,time:天数,这里就是我们的时间,tmpd:华氏温度,date:日期
在文章中,作者在图1分析了多个指标和糖尿病的关联,我这里先绘制这个时间相关的折线图
先把日期变量转换一下格式
bc$date<-as.Date(bc$date)
我们先绘制一个pm10的图作者是以间隔1年为坐标,我这里数据年份多一点,这里以2年为间隔
ggplot(bc, aes(x =date, y = pm10,color="red"))+geom_line(size=1)+scale_x_date(date_labels = "%Y",date_breaks = "2 year")+xlab("Year")+ ylab("pm10")+theme_bw()+theme( axis.title=element_text(size=10,face="plain",color="black"),axis.text = element_text(size=10,face="plain",color="black"),legend.position = c(0.8,0.8),legend.background = element_blank())
这样单变量的时间折线图就绘制好了,多个变量就把它组合起来,先把长数据转成宽数据,这里收集了"pm10",“o3”,"rhum"这3个指标rhum我也不知道是什么,
library(reshape2)
dat<-melt(bc,id=c("date","time"),measure.vars = (c("pm10","o3","rhum")),variable.name = "measure",value.name = "value")
转换好以后就可以绘图了
ggplot(dat, aes(x =date, y = value,color="red"))+geom_line(aes(group=measure,col=measure),size=1)+scale_x_date(date_labels = "%Y",date_breaks = "2 year")+facet_wrap(~measure)+xlab("Year")+ ylab("pm10")+theme_bw()+theme( axis.title=element_text(size=10,face="plain",color="black"),axis.text = element_text(size=10,face="plain",color="black"),legend.position = c(0.8,0.8),legend.background = element_blank())
稍微调整一下就生成好图片了
这样时间相关折线图就绘制好了,接下来作者图2绘制了一个在不同年龄分层分析中,T2DM死亡率与污染物浓度增加10μg/m 3相关的95%CI变化百分比(%),
这其实就是个带误差和可信区间的折线图,数据是以年龄来分层,我这里没有分层变量,我自己生成一个fage变量,0表示低龄,1表示高龄
set.seed(1234)
bc$fage<-sample(0:1,size=5114,replace=TRUE)
bc$fage<-as.factor(bc$fage)
在文章中
我们先来生成一个pm10lag01变量,赋值为过去两天大气PM10浓度的移动均值,表示图中的lag01
pm10lag<-runMean(bc.f$pm10,0:i)###产生 pm10lag01变量,赋值为过去两天大气PM10浓度的移动均值
搞好变量后就可以建立模型了,文章中作者已经给出它的模型,以及模型的详细解释,我这里就直接上代码了,dow这里是星期几,作者转成了一个分类变量,就是周末和非周末,我这里就不弄了
fit1<-gam(death~pm10lag01+ns(temp,3)+ns(o3,3)+ns(date,7*14),family = quasipoisson(),data=bc) #GAM 模型拟合
summary(fit1)
这部分的操作我在既往文章《R语言mgcv包时间序列分析在空气污染与健康领域的应用(1)》有介绍,有兴趣的可以自己看一下,我这里直接上代码了
b<-as.numeric(summary(fit1)$ coeff[2,1])#提取系数
se<-as.numeric(summary(fit1)$ coeff[2,2]) #提取标准误
ER<-(exp(b*10)-1)*100 ####计算 PM 10每升高 10μg /m 3 ,死亡的超额危险度ER
ERlp<-(exp((b-1.96*se)* 10)-1)*100 #计算95%CI
ERup<-(exp((b+1.96*se)* 10)-1)*100 #计算95%CI
这样我们我们就制作出了lag01的线段图数据,作者用了单日滞后模型和多日滞后模型,分别是lag0—lag7我这里制作多日滞后模型,制作8天需要写一个循环。
其实就是把前面的过程整合起来
for (i in 0:7) {dat<-NULLpm10lag<-runMean(bc$pm10,0:i)###产生 pm10lag01变量,赋值为过去两天大气PM10浓度的移动均值bc$pm10lag<-pm10lagfit<-gam(death~pm10lag+ns(temp,3)+ns(o3,3)+ns(date,7*14),family = quasipoisson(),data=bc) #GAM 模型拟合 b<-as.numeric(summary(fit)$coeff[2,1])#提取系数se<-as.numeric(summary(fit)$ coeff [2,2]) #提取标准误ER<-(exp(b*10)-1)*100 ####计算 PM 10每升高 10μg /m 3 ,死亡的超额危险度ERERlp<-(exp((b-1.96*se)* 10)-1)*100 #计算95%CIERup<-(exp((b+1.96*se)* 10)-1)*100 #计算95%CIlag<- paste0(i)d<-data.frame(lag=lag,ER=ER,se=se,ERlp=ERlp,ERup=ERup)dat<-rbind(dat,d)
}
最后的出dat就是lag0—lag7的数据
得出以后就可以绘图了
pd <- position_dodge(0.001)ggplot(dat, aes(x=lag, y=ER)) + geom_errorbar(aes(ymin=ERlp,ymax=ERup),width=.1,position=pd) +geom_line(position=pd) +geom_point(position=pd)
这样图就完成了,如果画分类怎么画呢,就是对数据取亚组就可以了,下面来演示一下,先把数据分成两个亚组
bc.m<-subset(bc,bc$fage==0)
bc.f<-subset(bc,bc$fage==1)
分好亚组后就可以得到两个数据,我们就可以像之前一样,分别对每个单组一样跑循环,最后生成两个数据dat1,dat2
我们需要把这两个图合并,生成一个绘图数据
datapolt<-rbind(dat1,dat2)
最后绘图
pd <- position_dodge(0.1)ggplot(datapolt, aes(x=lag, y=ER, colour=group)) + geom_errorbar(aes(ymin= ERlp, ymax= ERup), width=.1) +geom_line(position=pd) +geom_point(position=pd)
修饰一下
ggplot(datapolt, aes(x=lag, y=ER, colour=group)) + geom_errorbar(aes(ymin= ERlp, ymax= ERup), width=.1) +geom_line(position=pd) +geom_point(position=pd)+theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+xlab("lag天数")+ylab("ER")+ggtitle("lag天数与ER关系")
如果不想要连线
ggplot(datapolt, aes(x=lag, y=ER, colour=group)) + geom_errorbar(aes(ymin= ERlp, ymax= ERup), width=.1) +geom_point(position=pd)+theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+xlab("lag天数")+ylab("ER")+ggtitle("lag天数与ER关系")
这个分类变量是我自己生成的,绘制出图形肯定有点怪,但是方法就是这样了,如果想对多个指标进行分面,可以参照图一方法,我这里就不弄了,下节继续介绍下图的绘制
OK,本章结束觉得有用的话多多分享哟。
原创不易,需要文章数据和全部代码的朋友,请把本文章转发朋友圈集10个赞,截图发给我,嫌麻烦的给我打赏5元截图发给我也可以。
相关文章:
R语言广义可加模型在空气环境污染方面的应用(1)
粉丝私信我希望复制一篇文章的图片,图片来源于文章:Wu C, Yan Y, Chen X, Gong J, Guo Y, Zhao Y, Yang N, Dai J, Zhang F, Xiang H. Short-term exposure to ambient air pollution and type 2 diabetes mortality: A population-based time series st…...
CSDN 编程竞赛二十九期题解
竞赛总览 CSDN 编程竞赛二十九期:比赛详情 (csdn.net) 竞赛题解 题目1、订班服 小A班级订班服了!可是小A是个小糊涂鬼,整错了好多人的衣服的大小。小A只能自己掏钱包来补钱了。小A想知道自己至少需要买多少件衣服。 #include <cstdio…...
基于STM32采用CS创世 SD NAND(贴片SD卡)完成FATFS文件系统移植与测试
一、前言 在STM32项目开发中,经常会用到存储芯片存储数据。 比如:关机时保存机器运行过程中的状态数据,上电再从存储芯片里读取数据恢复;在存储芯片里也会存放很多资源文件。比如,开机音乐,界面上的菜单图…...
K_A12_007 基于STM32等单片机驱动AS608光学指纹识别模块 OLED0.96显示
K_A12_007 基于STM32等单片机驱动AS608光学指纹识别模块 OLED0.96显示一、资源说明二、基本参数参数引脚说明三、驱动说明对应程序:四、部分代码说明1、接线引脚定义1.1、STC89C52RCAS608光学指纹模块1.2、STM32F103C8T6AS608光学指纹模块五、基础知识学习与相关资料下载六、视…...
map和set介绍及其底层模拟实现
致努力前行的人: 要努力,但不要着急,繁花锦簇,硕果累累都需要过程! 目录 1.关联式容器 2.键值对 3.树形结构的关联式容器 3.1set的介绍 3.2set的使用 3.3multiset的使用 3.4map的使用 3.5multimap的使用 4.常见的面试题…...
实现一个比ant功能更丰富的Modal组件
普通的modal组件如下: 我们写的modal额外支持,后面没有蒙版,并且Modal框能够拖拽 还支持渲染在文档流里,上面的都是fixed布局,我们这个正常渲染到文档下面: render部分 <RenderDialog{...restState}visi…...
2023美赛F题思路数据代码分享
文章目录赛题思路2023年美国大学生数学建模竞赛选题&论文一、关于选题二、关于论文格式三、关于论文提交四、论文提交流程注意不要手滑美赛F题思路数据代码【最新】赛题思路 (赛题出来以后第一时间在CSDN分享) 最新进度在文章最下方卡片,加入获取一手资源 202…...
Flutter如何与Native(Android)进行交互
前言 上一篇文章《Flutter混合开发:Android中如何启动Flutter》中我们介绍了如何在Native(Android项目)中启动Flutter,展示Flutter页面。但是在开发过程中,很多时候并不是简单的展示一个页面即可,还会涉及…...
数据库主从复制和读写分离
主从数据库和数据库集群的一些问题 数据库集群和主从数据库最本质的区别,其实也就是data-sharing和nothing-sharing的区别。集群是共享存储的。主从复制中没有任何共享。每台机器都是独立且完整的系统。 什么是主从复制? 主从复制,是用来建立一个和主数…...
Java并发编程面试题——线程安全(原子性、可见性、有序性)
文章目录一、原子性高频问题1.1 Java中如何实现线程安全?1.2 CAS底层实现1.3 CAS的常见问题1.4 四种引用类型 ThreadLocal的问题?二、可见性高频问题2.1 Java的内存模型2.2 保证可见性的方式2.3 volatile修饰引用数据类型2.4 有了MESI协议,为啥还有vol…...
DialogFragment内存泄露问题能不能一次性改好
孽缘 自DialogFragment在Android3.0之后作为一种特殊的Fragment引入,官方建议使用DialogFragment代替Dialog或者AllertDialog来实现弹框的功能,因为它可以更好的管理Dialog的生命周期以及可以更好复用。 然而建议虽好,实用须谨慎,…...
java学习--多线程
多线程 了解多线程 多线程是指从软件或者硬件上实现多个线程并发执行的技术。 具有多线程能力的计算机因有硬件支持而能够在同一时间执行多个线程,提升性能。 并发和并行 并行:在同一时刻,有多个指令在CPU上同时执行并发࿱…...
90后阿里P7技术专家晒出工资单:狠补了这个,真香...
最近一哥们跟我聊天装逼,说他最近从阿里跳槽了,我问他跳出来拿了多少?哥们表示很得意,说跳槽到新公司一个月后发了工资,月入5万多,表示很满足!这样的高薪资着实让人羡慕,我猜这是税后…...
2023美赛C题:Wordle筛选算法
Wordle 规则介绍 Wordle 每天会更新一个5个字母的单词,在6次尝试中猜出单词就算成功。每个猜测必须是一个有效的单词(不能是不能组成单词的字母排列)。 每次猜测后,字母块的颜色会改变,颜色含义如下: 程…...
SpringBoot 集成 Kafka
SpringBoot 集成 Kafka1 安装 Kafka2 创建 Topic3 Java 创建 Topic4 SpringBoot 项目4.1 pom.xml4.2 application.yml4.3 KafkaApplication.java4.4 CustomizePartitioner.java4.5 KafkaInitialConfig.java4.6 SendMessageController.java5 测试1 安装 Kafka Docker 安装 Kafk…...
OpenCV 图像金字塔算子
本文是OpenCV图像视觉入门之路的第14篇文章,本文详细的介绍了图像金字塔算子的各种操作,例如:高斯金字塔算子 、拉普拉斯金字塔算子等操作。 高斯金字塔中的较高级别(低分辨率)是通过先用高斯核对图像进行卷积再删除偶…...
【自学Linux】Linux一切皆文件
Linux一切皆文件 Linux一切皆文件教程 Linux 中所有内容都是以文件的形式保存和管理的,即一切皆文件,普通文件是文件,目录是文件,硬件设备(键盘、监视器、硬盘、打印机)是文件,就连套接字&…...
CUDA C++扩展的详细描述
CUDA C扩展的详细描述 文章目录CUDA C扩展的详细描述CUDA函数执行空间说明符B.1.1 \_\_global\_\_B.1.2 \_\_device\_\_B.1.3 \_\_host\_\_B.1.4 Undefined behaviorB.1.5 __noinline__ and __forceinline__B.2 Variable Memory Space SpecifiersB.2.1 \_\_device\_\_B.2.2. \_…...
为什么重写equals必须重写hashCode
关于这个问题,看了网上很多答案,感觉都参差不齐,没有答到要点,这次就记录一下! 首先我们为什么要重写equals?这个方法是用来干嘛的? public boolean equals (Object object&#x…...
< 每日小技巧:N个很棒的 Vue 开发技巧, 持续记录ing >
每日小技巧:6 个很棒的 Vue 开发技巧👉 ① Watch 妙用> watch的高级使用> 一个监听器触发多个方法> watch 监听多个变量👉 ② 自定义事件 $emit() 和 事件参数 $event👉 ③ 监听组件生命周期常规写法hook写法ὄ…...
数据结构与算法之二分查找分而治之思想
决定我们成为什么样人的,不是我们的能力,而是我们的选择。——《哈利波特与密室》二分查找是查找算法里面是很优秀的一个算法,特别是在有序的数组中,这种算法思想体现的淋漓尽致。一.题目描述及其要求请实现无重复数字的升序数组的…...
训练自己的中文word2vec(词向量)--skip-gram方法
训练自己的中文word2vec(词向量)–skip-gram方法 什么是词向量 将单词映射/嵌入(Embedding)到一个新的空间,形成词向量,以此来表示词的语义信息,在这个新的空间中,语义相同的单…...
ubuntu系统环境配置和常用软件安装
系统环境 修改文件夹名称为英文 参考链接 export LANGen_US xdg-user-dirs-gtk-update 常用软件安装 常用工具 ping 和ifconfig工具 sudo apt install -y net-tools inetutils-ping 截图软件 sudo apt install -y net-tools inetutils-ping flameshot 录屏 sudo apt-get i…...
【1139. 最大的以 1 为边界的正方形】
来源:力扣(LeetCode) 描述: 给你一个由若干 0 和 1 组成的二维网格 grid,请你找出边界全部由 1 组成的最大 正方形 子网格,并返回该子网格中的元素数量。如果不存在,则返回 0。 示例 1&#…...
windows11安装sqlserver2022报错
window11安装SQL Server 2022 报错 糟糕… 无法安装SQL Server (setup.exe)。此 SQL Server安装程序介质不支持此OS的语言,或没有SQL Server英语版本的安装文件。请使用匹配的特定语言SQL Server介质;或安装两个特定语言MUI,然后通过控制面板的区域设置…...
Python快速上手系列--日志模块--详解篇
前言本篇主要说说日志模块,在写自动化测试框架的时候我们就需要用到这个模块了,方便我们快速的定位错误,了解软件的运行情况,更加顺畅的调试程序。为什么要用到日志模块,直接print不就好了!那得写多少print…...
【THREE.JS学习(1)】绘制一个可以旋转、放缩的立方体
学习新技能,做一下笔记。在使用ThreeJS的时候,首先创建一个场景const scene new THREE.Scene();接着,创建一个相机其中,THREE.PerspectiveCamera()四个参数分别为:1.fov 相机视锥体竖直方向视野…...
数仓实战 - 滴滴出行
项目大致流程: 1、项目业务背景 1.1 目的 本案例将某出行打车的日志数据来进行数据分析,例如:我们需要统计某一天订单量是多少、预约订单与非预约订单的占比是多少、不同时段订单占比等 数据海量 – 大数据 hive比MySQL慢很多 1.2 项目架…...
python虚拟环境与环境变量
一、环境变量 1.环境变量 在命令行下,使用可执行文件,需要来到可执行文件的路径下执行 如果在任意路径下执行可执行文件,能够有响应,就需要在环境变量配置 2.设置环境变量 用户变量:当前用户登录到系统,…...
BeautifulSoup文档4-详细方法 | 用什么方法对文档树进行搜索?
4-详细方法 | 用什么方法对文档树进行搜索?1 过滤器1.1 字符串1.2 正则表达式1.3 列表1.4 True1.5 可以自定义方法2 find_all()2.1 参数原型2.2 name参数2.3 keyword 参数2.4 string 参数2.5 limit 参数2.6 recursive 参数3 find()4 find_parents()和find_parent()5…...
外贸都是在哪些网站做/seo排名优化方式
在《Server 层混杂信息字典表 | 全方位认识 information_schema》中,我们详细介绍了information_schema下的状态变量、系统变量、进程状态、字符集和校对规则等字典表,本期我们将为大家带来系列第五篇《InnoDB 层系统字典表 | 全方位认识 information_sc…...
qq空间可以做网站吗/电商seo优化
这里是对docker compose 网络配置的一些说明,详细的文档参考: https://docs.docker.com/compose/networking/ 1 default network 如果不显式指定,Compose会为每一个app设置一个default网络。每个service的container会加入这个default网络并…...
张家港网站哪家做的好/电商网站设计论文
一、ACL介绍 权限控制(ACL)主要为RocketMQ提供Topic资源级别的用户访问控制。 用户在使用RocketMQ权限控制时,可以在Client客户端通过 RPCHook注入AccessKey和SecretKey签名;同时,将对应的权限控制属性(包…...
怎么把自己做的网页上传网站/百度快照优化排名
所遇问题:之前测试的项目中对“长图片”进行分享时,发现部分安卓手机无法调起微信……紧接着就被通知测试人员多测几部机型,找下规律……what……那么多手机难道我们要一部部都过?如果真的把所有手机过一部部测,不紧耗…...
成都市四方建设工程监理有限公司网站/怎么进行seo
翻译:疯狂的技术宅 原文:https://www.edureka.co/blog/interview-questions/react-interview-questions https://mp.weixin.qq.com/s/QxvUGi6li8sqQJFbFtgBpw 29. 你对受控组件和非受控组件了解多少? 受控组件非受控组件1. 没有维持自己的状…...
b站推广网站入口mmm/收录提交入口
题目题目描述贝茜在谷仓外的农场上,她想回到谷仓,在第二天早晨农夫约翰叫她起来挤奶之前尽可能多地睡上一觉.由于需要睡个好觉,贝茜必须尽快回到谷仓.农夫约翰的农场上有N(2≤N≤1000)个路标,每一个路标都有…...