Question

回归是研究自变量x对因变量y影响的一种数据分析方法。最简单的回归模型是一元线性回归（只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示），可以表示为Y=β₀+β₁x+ε，其中Y为因变量，x为自变量，β₁为影响系数，β₀为截距，ε为随机误差。

回归分析是广泛应用的统计分析方法，可用于分析自变量和因变量的影响关系（通过自变量求因变量），也可以分析自变量对因变量的影响方向（正向影响还是负向影响）。回归分析的主要应用场景是进行预测和控制，例如计划制定、KPI制定、目标制定等方面；也可以基于预测的数据与实际数据进行比对和分析，确定事件发展程度并给未来行动提供方向性指导。

常用的回归算法包括线性回归、二项式回归、对数回归、指数回归、核SVM、岭回归、Lasso等。

回归分析的优点是数据模式和结果便于理解，如线性回归用y=ax+b的形式表达，在解释和理解自变量与因变量关系式相对容易；在基于函数公式的业务应用中，可以直接使用代入法求解，因此应用起来比较容易。回归分析的缺点是只能分析少量变量之间的相互关系，无法处理海量变量间的相互作用关系，尤其是变量共同因素对因变量的影响程度。

4.2.1　注意回归自变量之间的共线性问题

在应用回归模型时，应注意识别和解决自变量间的共线性问题。有关该话题，请参照3.7节了解更多信息。

4.2.2　相关系数、判定系数和回归系数之间到底什么关系

在应用回归模型时经常会得到类似的回归方程：y=42.738x+169.94，其中R2=0.5252。如果对这两个变量做相关性分析，还会得到相关系数R=0.72468551874050。

上述结果中，42.738就是自变量x的回归系数，0.5252则是该方程的判定系数，0.72468551874050是两个变量的相关性系数。这三个值之间到底什么关系？

回归系数：在回归方程中表示自变量x对因变量y影响大小的参数，其绝对值的高低只能说明自变量与因变量之间的联系程度和变化量的比例。

自变量的单位不同会显著影响该系数，这意味着，量级越大的自变量，该系数越高。

判定系数：是自变量对因变量的方差解释程度的值，计算公式为：回归平方和与总离差平方和之比值。

相关系数：也称为解释系数，是衡量变量间的相关程度或密切程度的值，其本质是线性相关性的判断。

三者的相互关系：

判定系数是所有参与模型中自变量的对因变量联合影响程度，而非某个自变量的影响程度。本节示例中只有1个自变量，因此R²代表的是该自变量本身；如果有x1、x2两个自变量，那么R²代表的这两个自变量共同影响的结果。假如在线性回归中只有一个自变量，那么判定系数等于相关系数的平方。

回归系数与相关系数的关系：回归系数＞0，相关系数取值在（0，1），说明二者正相关；如果系数小于0，相关系数取值在（-1，0），说明二者负相关。

对于判定系数的高低没有统一的标准，笼统而言，值越高越能代表自变量对因变量的解释作用越大。普通数据通常达到0.6就已经很高，0.6意味着相关系数至少在0.8以上；而对于时间序列数据，判定系数达到0.9则很平常。

4.2.3　判定系数是否意味着相应的因果联系

在3.8节中我们提到了相关性和因果不是一回事，那么判定系数是因果联系吗？其实也不是，以一元线性回归方程（本节开始的方程）举例：

判定系数是相关系数的平方，既然相关系数不是因果，为什么其平方后就能成为因果？

判定系数的出发点是用来评估整个模型的拟合优度，直白点就是自变量引起的变动占总变动的百分比。那么二者有相同的变动趋势就意味着是因果关系吗？显然也不是。

当然，需要承认在数学关系上，确实反映了当自变量x变化后，y根据方程也会有相应的变化。但是，二者的关系严格上讲无法判断，原因是我们很可能无法知晓回归中的自变量是否真的是导致因变量发生变化的根本因素。只是从现象上，将二者一起/相关变化的规律，使用回归方程来进行解释。

4.2.4　注意应用回归模型时研究自变量是否产生变化

在应用回归模型做预测时，必须研究对因变量产生影响的自变量是否产生变化，主要考察两个方面：

一是是否有产生了新的对因变量影响更大自变量

