岭回归，Lasso回归及弹性网络回归原理及R实现

在线性大型数据集包含的变量数量大于样本数量的情况下，标准线性模型（或普通的最小二乘法）的效果较差。

更好的替代方法是``惩罚回归''，可以通过在方程中添加约束来创建线性惩罚模型，该模型因变量过多而受到惩罚（James et al. 2014,P. Bruce and Bruce (2017)）。这也称为收缩或调整方法。

施加此惩罚的结果是将系数值减小（即缩小）为零。这允许贡献较小的变量具有接近零或等于零的系数。

注意，收缩要求选择确定收缩量的调整参数（λ）。

在本章中，我们将介绍最常用的惩罚回归方法，包括岭回归，套索回归和弹性净回归。我们还将在R中提供实用示例。

Table of Contents

收缩方法

岭回归

Lasso回归

弹性网络回归

加载R包

准备数据

计算惩戒线性回归

准备数据

R函数

计算岭回归

计算lasso回归

计算弹性网络回归

模型比较

收缩方法

岭回归

岭回归缩小了回归系数，因此对结果影响较小的变量的系数接近零。系数的缩小是通过使用称为L2-正则化的惩罚项对回归模型进行惩罚来实现的，该惩罚项是系数平方的总和。

可以使用称为lambda（λ）的常量微调惩罚量。为λ选择一个合适的值至关重要。当λ=0时，惩罚项无效，岭回归将产生经典的最小二乘系数。但是，随着λ增大到无穷大，收缩罚分的影响增大，并且岭回归系数将接近于零。

请注意，与普通的最小二乘回归相比，岭回归受预测变量的规模影响很大。因此，最好在应用岭回归之前对预测变量进行标准化（即缩放）（James et al.2014），以使所有预测变量处于相同规模。

可以使用公式x'= x / sd（x）来实现预测变量x的标准化，其中sd（x）是x的标准偏差。这样的结果是，所有标准化预测变量的标准偏差都将为一，从而使最终拟合不依赖于预测变量的度量范围。

与普通的最小二乘方法（第@ref（linear-regression）章）相比，脊回归的一个重要优点是在具有大量预测变量（p ）大于观察次数（n）。

岭回归的一个缺点是它将包含最终模型中的所有预测变量，这与逐步回归方法（第@ref（stepwise-regression））不同，后者通常会选择涉及变量减少的模型。

岭回归将系数缩小为零，但不会将其中的任何一个都精确设置为零。套索回归是克服此缺点的替代方法。

Lasso回归

LASSO代表最小绝对收缩和选择算符。它通过使用称为L1-正则化的惩罚项惩罚回归模型来将回归系数缩小为零，这是绝对系数的总和。在套索回归的情况下，惩罚具有迫使一些系数估计（对于模型的贡献很小）正好等于零的效果。这意味着，套索还可被视为执行变量选择以降低模型复杂性的子集选择方法的替代方法。

就像在岭回归中一样，为套索选择一个好的λ值至关重要。

套索回归相对于岭回归的一个明显优势是，它可以生成更简单，更易解释的模型，其中仅包含一组减少的预测变量。但是，无论是岭回归还是套索都不会普遍主导对方。

通常，套索在某些预测变量系数较大而其余预测变量系数很小的情况下可能会表现更好。当结果是许多预测变量的函数时，岭回归将表现更好，所有预测变量的系数都大致相等（James et al.2014）。

交叉验证方法可用于识别这两种技术中的哪一种在特定数据集上更好。

弹性网络回归

Elastic Net生成了一个回归模型，该模型同时受L1范数和L2范数的影响。这样的结果是有效地缩小系数（如在岭回归中）并将某些系数设置为零（如在LASSO中）。

加载R包

tidyverse -数据集操作及可视化
caret -机器学习流程
glmnet-建立惩罚方程

library(tidyverse)
library(caret)
library(glmnet)

准备数据

# Load the data
data("Boston", package = "MASS")
# Split the data into training and test set
set.seed(123)
training.samples <- Boston$medv %>%createDataPartition(p = 0.8, list = FALSE)
train.data  <- Boston[training.samples, ]
test.data <- Boston[-training.samples, ]

