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

往期回顾

从零开始学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 pdimport numpy as npimport matplotlib.pyplot as pltfrom sklearn.cross_validation import train_test_splitfrom sklearn.linear_model import Ridge,RidgeCVfrom sklearn.linear_model import Lasso,LassoCVfrom 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

结语

---

      OK，今天关于岭回归和LASSO回归的实战部分就介绍到这里，希望对数据挖掘或机器学习的朋友有所帮助，同时，也希望读者能够静下心来好好的复现一遍。如果你有任何问题，欢迎在公众号的留言区域表达你的疑问。同时，也欢迎各位朋友继续转发与分享文中的内容，让更多的朋友学习和进步。

关注“每天进步一点点2015”，与小编同进步！

• 文内代码及数据

链接: https://pan.baidu.com/s/1qYOGuc8 密码: xb3r