在建立回归模型时，需要综合考虑自变量的选择问题。如果遗漏了重要的变量，那么模型很可能无法正确反映实际情况，而且参数估计是有偏的，此时的回归模型及其不稳定且方差较大。同样在应用回归模型时，仍然需要重新评估该问题是否发生。例如，在做用户订单金额预测时是基于正常销售状态下的变量实现的；但当发生大型促销活动且促销活动因素没有被纳入到回归模型中时，原来的回归模型则无法有效预测。

二是原有自变量是否仍然控制在训练模型时的范围之内

如果有自变量的变化超过训练模型的范围，那么原来的经验公式可能无法在新的值域范围下适用，这种情况通常需要重新做研究和建模。例如，假设我们建立了一个回归模型，可以广告投放费用预测广告点击率，在训练回归模型时的广告费用在[0，1000]区间内，如果在广告费用超过1000（例如2000）时，则无法保证最终预测效果是否有效。

判定系数经常作为拟合好坏的主要参照指标，当一个新的指标加入模型后发现模型不变，此时无法根据判定系数来反推该指标的重要性程度，例如该指标无效（或有效）。

4.2.5　如何选择回归分析算法

回归算法按照自变量的个数分为一元回归和多元回归，按照影响是否线性分为线性回归和非线性回归。在面对不同回归方法的选择，注意参考以下因素：

入门的开始——简单额线性回归。如果是学习为主，那么不需要选择多么强大的模型，基于最小二乘法的普通线性回归最合适；同时，它适合数据集本身结构简单、分布规律有明显线性关系的场景。

如果自变量数量少或经过降维后得到了可以使用的二维变量（包括预测变量），那么可以直接通过散点图发现自变量和因变量的相互关系，然后选择最佳回归方法。

如果经过基本判断发现自变量间有较强的共线性关系，那么可以使用对多重共线性（自变量高度相关）能灵活处理的算法，例如岭回归。

如果数据集噪音较多，推荐使用主成分回归，因为主成分回归通过对参与回归的主成分的合理选择，可以去掉噪音。另外，各个主成分间相互正交，能解决多元线性回归中的共线性问题。这些都能有效地提高模型的抗干扰能力。

如果在高维度变量下，使用正则化回归方法效果更好，例如Lasso，Ridge和ElasticNet，或者使用逐步回归从中挑选出影响显著的自变量来建立回归模型。

如果要同时验证多个算法，并想从中选择一个来做做好的拟合，使用交叉检验做多个模型的效果对比，并通过R-square、Adjusted R-square、AIC、BIC以及各种残差、误差项指标做综合评估。

如果注重模型的可解释性，那么容易理解的线性回归、指数回归、对数回归、二项或多项式回归要比核回归、支持向量回归机等更适合。

集成或组合回归方法。一旦确认了几个方法，但又不确定如何取舍时，可以将多个回归模型做集成或组合方法使用，即同时将多个模型的结果通过加权、均值等方式确定最终输出结果值。

4.2.6　代码实操：Python回归分析

sklearn中的回归有多种方法，广义线性回归集中在linear_model库下，例如普通线性回归、Lasso、岭回归等；另外还有其他非线性回归方法，例如核svm、集成方法、贝叶斯回归、K近邻回归、决策树回归等，这些不同回归算法分布在不同的库中。

本示例模拟的是针对一批训练集做多个回归模型的训练和评估，从中选择效果较好的模型并对新数据集做回归预测。本示例主要使用Sklearn的多个回归算法做回归分析、用Matplotlib做图形展示。数据源文件regression.txt位于“附件-chapter4”中，默认工作目录为“附件-chapter4”（如果不是，请cd切换到该目录下，否则会报“IOError:File regression.txt does not exist”）。完整代码如下：

# 导入库
import numpy as np  # numpy库
from sklearn.linear_model import BayesianRidge, LinearRegression, ElasticNet  # 批量导入要实现的回归算法
from sklearn.svm import SVR  # SVM中的回归算法
from sklearn.ensemble.gradient_boosting import GradientBoostingRegressor  # 集成算法
from sklearn.model_selection import cross_val_score  # 交叉检验
from sklearn.metrics import explained_variance_score, mean_absolute_error, mean_squared_error, r2_score  # 批量导入指标算法
import pandas as pd  # 导入pandas
import matplotlib.pyplot as plt  # 导入图形展示库

