使用 R 进行线性、套索和岭回归

2019-11-12 08:00:00 · 飞浪

介绍许多组织使用机器学习来识别和解决业务问题。监督式机器学习算法有两种类型：分类和回归。本指南将重点介绍预测连续结果的回归模型。您将学习如何使用R实现线性和正则化回归模型。我们将讨论

介绍

许多组织使用机器学习来识别和解决业务问题。监督式机器学习算法有两种类型：分类和回归。本指南将重点介绍预测连续结果的回归模型。您将学习如何使用 R 实现线性和正则化回归模型。

我们将讨论的主题包括：

线性回归
岭回归
套索回归
弹性网络回归

让我们首先看看现实生活中的情况和数据。

数据

失业是任何国家面临的重大社会经济和政治问题，因此，管理失业是任何政府的首要任务。在本指南中，我们将构建回归算法来预测经济体中的失业率。

数据来自美国经济时间序列数据，可从https://research.stlouisfed.org/fred2获取。数据包含 574 行和 5 个变量，如下所述：

psavert：个人储蓄率
pce：个人消费支出，单位：十亿美元
uempmed：失业持续时间中位数，以周为单位
pop：总人口，单位：百万
unemploy：失业人数，以百万计。这是因变量。

评估指标

我们将使用两个指标来评估模型的性能：R 平方值和均方根误差 (RMSE)。理想情况下，较低的 RMSE 和较高的 R 平方值表明模型良好。

让我们首先加载所需的库和数据。

      library(plyr)
library(readr)
library(dplyr)
library(caret)
library(ggplot2)
library(repr)

dat <- read_csv("reg_data.csv")
glimpse(dat)
    

输出：

      Observations: 574
Variables: 5
$ pce      <dbl> 507.4, 510.5, 516.3, 512.9, 518.1, 525.8, 531.5, 534.2, 544.9, 544.6, 5...
$ pop      <dbl> 198.7, 198.9, 199.1, 199.3, 199.5, 199.7, 199.8, 199.9, 200.1, 200.2, 2...
$ psavert  <dbl> 12.5, 12.5, 11.7, 12.5, 12.5, 12.1, 11.7, 12.2, 11.6, 12.2, 12.0, 11.6,...
$ uempmed  <dbl> 4.5, 4.7, 4.6, 4.9, 4.7, 4.8, 5.1, 4.5, 4.1, 4.6, 4.4, 4.4, 4.5, 4.2, 4...
$ unemploy <dbl> 2.9, 2.9, 3.0, 3.1, 3.1, 3.0, 2.9, 3.0, 2.9, 2.7, 2.7, 2.9, 2.9, 2.8, 2...
    

输出显示数据集中的所有变量都是数值变量（标记为“dbl”）。

数据分区

我们将在训练集上构建模型，并在测试集上评估其性能。这称为用于评估模型性能的保留验证方法。

下面的第一行代码设置了随机种子，以确保结果的可重复性。第二行创建索引，用于对数据分区进行随机抽样观察。接下来的两行代码创建训练集和测试集，而最后两行打印训练集和测试集的维度。训练集包含 70% 的数据，而测试集包含剩余的 30%。

      set.seed(100) 

index = sample(1:nrow(dat), 0.7*nrow(dat)) 

train = dat[index,] # Create the training data 
test = dat[-index,] # Create the test data

dim(train)
dim(test)
    

输出：

      1] 401   5

[1] 173   5

缩放数值特征

数字特征需要缩放；否则，它们可能会对建模过程产生不利影响。下面的第一行代码创建一个包含独立数字变量名称的列表。第二行使用caret 包中的preProcess函数完成缩放任务。预处理对象仅适用于训练数据，而缩放则应用于训练集和测试集。这是在下面的第三行和第四行代码中完成的。第五行打印缩放后的训练数据集的摘要。

输出显示，现在除目标变量unemploy未被缩放外，所有数值特征的平均值都为零。

      cols = c('pce', 'pop', 'psavert', 'uempmed')

pre_proc_val <- preProcess(train[,cols], method = c("center", "scale"))

train[,cols] = predict(pre_proc_val, train[,cols])
test[,cols] = predict(pre_proc_val, test[,cols])

