北京建筑大学理学院，北京

收稿日期：2022年7月22日；录用日期：2022年8月15日；发布日期：2022年8月25日

摘要

近年来，随着社会经济发展和人们生活水平的提高，环境污染越来越受到人们的重视。在评估某地区环境好坏时，空气质量是其中一个重要的指标。同时，空气质量也与每个人息息相关。为此，本文通过2020年全国各省市的空气质量数据及各省相应经济数据、人口数据、公共预算支出等数据，使用最小二乘法、逐步回归法、岭回归法、Lasso回归、主成分回归法五种方法构建了多元线性回归模型。利用得到的模型进行分析，结果表明：人均可支配收入、人口数、工业产品产量与空气质量呈负相关；公共预算支出中的科学技术预算支出与空气质量呈正相关。并且工业产品产量及科学技术预算支出对空气质量的影响较为显著。除此之外，本文对比了五种线性回归方法建模精度，得到在存在多重共线性时，逐步回归法与主成分法预测精度较高。其余四种方法都明显好于普通最小二乘法。

关键词

线性回归，空气质量，多重共线性

Application of Linear Regression Methods in the Analysis of Air Quality Influencing Factors

Yihe Liu, Weiyan Mu, Tong Jin

School of Science, Beijing University of Civil Engineering and Architecture, Beijing

Received: Jul. 22^nd, 2022; accepted: Aug. 15^th, 2022; published: Aug. 25^th, 2022

ABSTRACT

In recent years, with the development of social economy and the improvement of people’s living standard, more and more people pay attention to environmental pollution. Air quality is an important index in evaluating the environment of a certain area. At the same time, air quality is closely related to everyone. Therefore, based on the air quality data of all provinces and cities in 2020 and the corresponding economic data, population data and public budget expenditure data of each province, this paper constructed a multiple linear regression model using five methods including least square method, stepwise regression method, mountain regression method, Lasso regression method and principal component regression method. The results show that per capita disposable income, population and industrial product output are negatively correlated with air quality. Science and technology budget expenditure in public budget expenditure is positively correlated with air quality. The output of industrial products and budget expenditure of science and technology have a significant impact on air quality. In addition, this paper compares the modeling accuracy of five linear regression methods, and obtains that the prediction accuracy of stepwise regression method and principal component method is higher in the presence of multicollinearity. The other four methods are obviously better than ordinary least square method.

Keywords:Linear Regression, Air Quality, Multicollinearity

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随着我国经济的不断发展，党的十九大指出，我国社会的主要矛盾转化为人民日益增长的美好生活需要同不平衡不充分的发展之间的矛盾。习近平主席说过：“绿水青山就是金山银山”。因此，从国家与社会、个人层面，都对环境质量越来越重视。而环境保护中的空气质量与每个人都息息相关。为探讨各个指标对空气质量的影响，黄煜宁等人以福州市PDPI作为解释变量，空气综合污染指数作为被解释变量，得出2004~2018年期间，福州环境空气质量与经济增长之间的关系 [1]。金仁浩等人在对各区域空气质量数据进行描述分析的基础上，从整体上分析北京市空气质量与气象、社会经济因素之间的关系 [2]。刘亦文通过选取中国30个重点城市2003~2018年的面板数据，考察碳排放总量和强度双约束政策对城市空气质量的平均因果效应 [3]。基于此，本文通过各省空气质量指标与各省经济指标、工业产品产量等数据建立多元线性回归模型，探讨哪些因素对空气质量有影响。为改善空气质量、保护环境应进行哪些举措提供建议。同时探讨各种回归方法，在处理具有多重共线性的数据时的表现。

本文第一部分介绍了对空气质量影响因素的研究；第二部分对数据进行初步处理与分析；第三部分使用最小二乘法、岭回归、Lasso回归、主成分回归法对数据进行建模，得到相应的模型及参数。第四部分对模型精度进行比较，对各个模型的参数进行解释。通过分析得到的模型及参数，得出环境污染与工业产量、经济发展、政府预算的相关关系。

2. 数据录入及初步处理