计算惩戒线性回归

准备数据

y 结局变量
x 预测变量。这应该使用model.matrix（）函数创建，该函数可以将任何分类变量自动转换为虚拟变量，这很重要，因为glmnet（）仅可以接受数字，定量输入。创建模型矩阵后，我们删除索引= 1处的截距分量。

# Predictor variables
x <- model.matrix(medv~., train.data)[,-1]
# Outcome variable
y <- train.data$medv

R函数

使用 glmnet()函数构建惩戒方程 [glmnet包]

glmnet(x, y, alpha = 1, lambda = NULL)

x: 预测变量矩阵
y: 二分类结局变量
alpha: 弹性网络混合系数.
- 1: Lasso回归
- 0: 岭回归
- 0 ~ 1 ：弹性网络回归
lamba: 收缩系数，需自定

在惩罚回归中，您需要指定一个常数λ来调整系数收缩量。可以将最适合您数据的lambda定义为使交叉验证预测错误率最小化的lambda。这可以使用cv.glmnet（）函数自动确定。

在以下各节中，我们首先计算岭，套索和弹性网回归模型。接下来，我们将比较不同的模型，以便为我们的数据选择最佳模型。

最好的模型定义为预测误差最低的模型RMSE。

计算岭回归

# Find the best lambda using cross-validation
set.seed(123)
cv <- cv.glmnet(x, y, alpha = 0)
# Display the best lambda value
cv$lambda.min

## [1] 0.758

# Fit the final model on the training data
model <- glmnet(x, y, alpha = 0, lambda = cv$lambda.min)
# Display regression coefficients
coef(model)

## 14 x 1 sparse Matrix of class "dgCMatrix"
##                    s0
## (Intercept)  28.69633
## crim         -0.07285
## zn            0.03417
## indus        -0.05745
## chas          2.49123
## nox         -11.09232
## rm            3.98132
## age          -0.00314
## dis          -1.19296
## rad           0.14068
## tax          -0.00610
## ptratio      -0.86400
## black         0.00937
## lstat        -0.47914

# Make predictions on the test data
x.test <- model.matrix(medv ~., test.data)[,-1]
predictions <- model %>% predict(x.test) %>% as.vector()
# Model performance metrics
data.frame(RMSE = RMSE(predictions, test.data$medv),Rsquare = R2(predictions, test.data$medv)
)

##   RMSE Rsquare
## 1 4.98   0.671

请注意，默认情况下，函数glmnet（）标准化变量，以便它们的大小可比。但是，系数始终以原始比例返回。

计算lasso回归

用于岭回归的R代码之间的唯一区别是，对于套索回归，需要指定参数alpha = 1而不是alpha = 0（对于岭回归）。

# Find the best lambda using cross-validation
set.seed(123)
cv <- cv.glmnet(x, y, alpha = 1)
# Display the best lambda value
cv$lambda.min

## [1] 0.00852

# Fit the final model on the training data
model <- glmnet(x, y, alpha = 1, lambda = cv$lambda.min)
# Dsiplay regression coefficients
coef(model)

## 14 x 1 sparse Matrix of class "dgCMatrix"
##                    s0
## (Intercept)  36.90539
## crim         -0.09222
## zn            0.04842
## indus        -0.00841
## chas          2.28624
## nox         -16.79651
## rm            3.81186
## age           .
## dis          -1.59603
## rad           0.28546
## tax          -0.01240
## ptratio      -0.95041
## black         0.00965
## lstat        -0.52880

# Make predictions on the test data
x.test <- model.matrix(medv ~., test.data)[,-1]
predictions <- model %>% predict(x.test) %>% as.vector()
# Model performance metrics
data.frame(RMSE = RMSE(predictions, test.data$medv),Rsquare = R2(predictions, test.data$medv)
)

##   RMSE Rsquare
## 1 4.99   0.671

计算弹性网络回归

我们使用caret包自动选择最佳的调整参数alpha和lambda。 caret程序包测试了一系列可能的alpha和lambda值，然后选择了lambda和alpha的最佳值，从而生成了最终模型，即弹性网模型。

