关于Lasso回归的一 个例子
关于Lasso回归的一 个例子
#给一位朋友做的例子。
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变化图
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值一个标准误之内的最简洁(系数最少)的模型。
#给一位朋友做的例子。
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值一个标准误之内的最简洁(系数最少)的模型。
什么叫做Cp?
Mallows's Cp, same with AIC BIC
还没看懂。mark先。
求问lasso应该从何学起??本人本科生一枚,学过概统,希望自学lasso,跪求楼主指点
Elements of statistical learning
传说中的ESL
请问最后的glmnet方法选择的最优特征此处是如何解读的
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
> 我来回应