原始数据使用全国各省会城市空气质量及达到二级的天数作为因变量Y，其中自变量M代表各省市人均可支配收入；自变量G代表每省市人均生产总值；自变量P为各省人口数；自变量B为各省市当年环境保护预算支出；自变量S为各省当年科技预算支出；自变量I是各省当年部分工业产品产量。工业产品产量合计中，选择了部分对环境污染较高且数值较大的工业产品。上述数据均为2020年数据，来自《中国统计年鉴2021》。原始数据见附表1。

由于自变量的量纲不同，为了使数据具有参考性，首先对数据进行标准化。利用R软件中的scale函数，可以得到原始数据进行中心化标准化后的结果。数据标准化结果见附表2。

为了检验标准化后的数据是否具有多重共线性，使用三种方法。利用cor函数、vif函数，分别计算数据的相关系数、VIF值及条件数。得到结果见表1、表2。其中表1保留四位小数。

表1. 相关系数

表2. 条件数及VIF值

由表1可知，有些自变量间相关系数大于0.5，例如：自变量M与自变量G的相关系数约为0.95。因此，可以认为这些系数之间具有相关性。数据具有多重共线性。由表2可得，条件数小于100，可以认为共线程度较小；但自变量M与自变量G的VIF值均大于10，认为数据具有多重共线性。综合三种方法，认为我们所选数据具有多重共线性。因此，在后文中选择的岭回归法、Lasso回归法、主成分法都可以部分的消除数据多重共线性，使得到的结果更加准确。

接下来检验数据是否服从正态分布。本文使用了密度图与QQ图，直观检验数据是否服从或近似服从正态分布。除此之外，使用Shapiro-Wilk方法，从数值上准确判断数据是否通过正态性检验。密度图与QQ图见图1、图2。

通过图1，可以看到标准化后数据的曲线近似为钟形；从图2得到数据的样本点基本分布在45度参考线内。通过观察，可以认为数据近似服从正态分布。再利用R软件中shapiro.test函数，得到对数据进行Shapiro-Wilk检验的p值为：0.08118 > 0.05。综上，我们认为数据基本符合正态分布。因此，我们可以使用标准化后的数据进行下面的统计分析。

图1. 密度图

图2. QQ图

3. 最小二乘法

3.1. 普通最小二乘法及逐步回归法介绍

3.1.1. 最小二乘法介绍

若有回归方程记为：

$y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_m{x}_m+\epsilon$

其中 ${\beta }_1,{\beta }_2,\cdots ,{\beta }_m$ 为回归系数，也是待估参数。 $\epsilon$ 是随机向量。若 ${\hat{\beta }}_i$ 是参数估计量， $y_k$ 是观测值。那么最小二乘法满足：

$Q=\underset{k=1}{\overset{n}{{\displaystyle \sum }}}{\left[y_k-\left({\hat{\beta }}_0+{\hat{\beta }}_1{x}_1+\cdots +{\hat{\beta }}_m{x}_m\right)\right]}^2$

达到最小，就可以解出 $\hat{\beta }$。 $\hat{\beta }$ 就是 $\beta$ 的最小二乘估计。

3.1.2. 逐步回归法介绍

逐步回归法是利用最小二乘法原理，通过建立正规方程，来解回归系数的一种回归方法。他是向前选择变量法和向后剔除变量法的组合，每进行一步只从中选取一个变量，能够使得到的回归方程更加简化，变量之间的相关性增大 [4]。

3.2. 普通最小二乘法建模

尽管经过检验，本例不适合使用普通最小二乘法，但仍使用最小二乘法进行分析。将普通最小二乘的结果作为之后各种方法的参考，与之进行比较。

利用R软件中lm函数，可以得到模型为：

$Y=-0.58M+0.11G-0.42P-0.26B-0.41I+0.74S$

使用普通最小二乘法得到的描述性统计量见表3。

表3. 普通最小二乘回归的描述性统计量

*代表p值在0.01到0.05范围内，**代表p值在0.001到0.01范围内，***代表p值小于0.001，下同。