在这里，我们将测试10个不同的alpha和lambda值的组合。使用选项tuneLength指定，最佳的alpha和lambda值是那些使交叉验证误差最小的值。

# Build the model using the training set
set.seed(123)
model <- train(medv ~., data = train.data, method = "glmnet",trControl = trainControl("cv", number = 10),tuneLength = 10
)
# Best tuning parameter
model$bestTune

##   alpha lambda
## 6   0.1   0.21

# Coefficient of the final model. You need
# to specify the best lambda
coef(model$finalModel, model$bestTune$lambda)

## 14 x 1 sparse Matrix of class "dgCMatrix"
##                     1
## (Intercept)  33.04083
## crim         -0.07898
## zn            0.04136
## indus        -0.03093
## chas          2.34443
## nox         -14.30442
## rm            3.90863
## age           .
## dis          -1.41783
## rad           0.20564
## tax          -0.00879
## ptratio      -0.91214
## black         0.00946
## lstat        -0.51770

# Make predictions on the test data
x.test <- model.matrix(medv ~., test.data)[,-1]
predictions <- model %>% predict(x.test)
# Model performance metrics
data.frame(RMSE = RMSE(predictions, test.data$medv),Rsquare = R2(predictions, test.data$medv)
)

##   RMSE Rsquare
## 1 4.98   0.672

模型比较

不同模型的性能指标具有可比性。使用LASSO或弹性网络回归将预测变量年龄的系数设置为零，从而导致与包含所有预测变量的岭回归相比更简单的模型。在我们的示例中，我们可以选择套索或弹性净回归模型。

请注意，我们可以使用caret操作流程轻松计算和比较岭，套索和弹性净回归。caret函数将自动选择最佳调整参数值，计算最终模型并使用交叉验证技术评估模型性能。

使用 caret 包

设置Lambda值的范围：

lambda <- 10^seq(-3, 3, length = 100)

计算岭回归

# Build the model
set.seed(123)
ridge <- train(medv ~., data = train.data, method = "glmnet",trControl = trainControl("cv", number = 10),tuneGrid = expand.grid(alpha = 0, lambda = lambda))
# Model coefficients
coef(ridge$finalModel, ridge$bestTune$lambda)
# Make predictions
predictions <- ridge %>% predict(test.data)
# Model prediction performance
data.frame(RMSE = RMSE(predictions, test.data$medv),Rsquare = R2(predictions, test.data$medv)
)

计算LASSO回归

# Build the model
set.seed(123)
lasso <- train(medv ~., data = train.data, method = "glmnet",trControl = trainControl("cv", number = 10),tuneGrid = expand.grid(alpha = 1, lambda = lambda))
# Model coefficients
coef(lasso$finalModel, lasso$bestTune$lambda)
# Make predictions
predictions <- lasso %>% predict(test.data)
# Model prediction performance
data.frame(RMSE = RMSE(predictions, test.data$medv),Rsquare = R2(predictions, test.data$medv)
)

弹性网络回归

# Build the model
set.seed(123)
elastic <- train(medv ~., data = train.data, method = "glmnet",trControl = trainControl("cv", number = 10),tuneLength = 10)
# Model coefficients
coef(elastic$finalModel, elastic$bestTune$lambda)
# Make predictions
predictions <- elastic %>% predict(test.data)
# Model prediction performance
data.frame(RMSE = RMSE(predictions, test.data$medv),Rsquare = R2(predictions, test.data$medv)
)

比较模型性能

使用caret包比较不同模型性能，最佳模型定义为最小化预测误差的模型。

models <- list(ridge = ridge, lasso = lasso, elastic = elastic)
resamples(models) %>% summary( metric = "RMSE")

##
## Call:
## summary.resamples(object = ., metric = "RMSE")
##
## Models: ridge, lasso, elastic
## Number of resamples: 10
##
## RMSE
##         Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## ridge   3.10    3.96   4.38 4.73    5.52 7.43    0
## lasso   3.16    4.03   4.39 4.73    5.51 7.27    0
## elastic 3.13    4.00   4.37 4.72    5.52 7.32    0