# 数据准备
raw_data = np.loadtxt('regression.txt')  # 读取数据文件
X = raw_data[:, :-1]  # 分割自变量
y = raw_data[:, -1]  # 分割因变量

# 训练回归模型
n_folds = 6  # 设置交叉检验的次数
model_br = BayesianRidge()  # 建立贝叶斯岭回归模型对象
model_lr = LinearRegression()  # 建立普通线性回归模型对象
model_etc = ElasticNet()  # 建立弹性网络回归模型对象
model_svr = SVR()  # 建立支持向量机回归模型对象
model_gbr = GradientBoostingRegressor()  # 建立梯度增强回归模型对象
model_names = ['BayesianRidge', 'LinearRegression', 'ElasticNet', 'SVR', 'GBR']  # 不同模型的名称列表
model_dic = [model_br, model_lr, model_etc, model_svr, model_gbr]  # 不同回归模型对象的集合
cv_score_list = []  # 交叉检验结果列表
pre_y_list = []  # 各个回归模型预测的y值列表
for model in model_dic:  # 读出每个回归模型对象
    scores = cross_val_score(model, X, y, cv=n_folds)  # 将每个回归模型导入交叉检验模型中做训练检验
    cv_score_list.append(scores)  # 将交叉检验结果存入结果列表
    pre_y_list.append(model.fit(X, y).predict(X))  # 将回归训练中得到的预测y存入列表

# 模型效果指标评估
n_samples, n_features = X.shape  # 总样本量,总特征数
model_metrics_name = [explained_variance_score, mean_absolute_error, mean_squared_error, r2_score]  # 回归评估指标对象集
model_metrics_list = []  # 回归评估指标列表
for i in range(5):  # 循环每个模型索引
    tmp_list = []  # 每个内循环的临时结果列表
    for m in model_metrics_name:  # 循环每个指标对象
        tmp_score = m(y, pre_y_list[i])  # 计算每个回归指标结果
        tmp_list.append(tmp_score)  # 将结果存入每个内循环的临时结果列表
    model_metrics_list.append(tmp_list)  # 将结果存入回归评估指标列表
df1 = pd.DataFrame(cv_score_list, index=model_names)  # 建立交叉检验的数据框
df2 = pd.DataFrame(model_metrics_list, index=model_names, columns=['ev', 'mae', 'mse', 'r2'])  # 建立回归指标的数据框
print ('samples: %d \t features: %d' % (n_samples, n_features))  # 打印输出样本量和特征数量
print (70 * '-')  # 打印分隔线
print ('cross validation result:')  # 打印输出标题
print (df1)  # 打印输出交叉检验的数据框
print (70 * '-')  # 打印分隔线
print ('regression metrics:')  # 打印输出标题
print (df2)  # 打印输出回归指标的数据框
print (70 * '-')  # 打印分隔线
print ('short name \t full name')  # 打印输出缩写和全名标题
print ('ev \t explained_variance')
print ('mae \t mean_absolute_error')
print ('mse \t mean_squared_error')
print ('r2 \t r2')
print (70 * '-')  # 打印分隔线

# 模型效果可视化
plt.figure()  # 创建画布
plt.plot(np.arange(X.shape[0]), y, color='k', label='true y')  # 画出原始值的曲线
color_list = ['r', 'b', 'g', 'y', 'c']  # 颜色列表
linestyle_list = ['-', '.', 'o', 'v', '*']  # 样式列表
for i, pre_y in enumerate(pre_y_list):  # 读出通过回归模型预测得到的索引及结果
    plt.plot(np.arange(X.shape[0]), pre_y_list[i], color_list[i], label=model_names[i])  # 画出每条预测结果线
plt.title('regression result comparison')  # 标题
plt.legend(loc='upper right')  # 图例位置
plt.ylabel('real and predicted value')  # y轴标题
plt.show()  # 展示图像

