LASSO回归，全称：least absolute shrinkage and selection operator（最小绝对值收敛和选择算子算法），是一种压缩估计。在拟合广义线性模型的同时进行
（variable selection）和
（regularization），以提高所得统计模型的预测准确性和可解释性
LASSO回归复杂度调整的程度由参数 λ 来控制，λ 越大对变量较多的线性模型的惩罚力度就越大，从而最终获得一个变量较少的模型。 LASSO回归与
同属于一个被称为
的广义线性模型家族。 这一家族的模型除了相同作用的参数 λ之外，还有另一个参数 α来控制应对高相关性(highly correlated)数据时模型的性状。 LASSO回归 α=1，Ridge回归 α=0，一般Elastic Net模型 0<α<1
：加载数据
# install.packages("glmnet")
library(glmnet)
## 定义响应变量
y <- mtcars$hp
## 定义预测变量矩阵
x <- data.matrix(mtcars[, c('mpg', 'wt', 'drat', 'qsec')])
：拟合 Lasso 回归模型
library(glmnet)
# 执行k折交叉验证(https://www.statology.org/k-fold-cross-validation/)并确定产生最低测试均方误差 (MSE) 的 lambda 值
cv_model <- cv.glmnet(x, y, alpha = 1)
# 找寻最小化测试 MSE 的最佳 lambda 值
best_lambda <- cv_model$lambda.min #l
best_lambda
## [1] 2.215202
# produce plot of test MSE by lambda value
plot(cv_model)
glmnet-1
★ 图中有两条虚线，左边为
★ ambda.min是指在所有的λ值中，得到最小目标参量均值的那一个；
lambda.1se是指在 lambda.min一个方差范围内得到最简单模型的那一个λ值；
λ值到达一定大小之后，继续增加模型自变量个数即缩小λ值，并不能很显著的提高模型性能， lambda.1se 给出的就是一个具备优良性能但是自变量个数最少的模型
：分析最终模型
# 找寻最佳模型系数
best_model <- glmnet(x, y, alpha = 1 , lambda = best_lambda)
coef(best_model)
## 5 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 486.013958
## mpg -2.916499
## wt 21.989615
## drat .
## qsec -19.692037
print(best_model)
## Call: glmnet(x = x, y = y, alpha = 1, lambda = best_lambda)
## Df %Dev Lambda
## 1 3 80.47 2.215
glmnet()参数 family 规定了回归模型的类型：
family = "gaussian" 适用于一维连续因变量（univariate）
family = "mgaussian" 适用于多维连续因变量（multivariate）
family = "poisson" 适用于非负次数因变量（count）
family = "binomial" 适用于二元离散因变量（binary）
family = "multinomial" 适用于多元离散因变量（category）
★ 没有显示预测变量 drat的系数，因为套索回归将系数一直缩小到零。这意味着它完全从模型中，因为它的影响力不够。
★ 和
之间的关键区别。Ridge回归将所有系数缩小为零，但LASSO回归有可能通过将系数完全缩小为零来从模型中删除预测变量。
# 定义一个新矩阵
new = matrix(c(24, 2.5, 3.5, 18.5), nrow = 1, ncol = 4)
# 使用lasso回归模型预测响应值
predict(best_model, s = best_lambda, newx = new)
## s1
[1,] 106.6893
★ 根据输入值，模型预测这辆车的hp值为106.6893
# 使用拟合最佳模型进行预测
y_predicted <- predict(best_model, s = best_lambda, newx = x)
y_predicted
## s1
## Mazda RX4 158.24933
## Mazda RX4 Wag 152.82914
## Datsun 710 104.06486
## Hornet 4 Drive 111.48428
## Hornet Sportabout 171.96122
## Valiant 111.13639
## Duster 360 210.88907
## Merc 240D 91.15750
## Merc 230 37.83740
## Merc 280 145.29716
# find SST（离差平方和 ） and SSE（残差平方和）
sst <- sum((y - mean(y))^2)
sse <- sum((y_predicted - y)^2)
# find R-Squared
rsq <- 1 - sse/sst
rsq
## [1] 0.8047064
★ 最佳模型能够训练数据响应值中80.47%的变异
：拟合广义线性模型
从上面的模型可以看出，drat的系数因为影响力不够被删除，即，可把剩下的参数mpg、wt和qsec，拿出来组成广义线性方程
mod <- glm(y~mpg+wt+qsec, data = data.frame(x))
summary(mod)
## Call:
## glm(formula = y ~ mpg + wt + qsec, data = data.frame(x))
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -48.392 -14.285 -6.816 12.320 97.624
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 494.570 78.276 6.318 7.81e-07 ***
## mpg -2.700 2.269 -1.190 0.2442
## wt 25.042 12.892 1.942 0.0622 .
## qsec -20.966 3.865 -5.425 8.68e-06 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## (Dispersion parameter for gaussian family taken to be 1006.536)
## Null deviance: 145727 on 31 degrees of freedom
## Residual deviance: 28183 on 28 degrees of freedom
## AIC: 317.8
## Number of Fisher Scoring iterations: 2
★ 依据上面回归系数的结果，有2个指标（wt和qsec）入选，写出logistic回归的方程式：
Logit(y)=494.570 + 25.042 × wt - 20.966 × qsec
：Hosmer-Lemeshow拟合优度
# install.packages("ResourceSelection")
library(ResourceSelection)
hoslem.test(mod$y,fitted(mod))
## Hosmer and Lemeshow goodness of fit (GOF) test
## data: mod$y, fitted(mod)
## X-squared = -1.4297, df = 8, p-value = 1
★ p = 1，表明无法拒绝原假设
网友评论