关于Lasso回归的一个例子

转载自：http://site.douban.com/182577/widget/notes/10567181/note/285262818/

#给一位朋友做的例子。

 Lasso,套索。一种变量选择方法，使用罚约束来筛掉拟合模型中的系数。 可参考统计学习巨著ESL第2版（ESL这本书的主线可以说就是线性模型加罚约束）。

 这个例子来自吴喜之老师《复杂数据统计方法》p29。第一种方案采用lars包（吴老师书里的方法，细节略有修正）。这个包的方法包括：Least Angle Regression, Lasso 和 Forward Stagewise，作者是两位大牛:Trevor Hastie和 Brad Efron。

library(lars)

data(diabetes)
attach(diabetes)


w <- cbind(diabetes$x, diabetes$y, diabetes$x2)

x <- as.matrix(w[, 1:10])
y <- as.matrix(w[, 11])#响应变量
x2 <- as.matrix(w[, 12:75])
laa <- lars(x2, y) #lars函数默认方法为lasso
plot(laa)

 

summary(laa) #共104步
cva <- cv.lars(x2, y, K = 10, plot.it = TRUE) #10折交叉验证,并绘制cv的变化图

 

best <- cva$index[which.min(cva$cv)]
coef <- coef.lars(laa, mode = "fraction", s = best) #使得CV最小时的系数
coef[coef != 0] #通过CV选择的变量
## sex bmi map hdl ltg glu bmi^2 glu^2 
## -92.967 502.356 241.632 -174.935 465.399 10.829 33.919 62.356 
## age:sex age:map age:ltg age:glu bmi:map 
## 97.308 28.427 4.497 13.780 79.712



laa$Cp[which.min(laa$Cp)] #最小的Cp值及其步数
## 15 
## 18.2

coef1 <- coef.lars(laa, mode = "step", s = 15) #使得Cp最小时的系数（s=15）
coef1[coef1 != 0] #通过CV选择的变量
## sex bmi map hdl ltg glu age^2 bmi^2 
## -112.313 501.617 251.792 -187.828 467.779 18.130 7.609 38.748 
## glu^2 age:sex age:map age:ltg age:glu bmi:map 
## 69.810 107.728 30.045 8.613 11.600 85.688


第2种方案采用glmnet包。这个包的方法是Lasso and elastic-net regularized generalized linear models。这个包的作者是erome Friedman, Trevor Hastie, Rob Tibshiran。对，就是ESL的三位作者。

library(glmnet)

gla <- cv.glmnet(x2, y, nfolds = 10) #cv做交叉验证来确定模型，nfolds=10其实是默认值
plot(gla) #绘制cv变化图

左边线对应最佳lamda，右侧线对应一个SE内最佳模型

gla$lambda.1se #最佳lambda值
## [1] 6.401
gla.best <- gla$glmnet.fit #对应的最佳模型
gla.coef <- coef(gla$glmnet.fit, s = gla$lambda.1se) #系数
gla.coef[which(gla.coef != 0)] #选择的变量
## [1] 152.133 500.895 184.880 -108.360 443.883 6.644 25.214
## [8] 38.632 5.620 1.060 46.970

最终输出的模型系数，对应图三中有两条线，右侧线，就是最佳的lambda值一个标准误之内的最简洁（系数最少）的模型。