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<-NULL
  pm10lag<-runMean(bc$pm10,0:i)###产生 pm10lag01变量，赋值为过去两天大气PM10浓度的移动均值
  bc$pm10lag<-pm10lag
  fit<-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 ，死亡的超额危险度ER
  ERlp<-(exp((b-1.96*se)* 10)-1)*100 #计算95%CI
  ERup<-(exp((b+1.96*se)* 10)-1)*100 #计算95%CI
  lag<- 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元截图发给我也可以。