# 模型应用
print ('regression prediction')
new_point_set = [[1.05393, 0., 8.14, 0., 0.538, 5.935, 29.3, 4.4986, 4., 307., 21., 386.85, 6.58],
                 [0.7842, 0., 8.14, 0., 0.538, 5.99, 81.7, 4.2579, 4., 307., 21., 386.75, 14.67],
                 [0.80271, 0., 8.14, 0., 0.538, 5.456, 36.6, 3.7965, 4., 307., 21., 288.99, 11.69],
                 [0.7258, 0., 8.14, 0., 0.538, 5.727, 69.5, 3.7965, 4., 307., 21., 390.95, 11.28]]  # 要预测的新数据集
for i, new_point in enumerate(new_point_set):  # 循环读出每个要预测的数据点
    new_pre_y = model_gbr.predict(new_point)  # 使用GBR进行预测
    print ('predict for new point %d is:  %.2f' % (i + 1, new_pre_y))  # 打印输出每个数据点的预测信息

上述代码以空行分为6个部分。

第一部分导入库。本示例中用到的库比较多：Numpy用于文件读写和基本数据处理；sklearn.linear_model中的三个广义线性回归下的回归算法BayesianRidge、Linear-Regre-ssion、ElasticNet，sklearn.svm中的svr回归算法，sklearn.ensemble.gradient_boosting中的GradientBoostingRegressor回归算法；cross_val_score用来做交叉检验；sklearn.metrics下的多个回归模型指标评估方法用来做模型评估；Pandas用来做格式化的数据框，便于数据结果输出；Matplotlib用来做预测结果的展示和对比。这里导入多种回归算法，目的是检验在没有先验经验的条件下，通过多个模型的训练，从中找到最佳拟合算法做后续应用。

第二部分数据准备。使用Numpy的loadtxt方法读取数据文件，并使用索引切片功能从数据集中分割出自变量X和因变量y。自变量X拥有506个样本，13个特征变量。

第三部分训练回归模型。本过程使用交叉检验的方法评估不同模型的训练效果。

先设置交叉检验的次数为6（6折交叉检验），后续会在交叉检验模型训练中用到；

接着分别建立贝叶斯岭回归（BayesianRidge）、普通线性回归（LinearRegression）、弹性网络回归（ElasticNet）、支持向量机回归（SVR）、梯度增强回归（GradientBoost-ingRegressor）模型对象，前3个算法属于广义线性回归，后两个属于支持向量机和梯度增强算法的变体；

然后分别建立不同模型的名称列表、回归模型对象的集合、交叉检验结果空列表、回归模型预测的y值空列表，分别用于后续名称读取、模型代入交叉检验算法、交叉检验结果数据存储和回归预测的结果存储。

最后通过for循环读出每个回归模型对象，对每个回归模型导入交叉检验模型cross_val_score中做训练检验并设置检验6次，并将交叉检验结果通过append方法存入结果列表，将回归训练中得到的预测y通过append方法存入列表。

相关知识点：cross_val_score

cross_val_score是sklearn.model_selection中的交叉检验工具，可以对特定算法模型进行交叉检验。常用参数：

sklearn.model_selection.cross_val_score(estimator, X, y=None, groups=None, scoring= None, cv=None, n_jobs=1, verbose=0, fit_params=None, pre_dispatch='2*n_jobs')

estimator：要应用交叉检验的模型对象，该模型在交叉检验中会通过使用fit方法来训练模型

X，y：交叉检验的数据集X和目标y，其中X至少是2维数据

cv：要进行交叉检验的模式。如果cv值为空，那么默认使用3折交叉检验；如果cv值是整数，那么使用按照指定的数量做交叉检验。这两种情况下都会调用sklearn.model_selection.StratifiedKFold方法实现。如果cv值是用来交叉检验的生成器对象或可迭代的测试和训练集，那么将使用sklearn.model_selection.KFold方法。

cross_val_score交叉检验后返回的得分默认调用算法模型的score方法做得分估计，因此不同的算法模型其得分计算方法可能有差异，具体取决于模型本身score方法（scorer（estimator，X，y））的计算逻辑。当然，也可以通过设置scoring参数字符串来手动指定得分计算方法，例如cross_val_score（estimator，X，y，scor-ing='r2'，cv=6）。常用scoring方法如表4-3所示。

表4-3　手动指定的scoring方法

第四部分模型效果指标评估。

先通过自变量数据集的shape获得总样本量和总特征数；

