Python实现城市空气质量分析与预测

AQI分析与预测

背景信息

AQI指的是空气质量指数，用来衡量一个城市的空气清洁或污染的程度，数值越小则空气质量越好。

任务说明

• 哪些城市的空气质量较好 / 较差？【描述性统计分析】
• 空气质量在地理上的分布，是否存在着一定的规律？【描述性统计分析】
• 沿海城市与内陆城市的空气质量是否存在不同？【推断统计分析】
• 空气质量主要受到哪些因素的影响？【相关系数分析】
• 全国空气质量的整体情况是怎样的？【区间估计】

• 怎么来预测一个城市的空气质量？【线性回归】

  数据集描述

  现在获取了2015年的空气质量数据集，该数据集包含了全国主要城市的相关数据及空气质量指数，数据情况如下：

  列名	含义
  City	城市名
  AQI	空气质量指数
  Precipitation	降雨量
  GDP	城市生产总值
  Temperature	温度
  Longitude	纬度
  Latitude	经度
  Altitude	海拔高度
  PopulationDensity	人口密度
  Coastal	是否沿海
  GreenCoverageRate	绿化覆盖率
  Incineration(10000ton)	焚烧量(10000吨)

  最好/最差的5个城市

  data = pd.read_csv("data.csv")
  data.drop_duplicates(inplace=True)
  
  ## 最好的5个城市
  t = data[["City", "AQI"]].sort_values("AQI")
  display(t.iloc[:5])
  plt.xticks(rotation=45)  # x 轴刻度旋转度数（与水平方向的逆时针夹角度数）
  sns.barplot(x="City", y="AQI", data=t.iloc[:5])
  
  ## 最差的5个城市
  t = data[["City", "AQI"]].sort_values("AQI", ascending=False)
  display(t.iloc[:5])
  plt.xticks(rotation=45)
  sns.barplot(x="City", y="AQI", data=t.iloc[:5])

  全国城市空气质量

  城市空气质量等级统计

  对于AQI，可以对空气质量进行等级划分，划分标准如下：

  AQI指数	等级	描述
  0-50	一级	优
  51-100	二级	良
  101-150	三级	轻度污染
  151-200	四级	中度污染
  201-300	五级	重度污染
  >300	六级	严重污染

根据该标准，来统计下全国空气质量每个等级的数量。

## 将AQI转换为对应的等级
def value_to_level(AQI):
    if AQI >= 0 and AQI <= 50:
        return "一级"
    elif AQI >= 51 and AQI <= 100:
        return "二级"
    elif AQI >= 101 and AQI <= 150:
        return "三级"
    elif AQI >= 151 and AQI <= 200:
        return "四级"
    elif AQI >= 201 and AQI <= 300:
        return "五级"
    else:
        return "六级"

level = data["AQI"].apply(value_to_level)
display(level.value_counts())
sns.countplot(x=level, order=["一级", "二级", "三级", "四级", "五级"])

可见，城市空气质量主要以一级(优)与二级(良)为主，三级(轻度污染)占一部分，更高污染的城市占少数。

空气质量分布

接下来绘制下全国各城市的空气质量分布图。

## hue：分组依据的类别变量
## palette：在对数据进行分组时，设置不同组数据的显示颜色
## RdYlGn_r：红黄绿的渐变颜色，r是reverse的意思
sns.scatterplot(x="Longitude", y="Latitude", hue="AQI", palette=plt.cm.RdYlGn_r, data=data)

从结果上可以发现，从大致的地理位置，西部城市要好于东部城市，南部城市要好于北部城市。

临海空气质量优于内陆

数量统计

先来统计下临海城市与内陆城市的数量：

display(data["Coastal"].value_counts())
sns.countplot(x="Coastal", data=data)

分布统计

然后来观察临海城市与内陆城市的散点分布：

sns.stripplot(x="Coastal", y="AQI", data=data)  # 分类散点图
sns.swarmplot(x="Coastal", y="AQI", data=data)  # 分簇散点图（蜂群图）

