从零开始学Python【25】--岭回归及LASSO回归（实战部分）

发表: 2017-11-29 浏览: 8771

LASSO回归 Python 岭回归

往期回顾

从零开始学Python【24】--岭回归及LASSO回归（理论部分）

从零开始学Python【23】--线性回归诊断（第二部分）

从零开始学Python【22】--线性回归诊断（第一部分）

从零开始学Python【21】--线性回归（实战部分）

从零开始学Python【20】--线性回归（理论部分）

前言

在《从零开始学Python【24】--岭回归及LASSO回归（理论部分）》一文中我们详细介绍了关于岭回归和LASSO回归的理论知识，其实质就是在线性回归的基础上添加了2范数和1范数的惩罚项。这两个模型的关键点是找到一个合理的lambda系数，来平衡模型的方差和偏差，从而得到比较符合实际的回归系数。本期是基于之前讨论的理论部分，采用Python和R语言，完成对岭回归和LASSO回归的实战。文中所利用的数据集来自R语言ISLR包中的Hitters数据集，描述的是关于棒球运动员收入的相关信息，数据集和脚本可以之文末查看链接获得。

岭回归

原始数据集存在收入的缺失，我们不妨先把这样的观测删除，同时，数据集中含有离散变量，需要将这些变量转换为哑变量后方可建模，故第一步需要对原始数据集进行清洗。

# ===== Python3 =====

# 导入第三方包
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import Ridge,RidgeCV
from sklearn.linear_model import Lasso,LassoCV
from sklearn.metrics import mean_squared_error

# 读取数据
df = pd.read_csv('Hitters.csv')
# 哑变量处理
dummies = pd.get_dummies(df[['League', 'Division', 'NewLeague']])
# 将原始数据集与哑变量数据合并起来
mydf = df.join(dummies)
# 缺失值删除
mydf = mydf.dropna()
# 删除不必要的变量(字符串变量和各哑变量中的一个变量)
mydf = mydf.drop([ 'League', 'Division', 'NewLeague', 'League_N', 'Division_W', 'NewLeague_N'], axis = 1)

# 前5行展示

mydf.head()

上面的数据集清洗完毕，展现的是干净数据的前5行信息，下面要基于这个数据集进行建模。建模之前还需要将数据集拆分为两部分，一部分用于建模，另一部分用于模型的测试。

# 将数据集拆分成训练集和测试集
predictors = list(mydf.columns)

predictors.remove('Salary') # Salay变量为因变量，故需要排除

X_train, X_test, y_train, y_test = train_test_split(mydf[predictors],mydf['Salary'], train_size = 0.8, random_state = 1234 )

基于可视化，选择lambda参数

# 通过不确定的alphas值，生成不同的岭回归模型
alphas = 10**np.linspace(-3,3,100)

ridge_cofficients = []

for alpha in alphas:

    ridge = Ridge(alpha = alpha, normalize=True)

    ridge.fit(X_train, y_train)

    ridge_cofficients.append(ridge.coef_)

# 绘制alpha的对数与回归系数的关系
# 中文乱码和坐标轴负号的处理
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']

plt.rcParams['axes.unicode_minus'] = False
# 设置绘图风格
plt.style.use('ggplot')
plt.plot(alphas, ridge_cofficients)

plt.xscale('log')

plt.axis('tight')

plt.title('alpha系数与岭回归系数的关系')

plt.xlabel('Log Alpha')

plt.ylabel('Cofficients')

plt.show()

从上面的图形结果来看，alpha在10附近时，所有的自变量系数基本趋于稳定（但也不能完全确定是这个值）。接下来，我们采用交叉验证（CV）方法确定最佳的lambda值。

基于CV选择lambda值

# 为了找到最佳的lambda值，我们采用交叉验证方法

# 岭回归模型的交叉验证
ridge_cv = RidgeCV(alphas = alphas, normalize=True, scoring='mean_squared_error', cv = 10)

ridge_cv.fit(X_train, y_train)
# 取出最佳的lambda值ridge_best_alpha = ridge_cv.alpha_

ridge_best_alpha

不出所料，得到的lambda值确实在10附近，这里最佳的lambda值为10。下面，我们要基于这个最佳的lambda值进入岭回归模型的创建和模型验证的阶段。

# 基于最佳的lambda值建模
ridge = Ridge(alpha = ridge_best_alpha, normalize=True)

ridge.fit(X_train, y_train)
# 岭回归系数
ridge.coef_

# 预测
ridge_predict = ridge.predict(X_test)
# 预测效果验证
RMSE = np.sqrt(mean_squared_error(y_test,ridge_predict))

RMSE

经过模型的验证，得到的RMSE为319.9。接下来，我们利用同样的逻辑，对比一下LASSO回归模型的效果。

LASSO回归

基于可视化，选择lambda参数

# 通过不确定的alphas值，生成不同的LASSO回归模型
alphas = 10**np.linspace(-3,3,100)

