python数据分析实战：近地面臭氧浓度预测

1.背景

在上海隔离期间闲来无事，做了一个数据分析小项目，项目来源和鲸社区。以下是赛题描述：

近地面的臭氧是蓝天白云下的隐形杀手，高浓度的臭氧对人体健康有很大危害。
近年来，在全球变暖和城市化背景下，夏季极端高温频发，伴随着人为源排放的增加，为臭氧污染提供了有利的前提物和发生条件。臭氧污染存在非线性化学响应关系，其形成与其前体挥发性有机化合物(VOCs)和氮氧化物(NOx)的总量和比例密切相关，也可与颗粒物等其他污染物相互作用；臭氧污染具有明显的区域性特征，对气象因素极其敏感，受到局地的温度、相对湿度、风向、风速等气象条件影响较大。
本赛题需要建立基于气象要素和污染物浓度观测资料的逐小时臭氧污染预测模型。

提供数据：
气象站观测资料 train_weather.csv：

字段	数据类型	字段说明	单位
id	整型	数据标识
time	字符串	时间
pressure	浮点型	气压	hPa
wd	整型	10分钟平均风向	°
ws	浮点型	10分钟平均风速	m/s
tem	浮点型	气温	℃
rh	整型	相对湿度	%
rain	浮点型	降水量	mm

空气污染数据 train_air.csv：

字段	数据类型	字段说明	单位
id	整型	数据标识
PM2.5	整型	细颗粒物	µg/m³
PM10	整型	粗颗粒物	µg/m³
NO2	整型	二氧化氮	µg/m³
SO2	整型	二氧化硫	µg/m³
CO	浮点型	一氧化碳	µg/m³
O3	整型	臭氧	mg/m³

2.数据处理

在对数据进行处理之前，我们先要思考以下问题：

1. 原始数据中包含哪些信息？
2. 我们需要什么形式的数据？
3. 面对异常数据该如何处理？

查看一下原始数据：

原始数据中包含哪些信息？

train_weather中包括了站点类型、时间（精确到小时）、气象数据（气压、风向、风速、气温、相对湿度、降水量）。train_air的格式有所区别，但也包括了站点类型、时间（精确到小时），和空气污染数据（细颗粒物、粗颗粒物、二氧化氮、二氧化硫、一氧化碳、臭氧）。
我们将两组数据包含的信息合并且，得到了一个信息集：

离散	连续
year, month, day, hour, station	pressure, wd, ws, tem, rh, rain, PM2.5, PM10, NO2, SO2, CO, O3

我们需要什么形式的数据？

我们要的是特征集和结果集，特征集每一列代表一个特征向量，而这里的结果集就是O3的浓度。对于train_weather无须多言，已经是我们想要的形式，但train_air的形式特殊，需要进行转换。
在这里引入两个概念，长数据和宽数据：

长数据一般是指数据集中的变量没有做明确的细分，即变量中至少有一个变量中的元素存在值严重重复循环的情况（可以归为几类），表格整体的形状为长方形，即变量少而观察值多。
宽数据是指数据集对所有的变量进行了明确的细分，各变量的值不存在重复循环的情况也无法归类。数据总体的表现为变量多而观察值少。

我们最终要的是宽数据形式，即每列对应着一个特征向量。所以对于train_air整体上要做一个长—>宽的变换，具体为将‘type’列中‘PM2.5’、‘PM10’、‘SO2’等转化为新的特征列。但我们又发现‘station A’和‘station B’列要做一个宽—>长的操作，才能与train_weather相对应。
具体可以参考这个帖子：pandas长宽转换，df.melt进行宽转长，df.pivot进行长转宽。
具体代码如下：

df=df.rename(columns={'Station A':'A','Station B':'B'})
#melt函数宽数据转长数据
df=df.melt(['type','time'])
#pivot函数长数据转宽数据
df=df.pivot(index=['time','variable'],columns='type',values='value')
df=df.reset_index()
df=df.rename(columns={'variable':'station'})

其次一个问题，我们发现两组数据的时间格式并不一致，需要统一。在这里有两种做法，一是合并，二是拆分。时间格式详见这篇文章时间格式。

• 合并
  train_air的时间数据实际上有两列，‘yyyyMMdd’与int型表示的小时。train_weather的时间数据是一列‘yyyy/MM/dd HH:mm’。我将其都转化为‘yyyy-MM-dd HH:mm:ss’的格式。
  下面是我的操作，思路为，对于需要拼接的时间数据，先以字符串的格式对其组装，再用pd.to_datetime函数将其转换为标准格式。

#对train_air的操作
s=pd.to_datetime(df['date'],format='%Y%m%d').astype(str)+pd.to_datetime(df['hour'],format='%H').astype(str).str[10:]
df['time']=pd.to_datetime(s)
#对train_weather的操作
df['time']=pd.to_datetime(df['time'])

• 拆分
  把‘time’直接拆成‘year’，‘month’，‘day’，‘hour’四列，都以int的形式表示。这样做是为了尽可能精细化数据，提升后续模型训练的精确度。我是先进行合并，再进行的拆分，保证数据的一致性。