然后分组计算空气质量的均值：

display(data.groupby("Coastal")["AQI"].mean())
sns.barplot(x="Coastal", y="AQI", data=data)  # 柱状图

在柱状图（图中均值的竖线是均值的置信区间）中，仅显示了内陆城市与沿海城市的空气质量指数（AQI）的均值对比，可以用箱线图来显示更多的信息。

sns.boxplot(x="Coastal", y="AQI", data=data)

也可以绘制小提琴图，除了能够展示箱线图的信息外，还能呈现出分布的密度。

sns.violinplot(x="Coastal", y="AQI", data=data)

也可以同时绘制蜂群图和小提琴图：

sns.violinplot(x="Coastal", y="AQI", data=data, inner=None)
sns.swarmplot(x="Coastal", y="AQI", color="g", data=data)

两独立样本t检验

为了进一步验证结果，可以使用两独立样本t检验，来查看临海城市与内陆城市的均值差异是否显著。

在统计学中，两独立样本t检验的前提条件为:

1. 样本观测值之间是独立的
2. 两个变量服从正态分布

同时，还需要考虑方差是否相同，方差是否相同，检验结果也会有所差异。

正态分布检验

想判断一组数据是否来自于正态分布的总体，可以采用以下的几个方式进行验证：

1. 绘制数据的分布图
2. 绘制PP图与QQ图
3. 使用假设检验

首先，绘制两个样本的分布图：

from scipy import stats

coastal = data[data["Coastal"] == 1]["AQI"]
inland = data[data["Coastal"] == 0]["AQI"]

fig, ax = plt.subplots(1, 2)
fig.set_size_inches(15, 5)
sns.distplot(coastal, ax=ax[0])
sns.distplot(inland, ax=ax[1])

PP图与QQ图

PP图(Probability-Probability plot)与QQ图(Quantle-Quantile plot)本质上基本是相同的。用于检验样本数据的分布是否符合某个分布(默认为正态分布)。

import statsmodels.api as sm

def plot_pp_qq(d):
    '''绘制PP图与QQ图
    d: 要绘制的数值，array-like
    '''
    fig,ax=plt.subplots(1,2)
    fig.set_size_inches(15,5)
    scale_data=(d-d.mean())/d.std()
    p=sm.ProbPlot(data=scale_data,dist=stats.norm)
    p.ppplot(line='45',ax=ax[0])
    ax[0].set_title('PP图')
    p.qqplot(line='45',ax=ax[1])
    ax[1].set_title('QQ图')
    plt.show()
plot_pp_qq(coastal)
plot_pp_qq(inland)

从结果上可以发现，不太可能服从正态分布。

假设检验

## 原假设，观测值来自于正态分布的总体
print(stats.normaltest(coastal))
print(stats.normaltest(inland))

P<0.05，拒绝原假设，说明两个样本都不是来自于正态分布的总体。

非正态总体的检验

经过一系列的周折，结果令人失望，两样本同时服从正态分布的可能性几乎为0。这似乎意味着之前所作出的一切努力付诸东流。

在假设检验中，很多假设都是以总体为正态分布为前提的，例如Z检验与t检验。当总体并非正态分布时，可以使用非参数检验。

非参数检验

可以使用非参数检验，来实现两样本是否来自同一分布。

## 曼-惠特尼检验。应该仅在每个停本答量 > 20时使用
## 原假设，两个样本服从相同的分布
print(stats.mannwhitneyu(coastal, inland, alternative="two-sided"))  # 默认双边检验
## 威尔科克森秩和检验
## 原候设，两个样本服从相同的分布
print(stats.ranksums(coastal, inland))

P<0.05，拒绝原假设，说明两个样本是服从不同的分布的，意味着它们相应的均值中位数是不同的，也就意味着临海城市与内陆城市它们的空气质量并非是相同的。

主要受哪些因素影响

散点图矩阵

可以绘制散点图矩阵，该矩阵会显示两两变量之间的关系。因为变量较多，图像会变得很小，故这里仅选择部分变量进行绘制。

