多元线性回归之岭回归模型

前言

首先思考一个问题，为什么我们前面讲了一元线性回归及多元线性回归模型，还需要岭回归模型，岭回归模型具体解决什么问题，具体哪一类情况？

先说答案：为解决多元线性回归模型中可能存在的不可逆问题，统计学家提出了岭回归模型。该模型解决问题的思路就是在线性回归模型的目标函数之上增加l2正则项（也称为惩罚项）。

线性回归模型的短板

根据线性回归模型的参数估计公式$\beta =(X^{\prime }X{)}^{-1}{X}^{\prime }y$可知，得到$\beta$的前提是矩阵X’X可逆，但在实际应用中，可能会出现自变量个数多于样本量或者自变量间存在多重共线性的情况，即X’X的行列式为0。此时将无法根据公式计算回归系数的估计值$\beta$。

当列数比行数多

构造矩阵：

$\left[\begin{array}{ccc}1& 2& 5\\ 6& 1& 3\end{array}\right]$

计算X’X
$X^{\prime }X=\left[\begin{array}{cc}1& 6\\ 2& 1\\ 5& 3\end{array}\right]\left[\begin{array}{ccc}1& 2& 5\\ 6& 1& 3\end{array}\right]=\left[\begin{array}{ccc}37& 8& 23\\ 8& 5& 13\\ 23& 13& 34\end{array}\right]$

计算行列式：
$|X^{\prime }X|=37\ast 5\ast 34+8\ast 13\ast 23+23\ast 8\ast 13-37\ast 13\ast 13-8\ast 8\ast 34-23\ast 5\ast 23=0$

当列之间存在多重共线性

构造矩阵：

$\left[\begin{array}{ccc}1& 2& 2\\ 2& 5& 4\\ 2& 3& 4\end{array}\right]$

计算X’X
$X^{\prime }X=\left[\begin{array}{ccc}1& 2& 2\\ 2& 5& 3\\ 2& 4& 4\end{array}\right]\left[\begin{array}{ccc}1& 2& 2\\ 2& 5& 4\\ 2& 3& 4\end{array}\right]=\left[\begin{array}{ccc}9& 18& 18\\ 18& 38& 36\\ 18& 36& 36\end{array}\right]$

计算行列式：
$|X^{\prime }X|=9\ast 38\ast 36+18\ast 36\ast 18+18\ast 18\ast 36-9\ast 36\ast 36-18\ast 18\ast 36-18\ast 38\ast 18=0$

岭回归模型

$$J(\beta )=\sum (y-X\beta {)}^2+\lambda {\Vert \begin{array}{c}\beta \end{array}\Vert }_2^2=\sum (y-X\beta {)}^2+\sum \lambda {\beta }^2$$

在线性回归模型的目标函数之上添加l2正则项，其中$\lambda$为非负数；
当$\lambda =0$时，目标函数退化为线性回归模型的目标函数;
当$\lambda ⟶+\infty$时,通过缩减回归系数使$\beta$趋近于0;
$\lambda$是l2正则项平方的系数，用于平衡模型方差（回归系数的方差）和偏差

系数求解

展开岭回归模型中的平方项
$$\begin{gathered} J(\beta )=(y-X\beta {)}^{\prime }(y-X\beta )+\lambda {\beta }^{\prime }\beta \\ =(y^{\prime }-{\beta }^{\prime }{X}^{\prime })(y-X\beta )+\lambda {\beta }^{\prime }\beta \\ =y^{\prime }y-y^{\prime }X\beta -{\beta }^{\prime }{X}^{\prime }y+{\beta }^{\prime }{X}^{\prime }X\beta +\lambda {\beta }^{\prime }\beta \end{gathered}$$

计算导函数
$$\begin{gathered} \frac{\partial J(\beta )}{\partial \beta }=0-X^{\prime }y-X^{\prime }y+2X^{\prime }X\beta +2\lambda \beta \\ =2(X^{\prime }X+\lambda I)\beta -2X^{\prime }y \end{gathered}$$

参数求解
$$\begin{gathered} 2(X^{\prime }X+\lambda I)\beta -2X^{\prime }y=0 \\ \beta =(X^{\prime }X+\lambda I{)}^{-1}{X}^{\prime }y \end{gathered}$$

凸优化的等价命题

$J(\beta )={\sum }_{}(y-X\beta {)}^2+\lambda {\Vert \begin{array}{c}\beta \end{array}\Vert }_2^2={\sum }_{}(y-X\beta )+{\sum }_{}\lambda {\beta }^2$
等价于 $argmin{\sum }_{}(y-X\beta {)}^2$并附加约束${\sum }_{}{\beta }^2<=t$

$\lambda$值的确定:交叉验证法

交叉验证的思想

首先将数据集拆分成k个样本量大体相当的数据组(如图中的第一行)，并且每个数据组与其他组都没有重叠的观测；
然后从k组数据中挑选k-1组数据用于模型的训练，剩下的一组数据用于模型的测试(如图中的第二行)；
以此类推，将会得到k种训练集和测试集。在每一种训练集和测试集下，都会对应一个模型及模型得分(如均方误差)；

实战

# 导入第三方模块
import pandas as pd
import numpy as np
from sklearn import model_selection
from sklearn.linear_model import Ridge,RidgeCV
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error

# 读取糖尿病数据集
diabetes = pd.read_excel(r'diabetes.xlsx')
# 构造自变量（剔除患者性别、年龄和因变量）
predictors = diabetes.columns[2:-1]
# 将数据集拆分为训练集和测试集
X_train, X_test, y_train, y_test = model_selection.train_test_split(diabetes[predictors], diabetes['Y'],  test_size = 0.2, random_state = 1234 )
# 构造不同的Lambda值
Lambdas = np.logspace(-5, 2, 200)
# 岭回归模型的交叉验证
# 设置交叉验证的参数，对于每一个Lambda值，都执行10重交叉验证
ridge_cv = RidgeCV(alphas = Lambdas, scoring='neg_mean_squared_error', cv = 10)
# 模型拟合
ridge_cv.fit(X_train, y_train)
# 返回最佳的lambda值
ridge_best_Lambda = ridge_cv.alpha_
print('ridge_best_Lambda:',ridge_best_Lambda)



# 基于最佳的Lambda值建模
ridge = Ridge(alpha = ridge_best_Lambda)
ridge.fit(X_train, y_train)
# 返回岭回归系数
series = pd.Series(index = ['Intercept'] + X_train.columns.tolist(),data = [ridge.intercept_] + ridge.coef_.tolist())
print(series)
# 预测
ridge_predict = ridge.predict(X_test)
# 预测效果验证
RMSE = np.sqrt(mean_squared_error(y_test,ridge_predict))
print('RMSE:',RMSE)

ridge_best_Lambda: 0.6080224261649427
Intercept   -390.214200
BMI            6.250932
Intercept   -390.214200
BMI            6.250932
BMI            6.250932
BP             0.949270
S1            -1.117024
S2             0.766199
S3             0.795961
S4             6.479146
S5            67.100942
S6             0.379425
dtype: float64
RMSE: 53.38642682397508