summary(train)

输出：

      pce               pop             psavert            uempmed           unemploy     
 Min.   :-1.2135   Min.   :-1.6049   Min.   :-1.89481   Min.   :-1.1200   Min.   : 2.700  
 1st Qu.:-0.8856   1st Qu.:-0.8606   1st Qu.:-0.76518   1st Qu.:-0.6198   1st Qu.: 6.500  
 Median :-0.2501   Median :-0.1416   Median :-0.05513   Median :-0.2447   Median : 7.500  
 Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 7.791  
 3rd Qu.: 0.7797   3rd Qu.: 0.9060   3rd Qu.: 0.81629   3rd Qu.: 0.1055   3rd Qu.: 8.600  
 Max.   : 2.1244   Max.   : 1.8047   Max.   : 2.91417   Max.   : 4.1571   Max.   :15.300
    

线性回归

最简单的回归形式是线性回归，它假设预测变量与目标变量具有线性关系。假设输入变量具有高斯分布，并且彼此不相关（称为多重共线性问题）。

线性回归方程可以表示为以下形式：y = a1x1 + a2x2 + a3x3 + ..... + anxn + b

在上面的等式中：

y是目标变量。
x1，x2，x3，... xn是特征。
a1, a2, a3, ... an是系数。
b是模型的参数。

模型中的参数a和b通过普通最小二乘法 (OLS) 进行选择。该方法通过最小化残差（实际值 - 预测值）的平方和来实现。

为了拟合线性回归模型，第一步是使用lm()函数在下面的第一行代码中实例化算法。第二行打印训练模型的摘要。

      lr = lm(unemploy ~ uempmed + psavert + pop + pce, data = train)
summary(lr)
    

输出：

      Call:
lm(formula = unemploy ~ uempmed + psavert + pop + pce, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4262 -0.7253  0.0278  0.6697  3.2753 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.79077    0.04712 165.352  < 2e-16 ***
uempmed      2.18021    0.08588  25.386  < 2e-16 ***
psavert      0.79126    0.13244   5.975 5.14e-09 ***
pop          5.95419    0.37405  15.918  < 2e-16 ***
pce         -5.31578    0.32753 -16.230  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9435 on 396 degrees of freedom
Multiple R-squared:  0.8542,	Adjusted R-squared:  0.8527 
F-statistic: 579.9 on 4 and 396 DF,  p-value: < 2.2e-16
    

上述输出中的重要度代码“***”表明所有特征都是重要的预测因子。调整后的 R 平方值为 0.8527，也是一个不错的结果。让我们进一步评估该模型。

模型评估指标

第一步是创建一个函数来计算评估指标 R 平方和 RMSE。第二步是在训练数据上预测和评估模型，而第三步是在测试数据上预测和评估模型。

      #Step 1 - create the evaluation metrics function

eval_metrics = function(model, df, predictions, target){
    resids = df[,target] - predictions
    resids2 = resids**2
    N = length(predictions)
    r2 = as.character(round(summary(model)$r.squared, 2))
    adj_r2 = as.character(round(summary(model)$adj.r.squared, 2))
    print(adj_r2) #Adjusted R-squared
    print(as.character(round(sqrt(sum(resids2)/N), 2))) #RMSE
}

# Step 2 - predicting and evaluating the model on train data
predictions = predict(lr, newdata = train)
eval_metrics(lr, train, predictions, target = 'unemploy')

# Step 3 - predicting and evaluating the model on test data
predictions = predict(lr, newdata = test)
eval_metrics(lr, test, predictions, target = 'unemploy')
    

输出：

      1] "0.85"
      
[1] "0.94"
     
[1] "0.85"
     
[1] "1.1"

上面的输出显示，两个评估指标之一的 RMSE 对于训练数据为 0.94 百万，对于测试数据为 1.1 百万。另一方面，训练和测试数据的 R 平方值都在 85% 左右，这表明性能良好。

正则化

线性回归算法的工作原理是为每个独立变量选择系数，以最小化损失函数。但是，如果系数很大，则可能导致训练数据集的过度拟合，并且这种模型无法很好地推广到看不见的测试数据。为了克服这个缺点，我们将进行正则化，这会惩罚较大的系数。本指南的以下部分将讨论各种正则化算法。