lasso_cofficients = []

for alpha in alphas:

    lasso = Lasso(alpha = alpha, normalize=True, max_iter=10000)

    lasso.fit(X_train, y_train)

    lasso_cofficients.append(lasso.coef_)

# 绘制alpha的对数与回归系数的关系# 中文乱码和坐标轴负号的处理
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']

plt.rcParams['axes.unicode_minus'] = False
# 设置绘图风格
plt.style.use('ggplot')
plt.plot(alphas, lasso_cofficients)

plt.xscale('log')

plt.axis('tight')

plt.title('alpha系数与LASSO回归系数的关系')

plt.xlabel('Log Alpha')

plt.ylabel('Cofficients')

plt.show()

从图形结果来看，lambda值应该在1附近，此时LASSO回归的系数也基本趋于稳定（但也不能完全确定是这个值）。同样，我们利用CV方法，来寻找最佳的lambda值。

基于CV选择lambda值

# LASSO回归模型的交叉验证
lasso_cv = LassoCV(alphas = alphas, normalize=True, cv = 10, max_iter=10000)

lasso_cv.fit(X_train, y_train)
# 取出最佳的lambda值
lasso_best_alpha = lasso_cv.alpha_

lasso_best_alpha

通过CV方法得到的lambda结果是0.23，这与与我们看图得到的1这个值还是有一点差异的。下面，我们就基于交叉验证得到的最佳lambda值重新构造LASSO回归模型。

# 基于最佳的lambda值建模
lasso = Lasso(alpha = lasso_best_alpha, normalize=True, max_iter=10000)

lasso.fit(X_train, y_train)
# 岭回归系数lasso.coef_

# 预测
lasso_predict = lasso.predict(X_test)
# 预测效果验证
RMSE = np.sqrt(mean_squared_error(y_test,lasso_predict))

RMSE

对LASSO回归模型进行验证发现得到的RMSE更小，说明LASSO回归模型的拟合效果会更贴近于Hitters数据集的原貌。

上面的内容是基于Python工具对岭回归模型和LASSO回归模型的实战，接下来，我们再利用R语言对上面的过程再作一次复现，希望对R语言感兴趣的朋友能够有帮助。

R语言对比

由于上面的逻辑我们已经通过Python进行了一一说明，这里就不再赘述，只给出R语言代码仅供参考。

数据清洗

# 加载第三方包
library(caret)
library(glmnet)
library(ISLR)

# 哑变量处理
dummies <- dummyVars(~League+Division+NewLeague, data = Hitters)

dummies <- predict(dummies, newdata = Hitters)
# 数据合并
Hitters_dummy <- cbind(Hitters, dummies)

# 删除缺失值
Hitters_dummy <- na.omit(Hitters_dummy)
# 删除不必要的变量
Hitters_dummy <- subset(Hitters_dummy, 

                       select = -c(League,Division,NewLeague,League.N,Division.W,NewLeague.N))

head(Hitters_dummy)

岭回归—可视化选lambda

# 构建训练集和测试集set.seed(1)

index <- sample(1:nrow(Hitters_dummy), size = 0.8*nrow(Hitters_dummy))

train <- Hitters_dummy[index,]

test <- Hitters_dummy[-index,] 



# 绘制lambda值与岭回归系数的关系
fit_ridge <-  glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 0)

plot(fit_ridge,xvar = 'lambda',label = T)

岭回归—交叉验证选lambda

# 岭回归的交叉验证，确定最佳的lambda值
fit_ridge_cv <- cv.glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 0)

best_lambda_ridge <- fit_ridge_cv$lambda.min

best_lambda_ridge

岭回归—建模与验证

# 根据最佳lambda构建岭回归模型fit_ridge <- glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 0)coeff_ridge <- predict(fit_ridge, s = best_lambda_ridge, type = 'coefficients')coeff_ridge



# 模型评估pred_ridge <- predict(fit_ridge, s = best_lambda_ridge, newx = as.matrix(test[,-17]))RMSE <- sqrt(mean((test$Salary-pred_ridge)**2))RMSE

LASSO—可视化选lambda

# 绘制lambda值与LASSO回归系数的关系
fit_lasso <-  glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 1)

plot(fit_lasso,xvar = 'lambda',label = T)

LASSO—交叉验证选lambda

# 岭回归的交叉验证，确定最佳的lambda值
fit_lasso_cv <- cv.glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 1)

best_lambda_lasso <- fit_lasso_cv$lambda.min

best_lambda_lasso

LASSO—建模与验证

# 根据最佳lambda构建岭回归模型
fit_lasso <- glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 1)

coeff_lasso <- predict(fit_lasso, s = best_lambda_lasso, type = 'coefficients')

coeff_lasso

# 模型评估
pred_lasso <- predict(fit_lasso, s = best_lambda_lasso, newx = as.matrix(test[,-17]))

RMSE <- sqrt(mean((test$Salary-pred_lasso)**2))

RMSE