R實戰(zhàn) | Lasso回歸模型建立及變量篩選
Tibshirani(1996) 引入了 LASSO (Least Absolute Shrinkage and Selection Operator)模型峭沦，用于參數(shù)的選擇和收縮曹货。當(dāng)我們分析大數(shù)據(jù)時锅尘，這個模型非常有用烟具。在這篇文章中，我們學(xué)習(xí)如何使用R包glmnet 包建立LASSO 模型庐杨。
原文：https://mp.weixin.qq.com/s/Q-76vY6hr81Sh9Zs5KYcCg
這些回歸模型被稱為正則化或懲罰回歸模型咽笼。Lasso可以用于變量數(shù)量較多的大數(shù)據(jù)集。傳統(tǒng)的線性回歸模型無法處理這類大數(shù)據(jù)暖释。
雖然線性回歸估計器 (linear regression estimator)在偏-方差權(quán)衡關(guān)系方面是無偏估計器袭厂，但正則化或懲罰回歸，如Lasso, Ridge承認一些減少方差的偏倚球匕。這意味著后者的最小化問題有兩個組成部分:均方誤差(linear regression estimator)和懲罰參數(shù)()纹磺。Lasso的L1懲罰使變量選擇和收縮成為可能，而Ridge的L2懲罰使變量收縮成為可能亮曹。
[圖片上傳失敗...(image-af3b85-1651640667147)]
關(guān)于正則化詳見：零基礎(chǔ)"機器學(xué)習(xí)"自學(xué)筆記|Note8:正則化
X變量應(yīng)該用均值零和單位方差進行標準化橄杨，因為變量的尺度差異往往會使懲罰分配不均。
由上式可知照卦，第一部分是殘差平方和(Residual Sum of Squares,RSS)式矫，第二部分是懲罰項。該罰項由超參數(shù) λ 調(diào)整役耕。超參數(shù)由用戶通過人工搜索或交叉驗證的方式外源性給出采转。
當(dāng)Lasso中包含了某些變量，但RSS值的降低很小瞬痘，可以忽略不計時故慈，收縮懲罰的影響就會增加板熊。這意味著這個變量的系數(shù)是零(Lasso)或接近零(Ridge)。
本教程提供了一個循序漸進的例子察绷，如何在R中執(zhí)行套索回歸干签。
[TOC]
22
LASSO
示例數(shù)據(jù)準備及預(yù)處理
# install.packages('glmnet')
library(glmnet)
graphics.off() # clear all graphs
rm(list = ls())
# 示例數(shù)據(jù)準備
N = 500 # 觀測數(shù)
p = 20 # 變量數(shù)
# X variable
X = matrix(rnorm(N*p), ncol=p)
# 計算標準化前的均值和標準差
colMeans(X) # mean
apply(X,2,sd) # standard deviation
# 標準化
X = scale(X,center = T,scale = T)
# 計算標準化后的均值和標準差
colMeans(X) # mean
apply(X,2,sd) # standard deviation
#——————————————-
# Y variable
#——————————————-
beta = c( 0.15, -0.33, 0.25, -0.25, 0.05,rep(0, p/2-5),
-0.25, 0.12, -0.125, rep(0, p/2-3))
# Y variable, standardized Y
y = X%*%beta + rnorm(N, sd=0.5)
y = scale(y)
注：對實際數(shù)據(jù)，只需分別對
X拆撼，y進行標準化即可筒严，即scale()。
LASSO 模型建立
# Model
# 當(dāng)lambda = 0.01
lambda <- 0.01
# lasso
la.eq <- glmnet(X, y, lambda=lambda,
family='gaussian',
intercept = F, alpha=1)
# 當(dāng)alpha設(shè)置為0則為ridge回歸情萤，將alpha設(shè)置為0和1之間則為elastic net
# 系數(shù)結(jié)果 (lambda=0.01)
la.eq$beta[,1]
# Lasso篩選變量動態(tài)過程圖
la.eq <- glmnet(X, y, family="gaussian",
intercept = F, alpha=1)
# plot
plot(la.eq,xvar = "lambda", label = F)
# 也可以用下面的方法繪制
#matplot(log(la.eq$lambda), t(la.eq$beta),
# type="l", main="Lasso", lwd=2)
我們可以看到鸭蛙，當(dāng)lambda越大，各估計參數(shù)相應(yīng)的也被壓縮得更小筋岛，而當(dāng)lambda達到一定值以后娶视，一部分不重要的變量將被壓縮為0，代表該變量已被剔除出模型睁宰，圖中從左至左右斷下降的曲線如同被不斷增大的lambda一步一步壓縮肪获，直到壓縮為0。
變量篩選
模型已經(jīng)跑出來了柒傻，如何篩選變量呢孝赫？我們還得確定lambda，也就是在上圖中畫一條垂直于橫軸的直線红符，這樣我們才能知道哪些變量被壓縮為0青柄，以及未被壓縮為0的變量的系數(shù)的估計值究竟是多少。
那我們?nèi)绾芜x擇lambda呢预侯？我們可以使用R的glmnet中另一個函數(shù)cv.glmnet致开。這個函數(shù)使用的是“交叉驗證”挑選lambda。
# Run cross-validation & select lambda
#————————————————
mod_cv <- cv.glmnet(x=X, y=y, family="gaussian", # 默認nfolds = 10
intercept = F, alpha=1)
plot(mod_cv)
通過交叉驗證萎馅，我們可以選擇平均誤差最小的那個λ双戳，即mod_cv$lambda.min也可以選擇平均誤差在一個標準差以內(nèi)的最大的λ，即mod_cv$lambda.1se糜芳。
# lambda.min : the λ at which the minimal MSE is achieved.
# lambda.1se : the largest λ at which the MSE is within one standard error of the minimal MSE.
print(paste(mod_cv$lambda.min,
log(mod_cv$lambda.min)))
print(paste(mod_cv$lambda.1se,
log(mod_cv$lambda.1se)))
# 這里我們以lambda.min為最優(yōu) λ
best_lambda <- mod_cv$lambda.min
best_lambda
最后飒货，我們可以分析由最優(yōu)lambda值產(chǎn)生的最終模型。
# 最終模型的系數(shù)估計
#find coefficients of best model
best_model <- glmnet(X, y, alpha = 1, lambda = best_lambda)
coef(best_model)
> coef(best_model)
21 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) -9.499626e-18
V1 2.156750e-01
V2 -3.636329e-01
V3 3.247269e-01
V4 -2.554500e-01
V5 3.571598e-02
V6 -7.033558e-03
V7 .
V8 1.405279e-04
V9 .
V10 -1.285986e-03
V11 -3.187390e-01
V12 9.353756e-02
V13 -1.360398e-01
V14 .
V15 7.314506e-02
V16 .
V17 2.767650e-02
V18 3.150004e-03
V19 .
V20 -3.646384e-03
如變量沒有顯示系數(shù)峭竣，即lasso回歸收縮系數(shù)為零塘辅。這意味著它完全被排除在模型之外，因為它的影響力不夠邪驮。系數(shù)非0的變量即為我們篩選的重要特征莫辨。
使用最終模型進行預(yù)測
我們還可以使用最終的lasso回歸模型對新的觀測進行預(yù)測。
# 新觀測
new = matrix(rnorm(20), nrow=1, ncol=20)
#使用 lasso 回歸模型預(yù)測
predict(best_model, s = best_lambda, newx = new)
s1
[1,] 0.7519802
最后，我們可以在訓(xùn)練數(shù)據(jù)上計算模型的R-squared:
y_predicted <- predict(best_model, s = best_lambda, newx = X)
#find SST and SSE
sst <- sum((y - mean(y))^2)
sse <- sum((y_predicted - y)^2)
#find R-Squared
rsq <- 1 - sse/sst
rsq
> rsq
[1] 0.5444342
R2等于0.8047064沮榜。也就是說盘榨，最佳模型能夠解釋訓(xùn)練數(shù)據(jù)響應(yīng)值變化的80.47%。
示例數(shù)據(jù)和代碼領(lǐng)取
詳見：https://mp.weixin.qq.com/s/Q-76vY6hr81Sh9Zs5KYcCg
參考
- Lasso Regression in R (Step-by-Step) (statology.org)(https://www.statology.org/lasso-regression-in-r/)
- Lasso Regression Model with R code | R-bloggers(https://www.r-bloggers.com/2021/05/lasso-regression-model-with-r-code/)
往期內(nèi)容
- (免費教程+代碼領(lǐng)取)|跟著Cell學(xué)作圖系列合集
- Q&A | 如何在論文中畫出漂亮的插圖蟆融？
- Front Immunol 復(fù)現(xiàn) | 1. GEO數(shù)據(jù)下載及sva批次校正(PCA可視化)
- R繪圖 | 氣泡散點圖+擬合曲線
[圖片上傳失敗...(image-e5ad83-1651640667147)]