我们将使用glmnet ()包来构建正则化回归模型。glmnet函数不适用于数据框，因此我们需要为训练特征创建一个数字矩阵和一个目标值向量。

下面的代码行使用caret 包中的dummyVars函数执行创建模型矩阵的任务。然后应用predict函数创建用于训练和测试的数字模型矩阵。

      cols_reg = c('pce', 'pop', 'psavert', 'uempmed', 'unemploy')

dummies <- dummyVars(unemploy ~ ., data = dat[,cols_reg])

train_dummies = predict(dummies, newdata = train[,cols_reg])

test_dummies = predict(dummies, newdata = test[,cols_reg])

print(dim(train_dummies)); print(dim(test_dummies))

输出：

      1] 401   4
[1] 173   4
    

岭回归

岭回归是线性回归的扩展，其中损失函数被修改以最小化模型的复杂性。此修改通过添加等于系数幅度平方的惩罚参数来完成。

损失函数 = OLS + alpha * 总和（系数值的平方）

岭回归也称为l2 正则化。下面的代码行构建了一个岭回归模型。第一行加载库，而接下来的两行创建独立变量 (x) 和因变量 (y) 的训练数据矩阵。

第四行和第五行代码对测试数据集重复了同样的步骤。第六行创建了一个 lambda 值列表，供模型尝试，而第七行构建了岭回归模型。

该模型中使用的参数包括：

nlambda：确定要测试的正则化参数的数量。
alpha：确定要使用的权重。在岭回归的情况下，alpha 的值为零。
family：确定要使用的分布系列。由于这是一个回归模型，我们将使用高斯分布。
lambda：确定要尝试的 lambda 值。

最后一行代码打印模型信息。

      library(glmnet)

x = as.matrix(train_dummies)
y_train = train$unemploy

x_test = as.matrix(test_dummies)
y_test = test$unemploy

lambdas <- 10^seq(2, -3, by = -.1)
ridge_reg = glmnet(x, y_train, nlambda = 25, alpha = 0, family = 'gaussian', lambda = lambdas)

summary(ridge_reg)
    

输出：

      Length Class     Mode   
a0         51    -none-    numeric
beta      204    dgCMatrix S4     
df         51    -none-    numeric
dim         2    -none-    numeric
lambda     51    -none-    numeric
dev.ratio  51    -none-    numeric
nulldev     1    -none-    numeric
npasses     1    -none-    numeric
jerr        1    -none-    numeric
offset      1    -none-    logical
call        7    -none-    call   
nobs        1    -none-    numeric
    

线性回归模型和正则化回归模型之间的主要区别之一是后者涉及调整超参数 lambda。上面的代码针对不同的 lambda 值多次运行glmnet()模型。我们可以使用cv.glmnet()函数自动完成查找最佳 lambda 值的任务。这是使用下面的代码行执行的。

      cv_ridge <- cv.glmnet(x, y_train, alpha = 0, lambda = lambdas)
optimal_lambda <- cv_ridge$lambda.mi

_{免责声明：本内容来源于第三方作者授权、网友推荐或互联网整理，旨在为广大用户提供学习与参考之用。所有文本和图片版权归原创网站或作者本人所有，其观点并不代表本站立场。如有任何版权侵犯或转载不当之情况，请与我们取得联系，我们将尽快进行相关处理与修改。感谢您的理解与支持！}

_查看原文

技术指南

R语言

阅读全文

使用 R 进行线性、套索和岭回归

杭州电子商务研究院

5年前 · 面向社会、服务行业、政产学研结合、整合资源、和谐发展

