基于PM2.5浓度值数据集的相关分析
邹潇
摘要
本文通过对2010年到2014年北京市pm2.5数据进行相应的数据分析，首先将导入数据进行数据清洗，得到本论文想要数据集。其次，探究影响细微颗粒物pm2.5浓度值的关键因素，详细分析pm2.5与影响因素之间的相关性，同时又考虑到本数据是时间序列，前一小时会对后一小时数据集产生较大影响，所以在最后的lasso回归模型中我们加入了前期对后期影响，作为因变量，并做变量选择，建立回归模型。最后，通过对该数据分析得到一些pm2.5浓度值变化的规律以及可以实现前一小时相关指标来预测后一小时pm2.5浓度值。1. 绪论
1.1 研究背景
近年来，空气污染问题已经严重影响地球的气候系统，更加影响人们生活和工作，所以城市的空气污染问题越来越得到了人们的广泛关注，雾霾更是政府和人们关心的热点话题，其严重影响人们的身体健康和正常生活。 随着雾霜天气越来越大范围和频繁的发生，对人们的正常生活生产活动以及身体健康造成了不良影响。尤其在每年冬季全国大部分城市都是雾霾频发髙峰期，据相关报道称在雾霾髙发期，同时持续时间长的时候，各大医院呼吸道不适等症状的患者人数呈增加趋势，主要有３类疾病患者增多分别是呼吸道疾病患者、也脑血管疾病患者、皮炎等皮肤过敏患者。当在重雾疆天气下，空气中的pm2.5浓度会很高，由于与较粗的大气颗粒物相比，pm2.5粒径小，面积大，活性强，易附带有毒、有害物质（例如，重金属、微生物等），且在大气中的停留时间长、输送距离远，因而对人体健康和大气环境质量的影响更大。细微颗粒物pm2.5因其体积小，渗透能力强，而能最终进入到细支气管壁，影响肺中的气体交换，扰乱身体其他器官功能，对人体造成极大的伤害。 雾霾不仅对人体健康影响极大，并且对我们生活出行也造成了不少的烦恼，就在2015年十一二月，全国范围内尤其是京津冀周边等地区出现严重的雾霾天气，造成大多数中小学停课，机动车单双限行，所以在这样的情况下，能有效知道空气质量状况，并做好提前预警工作显得极为重要。
1.2 研究意义
对于上述问题的研究背景可以看出本文主要研究目的就是针对空气质量问题利用历史数据探讨影响其影响空气质量因素有哪些的因素有哪些，并建立预测模型。因为细微颗粒物pm2.5作为衡量空气质量的重要指标之一，对这项影响到人们正常生活的污染物浓度做出相关的分析，并提供未来某段时间内pm2.5的相关信息，当预测空气污染指标超标时，以使得人们可以减少室外暴露和出行，在一定程度上可以减少对人体的伤害，同时对那些在室外工作运动的人群，通过对污染指标的参考可以大大提高工作和生活的效率，故本文对pm2.5的影响以及大体趋势分析是具有研究意义的。
2. 数据来源及导入
本文数据来源于UCI网站上的Beijing PM2.5 Data Set，数据包括2010年1月1日至2014年12月31日间北京pm2.5指数以及相关天气指数数据，可查看其每个字段长度均相等，无肉眼缺失。rm(list = ls())
library(tidyverse)
-- Attaching packages ------------------ tidyverse 1.2.1 --
√ ggplot2 3.1.0 √ purrr 0.2.5
√ tibble 1.4.2 √ dplyr 0.7.8
√ tidyr 0.8.2 √ stringr 1.3.1
√ readr 1.2.1 √ forcats 0.3.0
-- Conflicts --------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
pm25_train=read.csv("D:/code/workplace/pm25_train.csv")
pm25_train%>%map(length)%>%str()
List of 13
$ date : int 35746
$ hour : int 35746
$ pm2.5 : int 35746
$ DEWP : int 35746
$ TEMP : int 35746
$ PRES : int 35746
$ Iws : int 35746
$ Is : int 35746
$ Ir : int 35746
$ cbwd_NE: int 35746
$ cbwd_NW: int 35746
$ cbwd_SE: int 35746
$ cbwd_cv: int 35746其中训练数据主要包括35746条记录，13个字段。其主要解释如下：
（1）date：观测数据发生的日期（年-月-日）
（2）hour：观测数据发生的时间点（时）
（3）pm2.5：观测时间点对应的pm2.5指数（(ug/m^3)
（4）DEWP：露点，空气中水气含量达到饱和的气温（℃）
（5）TEMP：温度，观测时间点对应的温度（℃）
（6）PRES：压强，观测时间点对应的压强（hPa）
（7）Iws：累积风速，观测时间点对应的累积风速（m/s）
（8）Is：累计降雪，到观测时间点为止累计降雪的时长（小时）
（9）Ir：累计降雨，到观测时间点为止累计降雨的时长（小时）
（10）cbwd_NE：观测时间点对应的风向为东北风(m/s)
（11）cbwd_NW：观测时间点对应的风向为西北风(m/s)
（12）cbwd_SE：观测时间点对应的风向为东南风(m/s)
（13）cbwd_cv：观测时间点对应的风向为静风(m/s) 3.数据整理
3.1 数据分拆与合并
library(lubridate)
Attaching package: 'lubridate'
The following object is masked from 'package:base':
date
pm25_train=pm25_train%>%
separate(date,into=c("year","month","day"),sep="-")
pm25_train$hour=sprintf("%02d",pm25_train$hour)
temp=pm25_train%>%
unite(id,year,month,day,hour,sep="")
pm25_train=pm25_train%>%
mutate(id=temp$id)%>%
select(id,everything())首先是把时间数据分拆成，方便后续根据年月日来分析数据特征；再通过对数据的年月日以及小时进行合并得到每个数据id，方便后续分析 通过对原始数据进行分拆和合并,使得数据更加利于分析。3.2 处理隐含缺失
pm25_new=pm25_train%>%
add_count(year,month,day)%>%
filter(n==24)
length(pm25_new$id)
[1] 32328从原始数据可以看到，数据集是按小时为单位进行统计的，所以每天应有24条数据，但有部分数据因为观测原因某些时点并未测得数据，为了方便以天为单位统计pm2.5的值，所以某一日中出现了缺失值，则删除该日的所有时刻的pm2.5的数据，故应剔除这部分数据，从结果中可以看到，每天不足24条数据的有3418条，剔除后剩余32238条记录。
按年份整理的数据集如下：pm25_new=pm25_new%>%
select(id,everything())
pm25_new=pm25_new%>%
mutate(date=make_date(year,month,day))
pm25_new%>%
group_by(year)%>%
summarize(count=n())%>%
knitr::kable(align = 'c',digits = 5)| year | count |
|---|---|
| 2010 | 6528 |
| 2011 | 5952 |
| 2012 | 6120 |
| 2013 | 6816 |
| 2014 | 6912 |
4. 数据探索
4.1数据可视化
我们通过对每天数据取均值，统计每天pm2.5变化情况，并按年份来如下图所示：library(lubridate)
library(grid)
daily=pm25_new%>%
mutate(date=make_date(year,month,day))%>%
group_by(date)%>%
summarize(mean=mean(pm2.5))
p1=daily%>%
ggplot(aes(date,mean))+
geom_line()
p2=daily%>%
ggplot(aes(x=year(date),y=mean,group=year(date)))+
geom_boxplot()+
guides(fill=FALSE)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2,1)))
vplayout <- function(x,y){
viewport(layout.pos.row = x, layout.pos.col = y)
}
print(p1, vp = vplayout(1,1))
print(p2, vp = vplayout(2,1)) 
可以看出pm2.5浓度的季节性变化非常明显。通常来说冬季pm2.5浓度最高，夏季pm2.5浓度最低。这也符合常识，冬季pm2.5浓度高的原因主要在于燃煤供暖，生物质燃烧和不利于空气扩散的气候；而夏季pm2.5浓度明显降低，这主要与人为排放和以供暖为目的的生物质的燃烧减少有关，当然海洋季风气候带来的洁净空气也对空气污染起着较强的稀释作用。
同时我们将pm2.5值按年度做箱线图看整体分布情况，可以看到近几年pm2.5浓度值稍有下降，箱线图上各分位数的值都是有所降低的，也均存在异常值，并且异常值的大小和数量并未有明显减少。为了更明显发现统计规律，根据2016年1月1日实施的国家标准GB3095-2012《环境空气质量标准》中对PM2.5等级的划分，我们可以对数据做以下划分：class=cut(daily$mean,breaks=
c(0,35,75,115,150,250,max(daily$mean)),
labels = c("优","良","轻度污染","中度污染","重度污染","严重污染"))
daily=daily%>%
mutate(class=class)
daily%>%
ggplot(aes(x=reorder(class,mean,fun=median),y=mean))+
geom_boxplot(aes(color=class))
将数据按照空气质量等级和颜色进行划分后，绘制分层柱状图，如下：daily%>%
ggplot(aes(x=year(date),fill=class),position="fill")+
geom_bar()
可以看出，从2010年至2014年，北京市pm2.5等级为优良的天数相比2010年都有所提升，空气质量也呈现逐年改善的趋势，但同时我们也注意到严重污染天数持续几年并未出现减少趋势，这也表明我们雾霾治理主要是针对轻中度污染天气效果明显，对严重污染天气效果不佳。
为了进一步探究pm2.5指数的季节性特征，我们需要对原数据进行季节划分，其中，三、四、五月划为春天，六、七、八月划为夏天，九、十、十一月划为秋天，十二、一、二月为冬天。我们的目标是将一年pm2.5浓度随季节的变化的特征、趋势提取出来加以分析。其具体划分如下：month=month(daily$date)
season=ifelse(
month>=3 & month<=5,"spring",
ifelse(month>=6 & month<=8,"summer",
ifelse(month>=9 & month<=11,"fall","winter")))
daily=daily%>%
mutate(season=season,month=month)
p3=daily%>%
ggplot(aes(x=reorder(season,mean,fun=median),y=mean))+
geom_boxplot(aes(color=season))+
theme(legend.title = element_blank(), legend.position = "top")
p4=daily%>%
ggplot(aes(x=mean))+
geom_freqpoly(aes(color=season),binwidth=30)+
theme(legend.title = element_blank(), legend.position = "top")
grid.newpage()
pushViewport(viewport(layout = grid.layout(1,2)))
print(p3, vp = vplayout(1,1))
print(p4, vp = vplayout(1,2)) 
首先我们可以看到，pm2.5浓度值趋势性极为明显，在冬季北京的pm2.5浓度不管是均值还是中位数均比其他季节要高，并且异常值点更多，浓度更高。在频数分布图中我们可以看到，pm2.5浓度值得变化在四季走势大致相同，但当出现pm2.5值超标的天气，在冬季出现的天数明显大于其他季节，同时可以看到夏季pm2.5整体在四季中整体较低，也符合我们生活常识。当然除了季节性，pm2.5浓度日昼夜变化幅度也较为明显，其具体情况如下图所示：
p5=pm25_new%>%
group_by(hour)%>%
summarize(pm25mean_hour=mean(pm2.5))%>%
ggplot(aes(x=hour,y=pm25mean_hour))+
geom_point()
p6=pm25_new%>%
ggplot(aes(x=hour,y=pm2.5,fill=hour))+
geom_boxplot()+
guides(fill=FALSE)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2,1)))
print(p5, vp = vplayout(1,1))
print(p6, vp = vplayout(2,1)) 
可以清晰看出，北京pm2.5浓度昼夜变化幅度较大，一天中pm2.5浓度最低的时候是在下午三点和四点之间，查阅相关文献可以知道，主要是因为这个时候大气边界层升高，风速也增加，但从下午4点后，pm2.5的浓度开始增加，大气边界层的高度不断降低，也逐渐到了下班时间，车辆排放以及二氧化氮的排放明显增加。而在凌晨1点左右，由于车辆量减少，pm2.5浓度逐渐降低，直到下午4点，如此循环。4.2 变量相关性分析
library(corrgram)
library(gclus)
Loading required package: cluster
corrgram(pm25_new[6:16],order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie,text.panel=panel.txt,
main="Correlogram of pm25_new intercorrelations")
通过对数据做相关关系图，可以看到与pm2.5线性相关关系均不是很强，但总体来说露点（DEWP）和累计风速（lws）对pm2.5有较强的线性关系，下面我们进一步来分析其相关性。standard=function(df){
(df-min(df))/(max(df)-min(df))
}
p7=pm25_new%>%
ggplot(aes(x=date,y=standard(Iws)))+
geom_line()+
labs(y="Iws")
p8=pm25_new%>%
ggplot(aes(x=date,y=standard(pm2.5)))+
geom_line()+
labs(y="pm2.5")
grid.newpage()
pushViewport(viewport(layout = grid.layout(2,1)))
print(p7, vp = vplayout(1,1))
print(p8, vp = vplayout(2,1))
通过对标准化的数据做线形图，可以看出整体趋势Iws与pm2.5有显著负相关，当累计风速增大的时候，必然使得pm2.5浓度有所下降；同理可以看到露点呈现周期性特征，与pm2.5浓度呈现正向相关。5. 数据建模
5.1 季节效应分析
从以上的分析中我们可以看出，pm2.5有强烈的季节效应，我们可通过对月份进行建模来试着消除这一效应
library(modelr)
mod=lm(mean~month,data=daily)
daily=daily%>%
add_residuals(mod)
daily%>%
ggplot(aes(date,resid))+
geom_ref_line(h=0)+
geom_line(color="grey50")+
geom_smooth()
`geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
从残差变化来看，通过对月份进行建模对季节性剔除效应还是比较明显，表示了每天的pm2.5浓度值与给定一月中每天预期pm2.5浓度值的偏差,可以看出去除大部分月份内的效应，我们可以从一年内的残差图可以看得更加清晰。daily%>%
filter(year(date)=="2014")%>%
ggplot(aes(date,resid))+
geom_ref_line(h=0)+
geom_line(color="grey50")+
geom_smooth()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
以2014年pm2.5数据为例，可以明显看出数据中还是存在大量的离群点。所以拟合的平均趋势与典型值之间差别加大，为了改善这一问题，我们采用稳健回归模型，同时我们加入季节因素来做回归看二者共同对pm2.5的影响。mod2=MASS::rlm(mean~month*season,data = daily)
daily%>%
add_residuals(mod2,"resid")%>%
filter(year(date)=="2014")%>%
ggplot(aes(date,resid))+
geom_ref_line(h=0)+
geom_line(color="grey50")+
geom_smooth()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
还是以2014年为例，总体来说更容易看出长期趋势，对异常值也看得更加清晰。5.2 LASSO回归
5.2.1 理论介绍
LASSO(The Least Absolute Shrinkage and Selectionator operator)是Robert Tibshirani 在1996年提出的特征选择算法，并能够准确且稳定地进行特征选择。通过构造一个惩罚函数来获得较为精炼的模型，最终使得一些变量的系数为零。其基本思想是在回归系数的绝对值之和小于某个指定常数的条件下，使得残差平方和最小化，从而使一些回归系数为0，从而达选择变量的目的。
lasso回归的特色就是在建立广义线型模型的时候，这里广义线型模型包含一维连续因变量、多维连续因变量、非负次数因变量、二元离散因变量、多元离散因变，除此之外，无论因变量是连续的还是离散的，lasso都能处理，总的来说，lasso对于数据的要求是极其低的，所以应用程度较广；除此之外，lasso还能够对变量进行筛选和对模型的复杂程度进行降低。这里的变量筛选是指不把所有的变量都放入模型中进行拟合，而是有选择的把变量放入模型从而得到更好的性能参数。更多的变量在拟合时往往可以给出一个看似更好的模型，但是同时也面临过度拟合的危险。
而对于本问题而言，为了比较不同变量对pm2.5浓度影响的大小，我们需要先对数据进行标准化，这里我们采用最常用的z-score标准化方法，即将数据转化为均值为0，方差为1的标准化数据。standard=function(df){
(df-mean(df))/sd(df)
}
pm25_standard=standard(as.matrix(pm25_new[6:16]))
pm25_standard=data.frame(pm25_standard)同时根据前面对一天各小时的pm2.5浓度值得分析，可以了解到pm2.5值与前一小时pm2.5浓度值关系十分密切，故把前一小时pm2.5值和各项指标一起作为自变量来预测后一小时pm2.5值，并用前9/10作为训练集，后1/10作为测试集，进行pm2.5值的预测。#divide train and test
N1=round(length(pm25_standard$pm2.5)*9/10)
pm25_data_before=pm25_standard[1:N1-1,]
pm25_after=pm25_standard$pm2.5[2:N1]
pm25_train1=data.frame(pm25_after,pm25_data_before)下面通过lasso对标准化的数据做变量选择library(glmnet)
Loading required package: Matrix
Attaching package: 'Matrix'
The following object is masked from 'package:tidyr':
expand
Loading required package: foreach
Attaching package: 'foreach'
The following objects are masked from 'package:purrr':
accumulate, when
Loaded glmnet 2.0-16
x=as.matrix(pm25_train1[2:11])
y=as.matrix(pm25_train1[1])
cvfit_lasso = cv.glmnet(x,y)
plot(cvfit_lasso)
cvfit_lasso$lambda.min
#> [1] 0.0009374997从结果我们可以看到最小的lambda=0.00094，进行如下变量选择，模型拟合优度R^2达到0.92，表明模型拟合效果较好。coefs_lasso = coef(cvfit_lasso, s = "lambda.min")
as.matrix(coefs_lasso)%>%
knitr::kable(caption="beta",align = 'c',digits = 5)| 1 | |
|---|---|
| (Intercept) | -4.48010 |
| pm2.5 | 0.94279 |
| DEWP | 0.00000 |
| TEMP | -0.11005 |
| PRES | 0.00000 |
| Iws | 0.00000 |
| Is | 0.00000 |
| Ir | -0.14543 |
| cbwd_NE | -2.99954 |
| cbwd_NW | -9.22538 |
| cbwd_SE | 0.00000 |
从lasso结果来看，系数为0的变量有：TEMP、PRES、Iws、Is、cbwd_cv，结果剩余5个变量。通过回归系数我们可以看出有无西北风对pm2.5值影响较大，都知道我们国家冬季大部分地区，包括北京冬季都是西北风，而刮风对于pm2.5浓度有很好的稀释作用，所以对其值的大小有很大影响。同时前一小时pm2.5值对后一小时pm2.5浓度也有较大影响。
并且从上面回归系数可以看出来，温度（TEMP）、累计降雨(Is)均与我们pm2.5指数呈现负相关，这也和我们做的季节性分析结论相一致，在夏季温度较高，降雨也较多，pm2.5浓度是四季中最低的，同理知道pm2.5在冬季的浓度最高。N1=round(length(pm25_standard$pm2.5)*9/10)
N2=length(pm25_standard$pm2.5)
pm25_data_after=pm25_standard[N1:(N2-1),]
pm25_true=pm25_standard$pm2.5[(N1+1):N2]
pm25_test_data=data.frame(pm25_data_after,pm25_true)
x1=as.matrix(pm25_train1[2:11])
pm25_pred=predict(cvfit_lasso, newx = x1[1:(N2-N1),], s = "lambda.min")
pm25_pred=pm25_pred*sd(pm25_new$pm2.5)+mean(pm25_new$pm2.5)
pm25_true=pm25_true*sd(pm25_new$pm2.5)+mean(pm25_new$pm2.5)
pm25_contrast=data.frame(pm25_pred,pm25_true)
pm25_contrast%>%
ggplot(aes(x=1:3233,y=pm25_true))+
geom_point()+
geom_line(aes(x=1:3233,y=pm25_pred),color="red")+
labs(title=paste("contrast the predict and true value"),
x="id",y="pm25_pred/pm25_true")
从上图可以看出，lasso模型预测结果大概能够反映pm2.5变化趋势且拟合效果较好，可以根据上述模型提前一小时知道pm2.5具体变化情况，判断其升高和降低情况。6.总结与不足
6.1 总结
本文主要通过对北京市2010至2014年pm2.5数据进行统计分析，并提供一步预测的lasso模型来对pm2.5浓度值提前预测。其主要分析结论如下：
从主要数据本身出发，做了相关的可视化分析，探究影响pm2.5浓度值变化趋势，分析pm2.5浓度值与其他影响因素之间的关系，通过对实验结果进行总结分析我们可以知道：从收集数据的年份来看，pm2.5浓度值有逐年有所下降，中度污染天数明显减少，以全年数据来看pm2.5值有很强的季节性特征，具体表现为冬季pm2.5浓度值平均最高，夏季最低；从一天数据分析知道，中下午3点到4点间pm2.5浓度值最低，而在凌晨1点左右达到一天中的顶峰其后有所缓和。
其次，分析了pm2.5浓度值与其他因素之间的相关性，得到与pm2.5影响因素较大的变量，可以发现累计风速与pm2.5相关性较强，并且呈明显负相关，而空气中水汽含量达到饱和的气温即露点与我们pm2.5呈现较强的正相关性。
最后，建立了基于前一小时pm2.5相关数据来预测下一小时pm2.5浓度值的lasso回归模型，从pm2.5浓度值预测效果来看，所采用的预测模型可以取得较好的预测效果。6.2 不足
本文的分析中，虽然对pm2.5浓度值做了定量的分析，但总体来说还是有一些不足之处。首先，样本选取是北京市2010至2014年数据，本身数据集比较单一，采集时间也偏早了，不具有时效性；同时，本文选择数据没有考虑其周围区域相关因素的影响，由于细微颗粒物易受风的影响，随着空气气流的传输，周围的pm2.5浓度也会影响目标区域pm2.5的变化；最后，本文采用的lasso回归模型，只能针对前一小时对后一小时的预测，所以有必要扩大时间区间，做更全面的分析。