train['year'] = train.time.dt.year
train['month'] = train.time.dt.month
train['day'] = train.time.dt.day
train['hour'] = train.time.dt.hour

面对异常数据该如何处理？

异常数据主要分为两种情况：一是空值、NaN或Null，在这里统称为缺失值；二是不符合常理的离群点，如过大或过小，或者inf。
说是两种情况，真正要处理的只有一，因为遇到第二种情况，直接把异常值转化为NaN即可。如何检测离群点又是一个知识体系了，自行查看：离群点检测。但离群点检测这个操作并非在所有情况都必要，尤其是使用一些鲁棒性较强的模型时。
如何处理空值？其实也是需要具体问题具体分析。但可以总结为三种操作：删除、填充、保留。

• 删除
  直接删除带空值的行或者缺失值过多的特征是一种简单粗暴的方式，优点是快捷且省心，因为做减法是无功但也无过，保留数据本来的面目，缺点是可能会损失太多信息。在缺失值较少的情况下，直接删除是个好方案，但空值过多时便要三思了。

• 填充
  空值填充
  常见的用均值、中位数或众数替代，或用插值法。也可用算法进行拟合填充。填充虽然尽可能保留了原有信息，但也会对数据的准确性造成一定影响。

• 保留
  先谈一下几类模型对缺失值的敏感程度，树模型（decision tree，random forest，xgboost，lightgbm）对于缺失值的敏感度较低，大部分时候可以在数据有缺失时使用。涉及到距离的算法（KNN，SVM）对缺失数据比较敏感。线性模型代价函数也往往涉及到距离的计算，所以也比较敏感。神经网络和贝叶斯模型对于缺失数据也不是非常敏感，但神经网络需要的数据量大，贝叶斯可以用于小样本。
  所以针对要使用的模型，自行决定是否保留。

3.模型训练

模型选择

各大回归模型的优缺点各大网站都能搜到，在这里仅谈我的个人经验。
我一般只从三种模型中选择：Random Forest以及boosting方法下的Xgboost、Lightgbm，主要原因就是经广泛认可效果好，在实际中大量应用。
这次直接选择了Lightgbm，相比于Xgboost，具备更优良的性能，比如计算快捷、消耗内存少、泛化能力也更强。
模型的选择见仁见智、具体问题具体分析。但如果像我一样，以解决实际问题为目的，不追求科研上的严谨性，那就直接选择应用最广泛的即可，不用在此过多纠结，如果效果不好再去换另一个。
python的Lightgbm方法有两个重要的优点：一是训练数据可以携带空缺值，二是不用在之前特意进行特征的选择。

交叉验证

在算力允许的情况下，一定要使用交叉验证，来避免过拟合。在sklearn.model_selection中，有train_test_split和cross_validate方法，前者用来快速分割训练集和测试集，后者就是用来交叉验证的。
现有的python集成框架已经相当完备，模型建立+交叉验证只需要两三行代码：

LGBM=lgb.LGBMRegressor(n_estimators=200)
scores=cross_validate(LGBM,X,Y,cv=5,scoring='r2',n_jobs=-1)
scores['test_score'].mean()

这里采用五折交叉验证，衡量回归好坏的指标使用的是R-square（也可以用均方根误差等进行衡量，详见scoring取值），产生五个数取平均。注意要使n_jobs=-1，它代表工作的core数量，等于-1的时候，表示cpu里的所有core进行工作，运算更快。

调参

关于Lightgbm的参数调整方式，这个文章说的很完备了：Lightgbm调参指南。总结起来就是：

按经验预先固定的参数：
learning_rate
n_estimators
min_split_gain
min_child_sample
min_child_weight
需要算法调节的参数 ：
max_depth
num_leaves
subsample
colsample_bytree
reg_alpha
reg_lambda

在这介绍一种调参的方式：网格调参（Grid Search）。实际上就是枚举法，把可能的参数列成表，看哪个模型评分最高。
我用sklearn.model_selection中的GridSearchCV方法对reg_alpha进行网格调参，最后的结果reg_alpha=0时最佳，也就是默认值。

param_test = {'reg_alpha':[0,2,5,10,15,20]}
gsearch = GridSearchCV(
lgb.LGBMRegressor(),
param_grid=param_test,
scoring='r2',
iid=False,
cv=5,
n_jobs=-1)
gsearch.fit(X,Y)
print(gsearch.best_params_, gsearch.best_score_)

从我的实际操作来看，调参对模型的提升比较微小，让R-square提升0.005就算是很不错了。如果想大幅度优化，还是得从结构上出发，换模型或者修改数据形式。

优化

调参也是优化的一环，但由于其提升甚微，在这里介绍一些我总结的其他方法。这些方法并不一定会带来正向效果，但至少会提供一些新思路。

