• R语言统计学

  • 请问R语言如何用lars包中的LASSO方法进行变量选择？

请问R语言如何用lars包中的LASSO方法进行变量选择？我在做的过程中出现Cp值为空，R-square为1，但是参数估计值并不接近真实值，正确的做法应该怎么做？

library(MASS) #generate random number
library(lars)#LASSO
cov_mat<-function(p)
{
covm<-matrix(0,p,p)
for (i in 1:p)
{
for (j in 1:p)
{
covm[i,j]<-0.5^(abs(i-j))
}
}
return(covm)
}
p<-10
n<-200
set.seed(666)
epsilon<-rnorm(n,0,1)
beta<-c(1,1.5,2,1,0,0,0,0,0,0)
set.seed(666)
x<-mvrnorm(n,rep(0,p),cov_mat(p))
y<-x%*%beta+epsilon
LASSOres<-lars(as.matrix(x),as.matrix(y),type="lasso")
结果如下
LASSOres

Call:
lars(x = as.matrix(x), y = as.matrix(y), type = "lasso")
R-squared: 1
Sequence of LASSO moves:


Var 3 2 1 4 7 9 8 6 5 10
Step 1 2 3 4 5 6 7 8 9 10

LASSOres$beta[10,]
[1] 0.81909408 1.28566343 1.71271692 0.52252613 -0.02578134 -0.19780722 -0.19525537 -0.12094906
[9] -0.16279040 0.00000000

6 天 后

先谈一个代码中的问题：你在生成epsilon和x时候两次都set.seed(666)。而mvrnorm本质上是一个rnorm然后再用协方差矩阵变换到你需要的协方差结构上，这就导致你的x本质上是epsilon“完全一致”（在一些线性变换之后）的，这导致了你的这个模拟的数据生成是有问题的。如果为了数据的可重复性的话，只需要最初的一次set.seed即可。

建议先去阅读一下lars的文档，里面其实写了lars是把整个solution path都计算了，从最开始一个变量都没选，到最后所有变量包含进来按照简单线性回归整个过程都给了。你如果plot(LASSOres)的话可以看到求解路径的图，另外这个打印结果其实也就是每一步到底选择进入模型哪一个变量。在你的函数中可以增加trace = TRUE参数要求lars打印每一次选择的结果。当然作为“上帝”的话你还可以加intercept = FALSE参数要求模型不需要拟合截距项。

至于如何在整个求解路径上选择合适的模型结果，也就是选择惩罚函数中的lambda的值，这个就是个相对玄学一点的问题。可以考虑使用AIC、BIC一类的准则，也可以使用CV这种数据导向的方法（CV之后也有cv.min或者cv.1se之类的选择区别），甚至也可以人为预先给定说我就是要前 x 个最先进入模型的变量。

关于估计精度问题，lasso本身就是一个有偏的估计量，lasso估计结果的估计相合性和变量选择相合性是不能同时达到的。也就是说选的准的Lasso估计值有偏，估计值（针对真实变量）准确的lasso选不准（一般是多选无关变量）。考虑一些非凸惩罚，如SCAD、MCP之类可以同时做到变量选择和估计结果的相合性。

当然以上的说法都建立在理论性质的讨论上，有各种对数据的假设条件，也有对样本量趋于无穷，以及其他一些参数相应的阶数要求。实际数据分析中并不一定能保证理论效果的完全一致。工具方面，对于lasso可以考虑glmnet包，在cv.glmnet中有cvmin和cv1se两种选择。对于MCP和SCAD可以考虑ncvreg包，也已经内含了cv.ncvreg。

楼上正解。这里用我的 msaenet 给出一个实现，大致对应 lasso 和 adaptive SCAD:

library(msaenet)

dat <- msaenet.sim.gaussian(
  n = 400,
  p = 10,
  rho = 0.5,
  coef = c(1, 1.5, 2, 1),
  snr = 2,
  p.train = 0.5,
  seed = 42
)

fit_lasso <- aenet(
  x = dat$x.tr,
  y = dat$y.tr,
  init = "enet",
  alphas = 1,
  tune = "cv",
  nfolds = 5L,
  rule = "lambda.1se",
  seed = 42
)

fit_lasso$beta.first

#> 10 x 1 sparse Matrix of class "dgCMatrix"
#> s0
#> V1  0.6096754
#> V2  0.9004872
#> V3  2.4194953
#> V4  0.5458725
#> V5  .
#> V6  .
#> V7  .
#> V8  .
#> V9  .
#> V10 .

fit_ascad <- suppressWarnings(asnet(
  x = dat$x.tr,
  y = dat$y.tr,
  init = "ridge",
  alphas = 1,
  tune = "cv",
  nfolds = 5,
  seed = 42
))

fit_ascad$beta

#> 10 x 1 sparse Matrix of class "dgCMatrix"
#>
#>  [1,] 1.0179300
#>  [2,] 1.1580254
#>  [3,] 2.5748324
#>  [4,] 0.9569498
#>  [5,] .
#>  [6,] .
#>  [7,] .
#>  [8,] .
#>  [9,] .
#> [10,] .

可以看到另外很重要的一点是推荐使用 signal-to-noise ratio (SNR) 来驱动模拟数据的生成.