R 利用回归分析与时间序列预测北京市PM2.5

注：代码全部在最后,数据来源UCI,链接如下：

https://archive.ics.uci.edu/ml/datasets/Beijing+PM2.5+Data

摘要

现代社会科技进步，人们的生活质量逐步提高，但伴随着各类工业和科技的发展，环境问题凸显，最初人们粗放式的经济发展方式在一定程度上对环境造成不可逆转的破坏。在各种环境污染问题，空气污染问题又是如今人们关注的重中之重。北京是我国首都，同时也是我国空气污染较为严重的几个北方城市之一，因此关注北京市空气污染情况是我国观测空气污染情况的重要关注对象之一。在本文中，首先，我就北京市PM2.5及其他污染物和气象情况的近几年观测数据做出相应回归分析，建立回归模型试图利用已知的其他污染物和气象情况信息对PM2.5做出预测；其次，利用现如今已观察到的数据做单变量时间序列分析，寻找其变化规律、做简单预测，为北京市人民规划出行做好提前空气质量预警。

关键词：PM2.5 多元回归分析 时间序列分析

 

一、数据集介绍

该数据集来自于北京市环境检测中心，包括了2013年3月1日到2017年2月28日，12个国家控制的空气质量监测站每小时的空气污染数据，且每个空气质量站点的气象数据都与中国气象局站相匹配。此空气污染数据集包括12个站点数据集，每个数据集包括了35064个时间数据，共有（35064*12）个样本。

变量包括：PM2.5、PM10、SO2（二氧化硫）、NO2（二氧化氮）、CO（一氧化碳）、O3（臭氧）为空气污染物；剩余变量TEMP（温度）、PRES（大气压）、DEWP（露点温度）、RAIN（降雨量）、WD（风向）、WSPM（风速）天气情况。

站点包括：万寿西宫（西城区）、官园（西城区）、万柳（海淀区）、天坛（东城区）、农展馆（朝阳区）、奥体中心（朝阳区）、怀柔、古城（顺义区）、顺义、东寺（平谷区）、定陵（昌平区）、昌平共十二个，其地理位置分布如下：

 

二、研究问题一

（一）研究问题

将SO2、NO2、CO、O3、TEMP（温度）、PRES（大气压）、DEWP（露点温度）、RAIN（降雨量）、WD（风向）、WSPM（风速）作为解释变量，Y作为被解释变量，以预测误差最小为原则，选择最优多元回归模型。

（二）研究方法

多元回归

（三）研究步骤

1.数据集选择

由于该数据集有北京12个观测点的完整时间数据，以下以PM2.5为例绘制出PM2.5分布的盒图和小提琴图，并按照均值从从低到高排列，发现12个观测点PM2.5的分布基本相似，以下选择其中一个地区奥体中心的数据集进行分析。

2.数据可视化

    在进行定量建模之前，先绘制图像对变量分布和变量之间的相关关系进行直观感受与定性分析。

观察SO2、NO2、CO和O3这四种污染物分别与PM2.5的相关关系：

可以看出PM2.5随着CO、NO2浓度的上升而呈线性上升趋势，O3和SO2和PM2.5的线性关系没有那么强烈。

接下来考虑气温、压强、露点温度、降雨量与被解释变量PM2.5的关系：

从图中可以看出气温、压强、露点温度、降雨量和PM2.5的相关性更弱，但具体会对PM2.5的预测有怎样的效果，需要在回归模型中进一步分析。

最后考虑风速和风向与被解释变量PM2.5的关系：

   在第一幅图可以看中，图例从上到下以PM2.5均值大小排序，可以看出西北方向的风PM2.5峰度高且整体分布靠左，东南方向的风峰度低而整体分布靠右，也就是PM2.5值整体偏大；风速和PM2.5的线性相关关系不是很明显。

3.回归分析（以奥体中心为例）

（1）前期准备

① 对待选入模型的变量进行平稳性检验，防止出现伪回归，使用ADF单位根检验，拒绝原假设，不存在单位根，序列平稳。

② 对所选的“奥体中心”数据做删失处理，由于数据量足够多，直接删除带有缺失值的行。

（2）变量选择：

先做全变量回归，发现有大量变量存在系数t检验结果不显著，并且存在VIF值较高的变量，需要进行变量选择，这里尝试了两种变量选择的方法，两种方法筛选出的变量相同。

① 根据BIC准则做逐步回归，依次从全模型剔除变量，直至BIC值最小。

② 利用10折交叉验证，将数据根据被解释变量均匀分为10份，每一份分别当作测试集，计算利用10个训练集训练出的模型在相应测试集上的MSE并取平均，建立循环寻找使得交叉验证所得MSE平均值最小的变量组合，即预测效果最好的变量组合。

（3）回归结果：

Coefficients:                    

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	232.812250	64.033810	3.636	0.000278	***
SO2	0.463988	0.019179	24.192	<2e-16	***
NO2	0.710602	0.015931	44.605	<2e-16	***
CO	0.035909	0.000440	81.608	<2e-16	***
O3	0.249223	0.009159	27.211	<2e-16	***
TEMP	-1.154538	0.086582	-13.335	<2e-16	***
PRES	-0.255064	0.062654	-4.071	4.7e-05	***
DEWP	1.612418	0.061157	26.365	<2e-16	***
RAIN	-1.605167	0.413330	-3.883	0.000103	***
WSPM	5.117095	0.378114	13.533	<2e-16	***
wdENE	-3.069006	1.125061	-2.728	0.006382	**
wdN	3.049390	1.506081	2.025	0.042913	*
wdNNW	4.275978	1.719299	2.487	0.012891	*
wdSE	5.081616	1.756569	2.893	0.003822	**
wdSSE	4.561950	2.014055	2.265	0.023523	*
wdSSW	3.442118	1.493166	2.305	0.021165	*
wdW	-5.690030	1.895453	-3.002	0.002687	**
wdWSW	-3.425869	1.428746	-2.398	0.016505	*

其拟合优度及显著性水平检验结果如下表：

Residual standard error	42.81 on 15891 degrees of freedom
Multiple R-squared	0.7164
Adjusted R-squared	0.7161
F-statistic	2361 on 17 and 15891 DF
p-value	< 2.2e-16

（4）模型检验：

a. 独立性检验：拒绝原假设rho==0，认为随机扰动项存在自相关；

b. 线性性检验：由图可以看出RAIN不是线性的；

 

c.同方差检验：拒绝原假设（方差恒定），模型存在异方差；

 

d.多重共线性检验:

下表计算了变量选择之后建立的回归模型的各解释变量的VIF

SO2	NO2	CO	O3	TEMP	PRES
1.658091	2.978911	2.513574	2.442158	8.321385	3.619437
DEWP	RAIN	WSPM	wdENE	wdN	wdNNW
6.036338	1.027168	1.801353	1.102301	1.064934	1.105588
wdSE	wdSSE	wdSSW	wdW	wdWSW	-
1.047875	1.042624	1.091834	1.033414	1.057227	-

 

e.鉴别强影响点与离群点：由结果得只有10个离群点和4个强影响点，相对于样本观测值数量而言极小，可忽略其影响；