接着创建一个列表用于存储接下来用到的回归模型评估指标对象，explained_variance_score，mean_absolute_error，mean_squared_error，r2_score，这些对象在第一步导入库时已经导入，因此无需做其他处理便可直接使用；创建用于存储不同回归评估指标的列表。

explained_variance_score：解释回归模型的方差得分，其值取值范围是[0，1]，越接近于1说明自变量越能解释因变量的方差变化，值越小则说明效果越差。

mean_absolute_error：平均绝对误差（Mean Absolute Error，MAE），用于评估预测结果和真实数据集的接近程度的程度，其值越小说明拟合效果越好。

mean_squared_error：均方差（Mean squared error，MSE），该指标计算的是拟合数据和原始数据对应样本点的误差的平方和的均值，其值越小说明拟合效果越好。

r2_score：判定系数，其含义是也是解释回归模型的方差得分，其值取值范围是[0，1]，越接近于1说明自变量越能解释因变量的方差变化，值越小则说明效果越差。

然后通过两层for循环来分别对每个模型和每种回归模型评估方法进行计算。需要注意的是里面有一个用于存储每个内循环的临时结果列表，需要在内循环结束之后，每次将其结果存储到外部循环中的回归评估指标列表。

下面是通过Pandas的Dataframe方法分别建立交叉检验和回归指标数据框，然后打印输出样本量和特征数量、交叉检验和回归指标数据框，并通过打印横线来做逻辑输出切分；为了更好的理解指标的含义，打印输出缩写和全名标题对照表。在打印过程中使用的方法跟聚类中的类似，但是对于格式化的输出，这里使用了Pandas的数据框，这样如果数据量比较多的情况下，更容易控制。输出结果如下：

samples: 506       features: 13
--------------------------------------------------------------------------
cross validation result:
                         0         1         2         3         4         5
BayesianRidge     0.662422  0.677079  0.549702  0.776896 -0.139738 -0.024448
LinearRegression  0.642240  0.611521  0.514471  0.785033 -0.143673 -0.015390
ElasticNet        0.582476  0.603773  0.365912  0.625645  0.437122  0.200454
SVR              -0.000799 -0.004447 -1.224386 -0.663773 -0.122252 -1.374062
GBR               0.747920  0.809555  0.767936  0.862972  0.375267  0.559758
--------------------------------------------------------------------------
regression metrics:
                        ev       mae        mse        r2
BayesianRidge     0.731143  3.319204  22.696772  0.731143
LinearRegression  0.740608  3.272945  21.897779  0.740608
ElasticNet        0.686094  3.592915  26.499828  0.686094
SVR               0.173548  5.447960  71.637552  0.151410
GBR               0.975126  1.151773   2.099835  0.975126
--------------------------------------------------------------------------
short name       full name
ev       explained_variance
mae       mean_absolute_error
mse       mean_squared_error
r2       r2

通过结果可以看出，增强梯度（GBR）回归是所有模型中拟合效果最好的，表现在能解释97%的方差变化，并且各个误差项的值都是最低的。另外，还有一个重要因素是，梯度增强算法在测试的6次中，其结果相对稳定性较高，这也说明了该算法在应对不同数据集的稳定效果。指标中r2跟ev的值是相同的，原因是r2作为判定系数，本身就是解释方差的比例，含义相同。

第五部分模型效果可视化。

该部分使用Matplotlib将原始数据y和其他5个回归模型预测到的y值通过趋势图进行对比。

先通过figure方法创建画布，然后画出原始数据y的分布，其中画图所需的x值域，使用自变量集X的形状得到一个自增数字列表。

接着画出其他5个回归预测结果展示所需要的颜色列表、样式列表，并通过for集合plot方法循环画出每个预测结果趋势线。

最后设置图形标题、图例位置（右上角）、y轴标题并展示图像，如图4-4所示。

不要忘记设置图像位置，否则即使在plot方法中设置了label值（图例标签值），也会由于缺少图例位置而无法显示图例。

图4-4　不同回归模型预测结果对比

相关知识点：Matplotlib基本样式控制

Matplotlib中对于样式的基本控制如下：

title：图形标题，通过字符串命名，用法为plt.title（'regression result comparison'）