sns.pairplot(data[["AQI", "PopulationDensity", "GreenCoverageRate"]])

协方差

协方差体现的是两个变量之间的离散程度以及两个变量变化步调是否一致。

相关系数

相关系数，用来体现两个数值变量之间的相关性，是协方差的标准化版本，最为常用的为皮尔逊相关系数，衡量两个数值变量的线性相关性。

相关系数的取值范围为[-1,1]，可以根据相关系数的取值(绝对值)来衡量两个变量的相关性：

• 0.8-1.0：极强相关
• 0.6-0.8：强相关
• 0.4-0.6：中等程度相关
• 0.2-0.4：弱相关

• 0.0-0.2：极弱相关或无相关

  x = data["AQI"]
  y = data["Precipitation"]
  ## 计算协方差
  a = (x - x.mean()) * (y - y.mean())
  cov = np.sum(a) / (len(a) - 1)
  print("协方差：", cov)
  ## 计算相关系数
  corr = cov / np.sqrt(x.var() * y.var())
  print("相关系数：", corr)
  
  ## series对象也提供了cov和corr方法
  print("协方差：", x.cov(y))
  print("相关系数：", x.corr(y))
  ## dataframe对象也提供了corr方法，可计算任意两个数值变量的相关系数
  data.corr()
  
  ## 用热图来展示相关系数
  plt.figure(figsize=(15, 10))
  ## annot：是否标记数值
  ax = sns.heatmap(data.corr(), cmap=plt.cm.RdYlGn, annot=True, fmt=".2f")

  AQI跟纬度有0.55的正相关，也就是说，随着纬度的增大，AQI也是增大的，即越北空气质量越差，说明南方城市的空气质量好一些。

  结果统计

  从结果中可知，空气质量指数主要受降雨量(-0.40)与纬度(0.55)影响。

• 降雨量越多，空气质量越好
• 纬度越低，空气质量越好

此外，还能够发现其他一些明显的细节:

• GDP (城市生产总值) 与Incineration (焚烧量)正相关(0.90)
• Temperature (温度)与Precipitation (降雨量)正相关(0.69)
• Temperature (温度)与Latitude (纬度)负相关(-0.81)
• Longitude (经度)与Altitude (海拔)负相关(-0.74)
• Latitude (纬度)与Precipitation (降雨量)负相关(-0.66)
• Temperature (温度)与Altitude (海拔)负相关(-0.46)

• Altitude (海拔)与Precipitation (降雨量)负相关(-0.32)

  关于空气质量的验证

  网上有传闻说，全国所有城市的空气质量指数均值在71左右，请问，这个消息可靠嘛？

  data["AQI"].mean()
  r = stats.ttest_1samp(data["AQI"], 71)  # 单样本t检验
  print("t值：", r.statistic)
  print("p值：", r.pvalue)

  P>0.05，故无法拒绝原假设，因此接受原假设。

要清楚，+-1.96倍的标准差，是正态分布在置信区间为95%下的临界值，严格来说，对 t 分布不是如此，只不过，当样本容量较大时，t分布近似正态分布。但当样本容量较小时，二者会有较大的差异。可以获取更准确的置信区间：

mean = data["AQI"].mean()
std = data["AQI"].std()
stats.t.interval(0.95, df=len(data) - 1, loc=mean, scale=std / np.sqrt(len(data)))

由此，计算出全国所有城市平均空气质量指数，有95%的可能在70.63-80.04之间。

空气质量预测

是否能够预测该城市的空气质量指数呢？答案是肯定的。可以通过对以往的数据，去建立一种模式，然后将这种模式去应用于未知的数据，进而预测结果。

数据转换

对于模型来说，内部进行的都是数学上的运算，在进行建模之前，我们需要首先进行转换，将类别变量转换为离散变量。

data["Coastal"] = data["Coastal"].map({"是": 1, "否": 0})
data["Coastal"].value_counts()