由表3可知，只有变量I与变量S的p值小于0.05，是显著的。因此，由普通最小二乘法得到的模型为：

$Y=-0.41I+0.74S$

由于数据进行了标准化，故将系数还原后得到：

$Y=292.9536-0.0005161216I+0.1612487S$

3.3. 逐步回归法建模

在本例中，初始模型包含所有变量。根据AIC值选出最优模型，做逐步回归。利用lm函数与step函数，可以得到逐步回归法建立的模型为：

$Y=-0.5M-0.57P-0.49I+0.69S$

此模型相应的描述性统计量见表4。

表4. 逐步回归法建造模型的描述性统计量

各个自变量均显著，其p值均小于0.05，因此最终模型即为上式。将标准化后的数据还原，得到的模型为：

$Y=394.8733-0.001802806M-0.008439665P-0.0006168282I+0.1503535S$

4. 岭回归

4.1. 岭回归介绍

岭回归估计是改进最小二乘估计的一种算法，可以解决最小二乘法在求解系数向量时矩阵无法求逆的缺点 [5]。岭回归估计放弃了无偏性和部分精确度，可以得到效果稍差但更符合实际的结果。岭回归非常灵活，使用时存在一定的主观性，但这种主观人为的性质正好可以使得定性分析与定量分析进行有机结合，在解决多重共线性问题时有特殊的作用 [6]。应用岭回归方法时，通常要使用标准化后的数据进行计算。首先用普通最小二乘法可以得到正规方程为 [7]：

$r_{xx}b=r_{yX}$

其中， $r_{xx}$ 是X的相关系数矩阵， $r_{yX}$ 是y与X的相关系数的向量。因此，在岭回归中可以使用 $k\left(k\ge 0\right)$ 作为估计量的偏差，得到岭回归的正规方程是：

$\left(r_{xx}+kI\right)b_R=r_{yX}$

k代表了岭回归中估计量的偏差值，这也是岭回归的特点。若 $k=0$，上式为普通最小二乘回归方法；若 $k>0$，则是岭回归方法，并且使得估计量 $b_R$ 的均方误差小于最小二乘回归估计量b的均方误差。如果选取了理想的k值，就可以使估计量的偏差与方差组合后达到最好效果。因此，在岭回归中选择合适的k值来进行计算是非常重要的。

4.2. 岭回归建模

在R软件中加载MASS包，通过lm.ridge函数对标准化后的数据做岭回归。在选取λ值时，首先在0~5范围内选取，精确到0.1，运算得到最优λ值为4.6。再设置范围为4~5，精确到0.01，得到最优λ值为4.65。因此，选择λ值为4.65。利用图片直观确认λ值时，可以使用岭迹图，并且做出λ值与残差平方和的关系图。如图3所示，从图3中左图得到，回归系数都相对稳定；由图3的右图得到在λ取值为4.65时，残差平方和相对较小。

图3. 岭迹图及λ值与残差平方和关系图

确定 $\lambda =4.65$ 后，使用岭回归方法得到的模型为：

$Y=-0.01M-0.04P-0.06B-0.1I+0.01S$

岭回归模型的描述性统计量见表5。

表5. 岭回归法建造模型的描述性统计量

由表5得只有变量B与变量I显著，因此使用岭回归方法对数据进行拟合的模型为：

$Y=-0.06B-0.1I$

岭回归使用标准化后数据进行建模，故还原自变量的系数为：

$Y=312.9246-0.02380766B-0.0001258833I$

5. Lasso回归

5.1. Lasso回归介绍

Lasso回归是一种“降维”思想的方法。在进行Lasso回归时，通过构造一个惩罚函数，可以减小变量的系数甚至将系数降为0，这是岭回归无法做到的。Lasso回归可以适用于线性与非线性两种情况，其原理是在普通的线性模型后，增加惩罚项L1，通过惩罚项来调节系数，降低维度。

假设因变量为： $y={\left(y_1,y_2,\cdots ,y_n\right)}^{\text{T}}$，自变量为： $X={\left(x_{1j},x_{2j},\cdots ,x_{nj}\right)}^{\text{T}},j=1,2,\cdots ,p$，$\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_n\right)}^{\text{T}}$ 是系数向量，考虑线性模型 [8]：