legend：设置图例，该方法有众多参数可用，经常使用的是其位置控制参数loc，通过该参数可控制图例显示的位置，值可以直接写定位字符串，也可以写定位代码。例如plt.legend（loc='upper right'）等价于plt.legend（loc=1）。常用location位置控制值如下：

表4-4　location位置控制值

xlabel：x轴标题，通过字符串命名，用法为plt.xlabel（'real and predicted value'）

ylabel：y轴标题，通过字符串命名，用法为plt.ylabel（'real and predicted value'）

xtickets：x轴刻度线标签，通过字符串与x轴刻度对应，用法为：plt.xticks（np.arange（5），（'Tom'，'Dick'，'Harry'，'Sally'，'Sue'））

ytickets：y轴刻度线标签，通过字符串与y轴刻度对应，用法为：plt.yticks（np.arange（5），（'Tom'，'Dick'，'Harry'，'Sally'，'Sue'））

xlim：设置x轴的范围，通过设置最大值最小值范围实现，用法为plt.xlim（（1，10））

ylim：设置y轴的范围，通过设置最大值最小值范围实现，用法为plt.ylim（（1，10））

text：增加文字到图形中，用法为plt.text（x，y，string），其中x和y为要添加文字的坐标，string为文字字符串

suptitle：图像居中主标题，通过字符串命名，用法为plt.suptitle（'classi-fication result'）

axhline：在图形中增加一条水平线，用法为plt.axhline（y=0，xmin=0，xmax=1），其中y=0控制水平线的垂直位置，xmin和xmax控制水平线长度

axvline：在图形中增加一条垂直线，用法为plt.axvline（x=0，ymin=0，ymax=1，其中x=0控制水平线的垂直位置，ymin和ymax控制水平线长度

样式控制小技巧：

关闭坐标轴刻度（不显示刻度）：plt.xticks（[]）、plt.yticks（[]）

关闭坐标轴：plt.axis（'off'）

通过legend可批量为不同的图形对象增加图例，例如plt.legend（[s1，s2]，['label1'，'label2']，loc=0），其中s1和s2为不同的图形对象，例如线图、散点图等。

从图中看出，梯度增强算法（青色）的拟合程度跟原始数据（黑色）最接近，这也印证了在上面指标中输出的结果。

第六部分模型应用。基于新给出的包含多个数据点的数据集，我们要预测未来值，直接使用for循环读出每个数据，然后使用最优回归模型梯度增强的predict方法做预测，得到如下结果：

predict for new point 1 is:  21.49
predict for new point 2 is:  16.84
predict for new point 3 is:  19.50
predict for new point 4 is:  19.16

上述过程中，主要需要考虑的关键点是：如何在样本维度较多的情况下，选择最佳的回归模型，这里我们使用了多个模型，并结合交叉检验的方法找到最优的回归方法。

代码实操小结：本节的代码看似较多，但大多数方法之前已经介绍过，主要应用了以下几个知识点：

通过Numpy的loadtxt读取数据文件。

对矩阵做切片操作。

使用sklearn.linear_model的BayesianRidge、LinearRegression、ElasticNet做广义线性回归分析。

使用sklearn.svm中的svr、sklearn.ensemble.gradient_boosting中的GradientBoosting-Regressor做非线性回归分析。

使用模型对象的fit方法训练模型，predict方法做预测，并可以将fit和predict一起使用（即model.fit（X_train）.predict（X_test））；通过列表将模型对象本身做集合以便于迭代循环处理。

使用sklearn.metrics中的explained_variance_score，mean_absolute_error，mean_squared_error，r2_score做回归模型效果评估。

对列表通过append追加元素。

通过矩阵的shape方法获得矩阵形状。

通过Pandas的DataFrame创建新的数据框，并指定列名、索引名。

通过enumerate方法结合for循环，读出目标可迭代对象的索引和迭代值，enume-rate方法在需要使用对象索引值非常有效。

使用Numpy的arange方法创建一个可迭代的连续区间。

使用matplotlib.pyplot的plot方法展示图形，并设置不同的展示样式，包括颜色、样式、图例等，针对展示图像设置标题、图例位置、y轴标题等。

使用print方法结合占位符输出特定字符串或变量，使用不同的占位符类型（d/f等）并指定宽度（使用\t做tab分隔）做打印数据格式化。