Observations: 574 Variables: 5 $ pce <dbl> 507.4, 510.5, 516.3, 512.9, 518.1, 525.8, 531.5, 534.2, 544.9, 544.6, 5... $ pop <dbl> 198.7, 198.9, 199.1, 199.3, 199.5, 199.7, 199.8, 199.9, 200.1, 200.2, 2... $ psavert <dbl> 12.5, 12.5, 11.7, 12.5, 12.5, 12.1, 11.7, 12.2, 11.6, 12.2, 12.0, 11.6,... $ uempmed <dbl> 4.5, 4.7, 4.6, 4.9, 4.7, 4.8, 5.1, 4.5, 4.1, 4.6, 4.4, 4.4, 4.5, 4.2, 4... $ unemploy <dbl> 2.9, 2.9, 3.0, 3.1, 3.1, 3.0, 2.9, 3.0, 2.9, 2.7, 2.7, 2.9, 2.9, 2.8, 2...

cols = c('pce', 'pop', 'psavert', 'uempmed') pre_proc_val <- preProcess(train[,cols], method = c("center", "scale")) train[,cols] = predict(pre_proc_val, train[,cols]) test[,cols] = predict(pre_proc_val, test[,cols]) summary(train)

pce pop psavert uempmed unemploy Min. :-1.2135 Min. :-1.6049 Min. :-1.89481 Min. :-1.1200 Min. : 2.700 1st Qu.:-0.8856 1st Qu.:-0.8606 1st Qu.:-0.76518 1st Qu.:-0.6198 1st Qu.: 6.500 Median :-0.2501 Median :-0.1416 Median :-0.05513 Median :-0.2447 Median : 7.500 Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.0000 Mean : 7.791 3rd Qu.: 0.7797 3rd Qu.: 0.9060 3rd Qu.: 0.81629 3rd Qu.: 0.1055 3rd Qu.: 8.600 Max. : 2.1244 Max. : 1.8047 Max. : 2.91417 Max. : 4.1571 Max. :15.300

Call: lm(formula = unemploy ~ uempmed + psavert + pop + pce, data = train) Residuals: Min 1Q Median 3Q Max -2.4262 -0.7253 0.0278 0.6697 3.2753 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.79077 0.04712 165.352 < 2e-16 *** uempmed 2.18021 0.08588 25.386 < 2e-16 *** psavert 0.79126 0.13244 5.975 5.14e-09 *** pop 5.95419 0.37405 15.918 < 2e-16 *** pce -5.31578 0.32753 -16.230 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9435 on 396 degrees of freedom Multiple R-squared: 0.8542, Adjusted R-squared: 0.8527 F-statistic: 579.9 on 4 and 396 DF, p-value: < 2.2e-16

#Step 1 - create the evaluation metrics function eval_metrics = function(model, df, predictions, target){ resids = df[,target] - predictions resids2 = resids**2 N = length(predictions) r2 = as.character(round(summary(model)$r.squared, 2)) adj_r2 = as.character(round(summary(model)$adj.r.squared, 2)) print(adj_r2) #Adjusted R-squared print(as.character(round(sqrt(sum(resids2)/N), 2))) #RMSE } # Step 2 - predicting and evaluating the model on train data predictions = predict(lr, newdata = train) eval_metrics(lr, train, predictions, target = 'unemploy') # Step 3 - predicting and evaluating the model on test data predictions = predict(lr, newdata = test) eval_metrics(lr, test, predictions, target = 'unemploy')

cols_reg = c('pce', 'pop', 'psavert', 'uempmed', 'unemploy') dummies <- dummyVars(unemploy ~ ., data = dat[,cols_reg]) train_dummies = predict(dummies, newdata = train[,cols_reg]) test_dummies = predict(dummies, newdata = test[,cols_reg]) print(dim(train_dummies)); print(dim(test_dummies))

library(glmnet) x = as.matrix(train_dummies) y_train = train$unemploy x_test = as.matrix(test_dummies) y_test = test$unemploy lambdas <- 10^seq(2, -3, by = -.1) ridge_reg = glmnet(x, y_train, nlambda = 25, alpha = 0, family = 'gaussian', lambda = lambdas) summary(ridge_reg)

Length Class Mode a0 51 -none- numeric beta 204 dgCMatrix S4 df 51 -none- numeric dim 2 -none- numeric lambda 51 -none- numeric dev.ratio 51 -none- numeric nulldev 1 -none- numeric npasses 1 -none- numeric jerr 1 -none- numeric offset 1 -none- logical call 7 -none- call nobs 1 -none- numeric