$y=\beta X+\epsilon$

Lasso方法的变量选择和参数估计可以由下式得到，其中λ是正则化参数：

$\beta \left(\text{Lasso}\right)=\arg \min {\displaystyle \underset{i=0}{\overset{n}{\sum }}{\left(y_i-{\displaystyle \underset{j=0}{\overset{p}{\sum }}{\beta }_j{x}_{ij}}\right)}^2}+\lambda {\displaystyle \underset{j=0}{\overset{p}{\sum }}\left|{\beta }_j\right|}$

在求解上式时，可以转化为带有惩罚项的优化问题，其中t与λ相对应，是调整参数：

$\beta \left(\text{Lasso}\right)=\arg \min {\displaystyle \underset{i=0}{\overset{n}{\sum }}{\left(y_i-{\displaystyle \underset{j=0}{\overset{p}{\sum }}{\beta }_j{x}_{ij}}\right)}^2}\text{s}\text{.t}.{\displaystyle \underset{j=0}{\overset{p}{\sum }}\left|{\beta }_j\right|}\le t$

5.2. Lasso回归建模

由于Lasso回归要求数据为矩阵形式，因此设x与y为数据标准化后的自变量与因变量数据的矩阵。使用glmnet函数进行回归拟合，并通过交叉验证选择λ值，得到图4、图5。图4展示了随着λ值的变化，变量系数的变化。图5中两条虚线代表了均方误差最小时的λ值及距离均方误差最小值一个标准误时得λ值。

图4. λ值与变量系数的变化

图5. λ值与均方误差

由图及计算可得，在Lasso回归中可取： ${\lambda }_1=0.0322$，${\lambda }_2=0.3964$。在 $\lambda =0.3964$ 时，模型含有五个变量，因此选择 ${\lambda }_2$ 作为 $\lambda$。得到的模型为：

$Y=-0.48M-0.42P-0.27B-0.4I+0.76S$

模型的描述性统计量见表6。

表6. Lasso回归法建造模型的描述性统计量

由表6可知，变量M、I、S是显著的，其p值小于0.05。因此，使用Lasso回归方法建造的模型为：

$Y=-0.48M-0.4I+0.76S$

使用标准化后的数据进行计算，将系数还原为：

$Y=347.192-0.001730694M-0.0005035332I+0.1656068S$

6. 主成分回归法

6.1. 主成分回归法介绍

主成分回归也是一种“降维”的方法。主成分回归可以使得原始数据信息损失最少，在这种前提下，通过线性变换将原始自变量的集合从高维空间映射到低维空间 [9]。通过对原始自变量的线性组合组成新的主成分变量，利用新的主成分对数据进行普通最小二乘回归，然后再回到原来的自变量。主成分回归适用于数据存在多重共线性的情况。但主成分回归有时难以解释每个主成分的具体含义，有时需要对得到的主成分进行因子旋转，来解释得到的主成分。

主成分回归的步骤为：

1) 首先对数据进行标准化，得到标准化后的矩阵X，及其对应的协方差阵Z。

2) 求Z的前m个特征值 ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_m$，及对应的特征向量 $a_1,a_2,\cdots ,a_m$，它们是标准正交的。

3) 求第i主成分 $F_i$，并且满足 $F_i=X{a}_i,i=1,2,\cdots ,m$。

4) 确认在第n个主成分时，前n个主成分可以完全解释要求的信息量。再用最小二乘法对前n个主成分进行回归。

6.2. 主成分回归法建模

进行主成分回归时，首先需要选择主成分。使用validationplot函数可视化RMSE (均方根误差)，并利用交叉验证检验RMSE。可以得到表7、表8，表7中展示了交叉检验的RMSE，表8中为响应变量能被主成分解释方差的百分比。图6为可视化RMSE。图7为碎石图。

表7. 交叉验证的RMSE

表8. 响应变量能被主成分解释方差的百分比