Coastal变量是一个二分类变量，在对变量进行转换时，可以将“是”与“否”转换为其他离散整数值，例如1和0。

假设某变量是一个多分类变量(三个或是三个以上的分类)，在对变量进行转换时，不可以将n个类别映射为0 1 2 ... n-1。

这时，可以采用独热编码即 One-Hot 编码，又称一位有效编码，其方法是把一个特征拆分为多个特征。

基模型

首先，不进行任何处理，建立一个基模型，作为基线（baseline）。后续的操作，可以在此基础上进行改进。

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

## city（城市名称）对预测毫无用处，删掉
x = data.drop(["City", "AQI"], axis=1)
y = data["AQI"]
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)

lr = LinearRegression()
lr.fit(x_train, y_train)
print(lr.score(x_train, y_train))
print(lr.score(x_test, y_test))

y_hat = lr.predict(x_test)
plt.figure(figsize=(15, 5))
plt.plot(y_test.values, "-r", label="真实值", marker="o")
plt.plot(y_hat, "-g", label="测试值", marker="D")
plt.legend(loc="upper left")
plt.title("线性回归预测结果", fontsize=20)

异常值处理

可以依据箱线图判断离群点的原则去探索异常值，然后使用临界值替换掉异常值。

## Coastal是类别变量，映射为离散变量，不会有异常值
for col in x.columns.drop("Coastal"):
    if pd.api.types.is_numeric_dtype(x_train[col]):
        quartile = np.quantile(x_train[col], [0.25, 0.75])
        IQR = quartile[1] - quartile[0]
        lower = quartile[0] - 1.5 * IQR
        upper = quartile[1] + 1.5 * IQR
        x_train[col][x_train[col] < lower] = lower
        x_train[col][x_train[col] > upper] = upper
        x_test[col][x_test[col] < lower] = lower
        x_test[col][x_test[col] > upper] = upper

去除异常值后，用新的训练集与测试集来评估模型的效果。

lr.fit(x_train, y_train)
print(lr.score(x_train, y_train))
print(lr.score(x_test, y_test))

结果跟之前相比有了一定的少许提升。

数据标准化

数据标准化不会对线性回归的评估效果造成影响，不过，会影响模型的权重系数。如果需要考虑线性回归的权重，则需要对数据首先进行标准化处理。

from sklearn.preprocessing import StandardScaler

s = StandardScaler()
x_train_scale = s.fit_transform(x_train)
x_test_scale = s.transform(x_test)

特征选择

上面我们使用所有可能的原始数据作为特征，建立模型，然而，特征并非越多越好，有些对模型质量没有什么改善的特征，可以进行删除，以提高模型质量。

RFECV模型

特征选择的方式有很多，这里使用RFECV方法来实现特征选择。RFECV分成两个部分，分别是:

• RFE( Recursive feature elimination )：递归特征消除，用来对特征进行重要性评级
• CV( Cross Validation )：交叉验证，在特征评级后，通过交叉验证，选择最佳数量的特征

具体过程如下：

RFE阶段

1. 初始的特征集为所有可用的特征
2. 使用当前特征集进行建模，然后计算每个特征的重要性
3. 删除最不重要的一个(或多个)特征，更新特征集
4. 跳转到步骤2，直到完成所有特征的重要性评级

CV阶段

1. 根据RFE阶段确定的特征重要性，依次选择不同数量的特征
2. 对选定的特征集进行交又验证