1. 充分挖掘原始数据中的信息
   在一开始时，我是将A、B两个站点分开预测的，R-suqare都仅徘徊在0.67左右，但我将A、B的数据合并进行训练，R-square直接提升到了0.78。为什么产生了1+1>2的效果？我想是因为A、B两个站点本质差异并不大，换句话说，二者并非相互独立关系，都能给对方提供信息上的补充。所以这说明模型训练不好，可能就是因为提供的信息太少，巧妇难为无米之炊！
   同理我还做了如下操作：保留原数据中的有空缺值的样本、将时间数据更加细化。都是在提供更多信息。
2. 尝试数据的变换
   尝试特征集归一及标准化，以及查看因变量是否正态，如果不是正态尝试变换调整。

$Y$	$\sqrt{Y}$

我发现了根号下的Y是更接近正态的，用其代替Y，效果并没变好但值得尝试。

结果

最后汇报一下我的训练结果：

R-square最终到了0.8以上，这个结果说明拟合效果是不错的。虽然R-square会随着特征数量增加而提升，但从0.67提升到了0.8以上，足以反映我们的优化是很有成效的。
下面是利用模型进行的预测：

代码

import math
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split,cross_validate,GridSearchCV
import lightgbm as lgb
import seaborn as sns
train_air=pd.read_csv('/home/mw/input/ozone3557/train_air.csv')
train_weather=pd.read_csv('/home/mw/input/ozone3557/train_weather.csv')
test_air=pd.read_csv('/home/mw/input/ozone3557/test_air.csv')
test_weather=pd.read_csv('/home/mw/input/ozone3557/test_weather.csv')
def clean(df):
    #转化为统一的时间格式
    if 'date' in df.columns:
        #对train_air的操作
        s=pd.to_datetime(df['date'],format='%Y%m%d').astype(str)+pd.to_datetime(df['hour'],format='%H').astype(str).str[10:]
        df['time']=pd.to_datetime(s)
        del df['date']
        del df['hour']
        df=df.rename(columns={'Station A':'A','Station B':'B'})
        #melt函数宽数据转长数据
        df=df.melt(['type','time'])
        #pivot函数长数据转宽数据
        df=df.pivot(index=['time','variable'],columns='type',values='value')
        df=df.reset_index()
        df=df.rename(columns={'variable':'station'})
    else:
        df['time']=pd.to_datetime(df['time'])
        f1=df['pressure']<1e5
        f2=df['wd']<1e5
        f3=df['ws']<1e5
        f4=df['tem']<1e5
        f5=df['rh']<1e5
        f6=df['rain']<1e5
        df=df.where(f1&f2&f3&f4&f5&f6)
    return df
def scale(df):
    if 'O3' in df.columns:
        df_y = df['O3']
        df_x_cate=df[['year','month','day','hour','station']].astype('category')#规定分类特征
        df_x_num = df.drop(['time', 'O3','year','month','day','hour','station'], axis=1)
        zscore = StandardScaler()
        zscore = zscore.fit_transform(df_x_num)
        df_x_num=pd.DataFrame(zscore,index=df_x_num.index,columns=df_x_num.columns)
        df_x=pd.concat([df_x_cate,df_x_num],axis=1)
        return df_x, df_y
    else:
        df_x_cate=df[['year','month','day','hour','station']].astype('category')#规定分类特征
        df_x_num = df.drop(['time','year','month','day','hour','station'], axis=1)
        zscore = StandardScaler()
        zscore = zscore.fit_transform(df_x_num)
        df_x_num=pd.DataFrame(zscore,index=df_x_num.index,columns=df_x_num.columns)
        df_x=pd.concat([df_x_cate,df_x_num],axis=1)
        return df_x
 
train_air=clean(train_air)
test_air=clean(test_air)
train_weather=clean(train_weather)
test_weather=clean(test_weather)
#合并两组数据并清除结果列有空缺值的行
train=pd.merge(train_weather,train_air,on=['station','time'],how='left')
train=train[~train.O3.isna()]
train['year'] = train.time.dt.year
train['month'] = train.time.dt.month
train['day'] = train.time.dt.day
train['hour'] = train.time.dt.hour
test=pd.merge(test_weather,test_air,on=['station','time'],how='left')
test['year'] = test.time.dt.year
test['month'] = test.time.dt.month
test['day'] = test.time.dt.day
test['hour'] = test.time.dt.hour
X,Y=scale(train)
test_X=scale(test)

#交叉验证
LGBM=lgb.LGBMRegressor(n_estimators=200)
scores=cross_validate(LGBM,X,Y,cv=5,scoring='r2',n_jobs=-1)
scores['test_score'].mean()
#网格调参
param_test = {'reg_alpha':[0,2,5,10,15,20]}
gsearch = GridSearchCV(
    lgb.LGBMRegressor(),
    param_grid=param_test,
    scoring='r2',
    iid=False,
    cv=5,
    n_jobs=-1)
gsearch.fit(X,Y)
print(gsearch.best_params_, gsearch.best_score_)
#预测
LGBM=lgb.LGBMRegressor(n_estimators=200,n_jobs=-1)
LGBM.fit(X,Y)
y_pred = LGBM.predict(test_X)
print(y_pred)