图6. 可视化RMSE

图7. 碎石图

综合上述表格及图片，选择三个主成分较为合适。此时，误差较小，同时再增加主成分个数对方差的解释量不会增加太多。因此，选择三个主成分进行最小二乘回归。通过princomp函数得到，前三个主成分的累计贡献率为94.6%。用Y对前三个主成分进行普通最小二乘回归，可以得到主成分回归的模型为：

$Y=-0.13Z_1+0.25Z_2-0.51Z_3$

其中

$Z_1=0.40M+0.40G+0.42P+0.45B+0.12I+0.54S$

$Z_2=0.48M+0.48G-0.41P-0.38B-0.48I$

$Z_3=0.26M+0.27G-0.35P+0.80I-0.30S$

同时，主成分回归的描述性统计量见表9。

表9. 主成分回归法构建模型的描述性统计量

$Z_1$ 不显著，因此，最终主成分回归得到的模型为：

$Y=-0.0126M-0.0177G+0.076P-0.095B-0.528I+0.153S$

由于数据进行了标准化，应将模型的系数进行还原。使用原始数据的标准差及均值等，将数据代入模型中，最终得到：

$Y=328.0003+0.001125289P-0.03769546B-0.0006646639I+0.03333925S$

7. 各种回归方法精度的比较

通过上述计算及建模，利用五种方法，将空气质量数据作为原始数据，我们得到了五个模型。将数据代入模型表达式，可以得到每种方法对应的预测结果，见附表3。通过与真实值进行比较，可以得到表10。表10中数据是通过计算得到模型的平均绝对误差与RMSE (均方根误差)，可以对各个模型的精度进行比较。

表10. 各个模型精度的比较

由表中可以得到，从平均绝对误差和RMSE两方面，逐步回归法与主成分回归法表现最好。普通最小二乘法表现最差，这也验证了最小二乘法在存在多重共线性及实际问题中，不如一些对最小二乘法进行改进的方法精度更高、适应性更强。岭回归和Lasso回归的表现比较接近。我们知道岭回归、Lasso回归、主成分回归在应对多重共线性时，可以有效减弱多重共线性的影响。在本文结果中也得到体现。

8. 结论

通过对模型精度的比较，逐步回归方法得到的模型较为准确。我们先对逐步回归方法获得的模型进行分析。逐步回归法得到的模型为：

$Y=394.8733-0.001802806M-0.008439665P-0.0006168282I+0.1503535S$

结合自变量的含义，我们得到空气质量与人均可支配收入、人口数、工业产品产量呈负相关，与政府科学技术预算支出呈正相关。这与常识相吻合。尽管有些省份人均可支配收入较高，但是空气质量不是很好；相反，在一些经济欠发展的省份，空气质量较好。同样地，空气质量与工业产品产量息息相关，这是因为进行工业生产时，往往会产生大量的有毒有害气体，对空气造成污染，例如：炼钢等。

结合其他方法生成的模型，可以得到其他模型自变量系数的分析与上述逐步回归模型的分析基本相同。但是注意到在主成分回归模型与岭回归模型中，环境保护预算支出的系数为负值，这与常识不符。我认为是本文中建模不够准确，所使用的数据量较小造成的失误，还可以进一步改善。

除此之外，综合来看所有模型，可以发现自变量I与自变量S基本出现在所有模型中，这说明工业产品产量与科学技术预算支出对空气质量的影响比较大。尤其是工业产品产量，这提示我们，要想改善环境质量、提高空气质量，进行产业结构的调整比较重要。应减小经济对第二产业的依赖，才能更好地做到经济与环境协调发展。同时，提高对科学技术预算的支出，加大对高新产业的支持，也能达到改善环境与空气质量的目的。

文章引用

刘一鹤,牟唯嫣,金 童. 线性回归方法在空气质量影响因素分析中的应用
Application of Linear Regression Methods in the Analysis of Air Quality Influencing Factors[J]. 应用数学进展, 2022, 11(08): 5936-5950. https://doi.org/10.12677/AAM.2022.118625

参考文献