3. 确定平均分最高的特征数量，完成特征选择

   from sklearn.feature_selection import RFECV
   
   ## estimator：要操作的模型
   ## step：每次删除的变量数
   ## cv：使用的交叉验证折数（就是把数据分成几份）
   ## scoring：评估方式
   rfecv = RFECV(estimator=lr, step=1, cv=5, scoring="r2")
   rfecv.fit(x_train, y_train)
   
   ## 返回经过选择之后，剩余的特征数量
   print(rfecv.n_features_)
   ## 返回经过特征选择后，使用缩减特征训练后的模型
   print(rfecv.estimator_)
   ## 返回每个特征的等级，数值越小，特征越重要
   print(rfecv.ranking_)
   ## 返回布尔数组，用来表示特征是否被选择
   print(rfecv.support_)
   ## 返回对应数量特征时，模型交叉验证的评分
   print(rfecv.cv_results_)
   
   ## 绘制特征选择过程中，使用交叉验证获取的R^2
   plt.plot(range(1, len(rfecv.cv_results_["mean_test_score"]) + 1), rfecv.cv_results_["mean_test_score"], marker="o")
   plt.xlabel("特征数量")
   plt.ylabel("交叉验证R^2值")

   然后，对测试集应用这种特征选择（变换），进行测试，获取测试结果。

   print("剔除的变量：", x_train.columns[~rfecv.support_])
   x_train_eli = rfecv.transform(x_train)
   x_test_eli = rfecv.transform(x_test)
   print(rfecv.estimator_.score(x_train_eli, y_train))
   print(rfecv.estimator_.score(x_test_eli, y_test))

   可以发现，经过特征选择后，使用剩余特征训练的模型，与之前模型表现上几乎没有变化，也就是说，我们成功删除了一些特征。

   绘制特征权重

   ## 获取列名与对应的权重，构成一个元组，作为列表的元素
   li = list(zip(x.columns.values[rfecv.support_], rfecv.estimator_.coef_))
   ## 根据权重的绝对值，对列表进行降序排列
   li.sort(key=lambda x: abs(x[1]), reverse=True)
   s = pd.Series(dict(li))
   plt.figure(figsize=(15, 5))
   ax = sns.barplot(y=s.index, x=s.values)
   for y, x in enumerate(s):
    # 在条形图上绘制数值
    t = ax.text(x / 2, y, round(x, 3))
    # 教值居中对齐
    t.set_ha("center")
   plt.show()

   分箱操作

   分箱后，我们不能将每个区间映射为离散数值，而是应当使用One-Hot编码。

   from sklearn.preprocessing import KBinsDiscretizer
   
   ## KBinsDiscretizer K个分箱的离散器，用于将数值（通常是连续变量）变量进行区间离散化操作
   ## n_bin：分箱（区间）个数
   ## encode：离散化编码方式，分为：onehot、onehot-dense与ordinal
   ##    onehot：使用独热编码，返回稀疏矩阵
   ##    onehot-dense：使用独热编码，返回稠密矩阵
   ##    ordinal：使用序数编码（0,1,2……）
   ## strategy：分箱的方式，分为：uniform、quantile、kmeans
   ##    uniform：每个区间的长度范围大致相同
   ##    quantile：每个区间包含的元素个数大致相同
   ##    kmeans：使用一维kmeans方式进行分箱
   
   k = KBinsDiscretizer(n_bins=[4, 5, 10, 6], encode="onehot-dense", strategy="uniform")
   ## 定义离散化的特征
   discretize = ["Longitude", "Temperature", "Precipitation", "Latitude"]
   
   x_train_eli = pd.DataFrame(data=x_train_eli, columns=x.columns[rfecv.support_])
   x_test_eli = pd.DataFrame(data=x_test_eli, columns=x.columns[rfecv.support_])
   
   r = k.fit_transform(x_train_eli[discretize])
   r = pd.DataFrame(r, index=x_train_eli.index)
   
   ## 获取除离散化特征之外的其他特征
   x_train_dis = x_train_eli.drop(discretize, axis=1)
   ## 将离散化后的特征与其他特征进行重新组合
   x_train_dis = pd.concat([x_train_dis, r], axis=1)
   ## 对测试集进行同样的离散化操作
   r = pd.DataFrame(k.transform(x_test_eli[discretize]), index=x_test_eli.index)
   x_test_dis = x_test_eli.drop(discretize, axis=1)
   x_test_dis = pd.concat([x_test_dis, r], axis=1)
   ## 查看转换之后的格式
   display(x_train_dis.head())