线性回归与岭回归以及LASSO回归结果比较

为了能够比较直观地了解到三种回归方法的区别，本文基于李子奈、潘文卿的《计量经济学》（第四版）第四章的同一案列进行三种方法的结果比较，具体推导过程可参见文后提供的参考链接。

首先导入需要用到的库

import warnings
import numpy as np
import pandas as pd
import seaborn as sns
import statsmodels.api as sm
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge,RidgeCV
from sklearn.linear_model import Lasso,LassoCV
from sklearn.model_selection import train_test_split
from statsmodels.stats.outliers_influence import variance_inflation_factor

导入需要使用的数据，此数据可在文后提供的链接进行免费下载，其中包含了《计量经济学》（第四版）中所有的数据。

data = pd.read_excel("./计量经济学数据.xlsx",sheet_name="Sheet6")
data.head()

数据预处理，对数据进行去对数处理

# 数据全部取对数
log_data = np.log10(data.iloc[:,1:])
log_data.head()

多重共线性检验

计算变量间的相关性

# 计算相关系数
corr = log_data.corr()
corr

绘制变量相关性热力图

plt.figure(figsize=(12,8),dpi=80)
# ["cyan","orange","green","blue","pink","black","red"]
sns.heatmap(corr,cmap="rainbow")
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
plt.title("变量相关性热力图")

通过以上图表可知，数据集中的变量间存在较高的相关性，下面再计算变量的方差膨胀因子。

计算方差膨胀因子

x = log_data.iloc[:,1:]
y = log_data.iloc[:,0]
col = list(range(len(x.columns)))
# 所有自变量的值与每个自变量逐个计算方差膨胀因子
vif =  [variance_inflation_factor(x.iloc[:,col].values, ix) for ix in range(len(x.columns))]
vif

线性回归（linear regression）

进行线性回归拟合

linear = sm.OLS(y,x).fit()
linear.summary()

                                    OLS Regression Results                            
======================================================================================
Dep. Variable:           粮食产业Y(万吨)  R-squared:                       0.985
Model:                            OLS   Adj. R-squared:                  0.981
Method:                 Least Squares   F-statistic:                     262.3
Date:                Sun, 10 May 2020   Prob (F-statistic):           1.17e-20
Time:                        21:13:27   Log-Likelihood:                 40.002
No. Observations:                  31   AIC:                            -66.00
Df Residuals:                      24   BIC:                            -55.97
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
=======================================================================================
                         coef      std err      t         P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------------
const                   -0.3332     0.160     -2.088      0.048      -0.662      -0.004
粮食播种面积X1(千公顷)     0.7569      0.092      8.196      0.000       0.566       0.947
有效灌溉面积X2（千公顷）    0.2463      0.097      2.529      0.018       0.045       0.447
化肥施用量X3（万吨）       0.0002      0.108      0.002      0.999      -0.223       0.224
大型拖拉机X4（千台）       0.0298      0.032      0.918      0.368      -0.037       0.097
小型拖拉机X5（千台）       -0.0325     0.034     -0.958      0.348      -0.102       0.037
农用排灌柴油机X6（千台）    0.0508      0.042      1.222      0.234      -0.035       0.137
========================================================================================
Omnibus:                        7.477   Durbin-Watson:                  1.275
Prob(Omnibus):                  0.024   Jarque-Bera (JB):               6.105
Skew:                           0.784   Prob(JB):                       0.0472
Kurtosis:                       4.507   Cond. No.                       89.8
========================================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

计算均方根误差

from sklearn.metrics import mean_squared_error
linear_predict = linear.predict(x)
RMSE = np.sqrt(mean_squared_error(y,linear_predict))
RMSE

由以上计算结果可知，调整后的$R^2$为0.981，同时有三个变量的方差膨胀因子都大于50，即存在多重共线性。

岭回归（ridge regression）

自变量数据重构

# 去除上一个模型加入的常量
x = x.iloc[:,1:]
# 选择训练集和测试集
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3,random_state=42)
alphas = np.linspace(-1,5,80)

绘制回归系数收敛曲线

ridge_coefs = []
for alpha in alphas:
    linear_R = Ridge(alpha=alpha, normalize=True)
    linear_R.fit(x_train,y_train)
    ridge_coefs.append(linear_R.coef_)
# 先设置画布大小，再画图
plt.figure(figsize=(12,8),dpi=80)   
plt.plot(alphas,ridge_coefs)
plt.title("回归系数收敛曲线")
plt.xlabel("lamdas")
plt.ylabel("coefficient")
plt.show()

进行10重交叉验证，确定最优lamda值

# 岭回归模型的交叉验证
# 设置交叉验证的参数，对于每一个lamda值，都执行10重交叉验证
ridge_cv = RidgeCV(alphas = alphas, normalize=True, scoring='neg_mean_squared_error', cv = 10)
ridge_cv.fit(x_test,y_test)
best_alpha = ridge_cv.alpha_
# 使用最优的lamda值进行拟合
linear_R = Ridge(alpha=best_alpha,normalize=True)
linear_R.fit(x_train,y_train)
ridge_predict = linear_R.predict(x_test)
linear_R.score(x_test,y_test)

计算均方根误差

# 预测效果验证
from sklearn.metrics import mean_squared_error
RMSE = np.sqrt(mean_squared_error(y_test,ridge_predict))
RMSE

LASSO回归

绘制回归系数收敛曲线

alphas = np.linspace(-1,1,80)
lasso_coefs = []
for alpha in alphas:
    linear_L = Lasso(alpha=alpha, normalize=True)
    linear_L.fit(x_train,y_train)
    lasso_coefs.append(linear_L.coef_)
# 先设置画布大小，再画图
plt.figure(figsize=(8,5),dpi=100)   
plt.plot(alphas,lasso_coefs)
plt.title("回归系数收敛曲线")
plt.xlabel("lamdas")
plt.ylabel("coefficient")
plt.show()

进行10重交叉验证，确定最优lamda值

lasso_cv = LassoCV(alphas = alphas, normalize=True, cv = 10, max_iter=10000)
lasso_cv.fit(x_train,y_train)
best_alpha = lasso_cv.alpha_
# 使用最优lamda值进行拟合
linear_L = Lasso(alpha=best_alpha,normalize=True,max_iter=10000)
linear_L.fit(x_train, y_train)
score = linear_L.score(x_test,y_test)
score

计算均方根误差

RMSE = np.sqrt(mean_squared_error(y_test,lasso_predict))
RMSE

综合比较

从本案例的三个拟合结果来看，首先，就拟合参数的截距项而言，从线性回归到岭回归再到LASSO回归，参数值在不断减小，其中LASSO的截距项最小，尽可能的把被解释变量都都交由自变量来解释。其次，再看自变量的，线性回归和岭回归的$X5$的参数都小于零，虽然保留了全部自变量的信息，但是在此案例中不符合实际经济意义，而LASSO回归则对自变量的数量进行了压缩，只保留了两个自变量。
从本案例的三个拟合效果来看，虽然线性回归的$R^2$最高，但是由于此案例的数据集存在多重共线性，因此极有可能存在过拟合的情况。虽然LASSO回归的$R^2$最低，同时均方根误差$(RMSE)$也相对较高，但是相对来说拟合效果还是比较好的，虽然只保留了两个解释变量，损失了一部分解释变量的信息，但是它是比较符合经济意义的。
所以就此案例拟合的三个模型来看，LASSO还是比较可取的。

数据集：
【1】计量经济学数据

参考资料：
【1】机器学习中的正则化
【2】机器学习十大经典算法之岭回归和LASSO回归（学习笔记整理）