R實(shí)戰(zhàn) | Lasso回歸模型建立及變量篩選
Tibshirani(1996) 引入了 LASSO (Least Absolute Shrinkage and Selection Operator)模型，用于參數(shù)的選擇和收縮。當(dāng)我們分析大數(shù)據(jù)時(shí)，這個(gè)模型非常有用。在這篇文章中，我們學(xué)習(xí)如何使用R包glmnet 包建立LASSO 模型。
原文：https://mp.weixin.qq.com/s/Q-76vY6hr81Sh9Zs5KYcCg
這些回歸模型被稱為正則化或懲罰回歸模型。Lasso可以用于變量數(shù)量較多的大數(shù)據(jù)集。傳統(tǒng)的線性回歸模型無(wú)法處理這類(lèi)大數(shù)據(jù)。
雖然線性回歸估計(jì)器 (linear regression estimator)在偏-方差權(quán)衡關(guān)系方面是無(wú)偏估計(jì)器，但正則化或懲罰回歸，如Lasso, Ridge承認(rèn)一些減少方差的偏倚。這意味著后者的最小化問(wèn)題有兩個(gè)組成部分:均方誤差(linear regression estimator)和懲罰參數(shù)()。Lasso的L1懲罰使變量選擇和收縮成為可能，而Ridge的L2懲罰使變量收縮成為可能。
[圖片上傳失敗...(image-af3b85-1651640667147)]
關(guān)于正則化詳見(jiàn)：零基礎(chǔ)"機(jī)器學(xué)習(xí)"自學(xué)筆記|Note8:正則化
X變量應(yīng)該用均值零和單位方差進(jìn)行標(biāo)準(zhǔn)化，因?yàn)樽兞康某叨炔町愅鶗?huì)使懲罰分配不均。
由上式可知，第一部分是殘差平方和(Residual Sum of Squares,RSS)，第二部分是懲罰項(xiàng)。該罰項(xiàng)由超參數(shù) λ 調(diào)整。超參數(shù)由用戶通過(guò)人工搜索或交叉驗(yàn)證的方式外源性給出。
當(dāng)Lasso中包含了某些變量，但RSS值的降低很小，可以忽略不計(jì)時(shí)，收縮懲罰的影響就會(huì)增加。這意味著這個(gè)變量的系數(shù)是零(Lasso)或接近零(Ridge)。
本教程提供了一個(gè)循序漸進(jìn)的例子，如何在R中執(zhí)行套索回歸。
[TOC]
22
LASSO
示例數(shù)據(jù)準(zhǔn)備及預(yù)處理
# install.packages('glmnet')
library(glmnet)
graphics.off() # clear all graphs
rm(list = ls())
# 示例數(shù)據(jù)準(zhǔn)備
N = 500 # 觀測(cè)數(shù)
p = 20 # 變量數(shù)
# X variable
X = matrix(rnorm(N*p), ncol=p)
# 計(jì)算標(biāo)準(zhǔn)化前的均值和標(biāo)準(zhǔn)差
colMeans(X) # mean
apply(X,2,sd) # standard deviation
# 標(biāo)準(zhǔn)化
X = scale(X,center = T,scale = T)
# 計(jì)算標(biāo)準(zhǔn)化后的均值和標(biāo)準(zhǔn)差
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)
注：對(duì)實(shí)際數(shù)據(jù)，只需分別對(duì)
X，y進(jìn)行標(biāo)準(zhǔn)化即可，即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篩選變量動(dòng)態(tài)過(guò)程圖
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越大，各估計(jì)參數(shù)相應(yīng)的也被壓縮得更小，而當(dāng)lambda達(dá)到一定值以后，一部分不重要的變量將被壓縮為0，代表該變量已被剔除出模型，圖中從左至左右斷下降的曲線如同被不斷增大的lambda一步一步壓縮，直到壓縮為0。
變量篩選
模型已經(jīng)跑出來(lái)了，如何篩選變量呢？我們還得確定lambda，也就是在上圖中畫(huà)一條垂直于橫軸的直線，這樣我們才能知道哪些變量被壓縮為0，以及未被壓縮為0的變量的系數(shù)的估計(jì)值究竟是多少。
那我們?nèi)绾芜x擇lambda呢？我們可以使用R的glmnet中另一個(gè)函數(shù)cv.glmnet。這個(gè)函數(shù)使用的是“交叉驗(yàn)證”挑選lambda。
# Run cross-validation & select lambda
#————————————————
mod_cv <- cv.glmnet(x=X, y=y, family="gaussian", # 默認(rèn)nfolds = 10
intercept = F, alpha=1)
plot(mod_cv)
通過(guò)交叉驗(yàn)證，我們可以選擇平均誤差最小的那個(gè)λ，即mod_cv$lambda.min也可以選擇平均誤差在一個(gè)標(biāo)準(zhǔn)差以內(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ù)估計(jì)
#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
如變量沒(méi)有顯示系數(shù)，即lasso回歸收縮系數(shù)為零。這意味著它完全被排除在模型之外，因?yàn)樗挠绊懥Σ粔颉?strong>系數(shù)非0的變量即為我們篩選的重要特征。
使用最終模型進(jìn)行預(yù)測(cè)
我們還可以使用最終的lasso回歸模型對(duì)新的觀測(cè)進(jìn)行預(yù)測(cè)。
# 新觀測(cè)
new = matrix(rnorm(20), nrow=1, ncol=20)
#使用 lasso 回歸模型預(yù)測(cè)
predict(best_model, s = best_lambda, newx = new)
s1
[1,] 0.7519802
最后，我們可以在訓(xùn)練數(shù)據(jù)上計(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。也就是說(shuō)，最佳模型能夠解釋訓(xùn)練數(shù)據(jù)響應(yīng)值變化的80.47%。
示例數(shù)據(jù)和代碼領(lǐng)取
詳見(jiàn)：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)容
- (免費(fèi)教程+代碼領(lǐng)取)|跟著Cell學(xué)作圖系列合集
- Q&A | 如何在論文中畫(huà)出漂亮的插圖？
- Front Immunol 復(fù)現(xiàn) | 1. GEO數(shù)據(jù)下載及sva批次校正(PCA可視化)
- R繪圖 | 氣泡散點(diǎn)圖+擬合曲線
[圖片上傳失敗...(image-e5ad83-1651640667147)]
