Usando o LASSO do pacote lars (ou glmnet) no R para seleção de variáveis

Desculpe se esta pergunta é um pouco básica.

Eu estou procurando usar a seleção de variáveis ​​do LASSO para um modelo de regressão linear múltipla em R. Eu tenho 15 preditores, um dos quais é categórico (isso causará um problema?). Após definir e , utilizo os seguintes comandos:$x$$y$

model = lars(x, y)
coef(model)

Meu problema é quando eu uso coef(model). Isso retorna uma matriz com 15 linhas, com um preditor extra adicionado a cada vez. No entanto, não há sugestão de qual modelo escolher. Perdi alguma coisa? Existe uma maneira de obter o pacote lars para retornar apenas um modelo " melhor "?

Há outras postagens sugerindo o uso, glmnetmas isso parece mais complicado. Uma tentativa é a seguinte, usando o mesmo e . Perdi alguma coisa aqui ?: $x$$y$

cv = cv.glmnet(x, y)
model = glmnet(x, y, type.gaussian="covariance", lambda=cv$lambda.min)
predict(model, type="coefficients")

O comando final retorna uma lista de minhas variáveis, a maioria com um coeficiente, embora algumas sejam = 0. Essa é a escolha correta do " melhor " modelo selecionado pelo LASSO? Se eu encaixar um modelo linear com todas as minhas variáveis ​​com coeficientes not=0, fico com estimativas de coeficientes muito semelhantes, mas um pouco diferentes. Existe uma razão para essa diferença? Seria aceitável refazer o modelo linear com essas variáveis ​​escolhidas pelo LASSO e tomá-lo como meu modelo final? Caso contrário, não vejo valores de p para significância. Perdi alguma coisa?

Faz

type.gaussian="covariance" 

garantir que isso glmnetusa regressão linear múltipla?

A normalização automática das variáveis ​​afeta os coeficientes? Existe alguma maneira de incluir termos de interação em um procedimento LASSO?

Estou procurando usar esse procedimento mais como uma demonstração de como o LASSO pode ser usado do que para qualquer modelo que realmente seja usado para qualquer inferência / previsão importante, se isso mudar alguma coisa.

Obrigado por separar um tempo para ler isso. Quaisquer comentários gerais sobre o LASSO / lars / glmnet também serão muito apreciados.

O uso glmneté realmente fácil quando você o compreende, graças à excelente vinheta em http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html (você também pode verificar a página do pacote CRAN). Quanto ao melhor lambda para glmnet, a regra geral é usar

cvfit <- glmnet::cv.glmnet(x, y)
coef(cvfit, s = "lambda.1se")

em vez de lambda.min.

Para fazer o mesmo, larsvocê deve fazê-lo manualmente. Aqui está a minha solução

cv <- lars::cv.lars(x, y, plot.it = FALSE, mode = "step")
idx <- which.max(cv$cv - cv$cv.error <= min(cv$cv))
coef(lars::lars(x, y))[idx,]

Lembre-se de que isso não é exatamente o mesmo, porque isso é interrompido em um nó de laço (quando uma variável entra) em vez de a qualquer momento.

Observe que agora glmneté o pacote preferido, ele é mantido ativamente, mais do que isso lars, e que houve perguntas sobre glmnetvs larsrespondidas anteriormente (os algoritmos usados ​​diferem).

Quanto à sua questão de usar o laço para escolher variáveis ​​e depois ajustar o OLS, é um debate em andamento. O Google para OLS postou Lasso e existem alguns artigos discutindo o tópico. Até os autores de Elements of Statistical Learning admitem que é possível.

Edit : Aqui está o código para reproduzir com mais precisão o que glmnetfazlars

  cv <- lars::cv.lars(x, y, plot.it = FALSE)
  ideal_l1_ratio <- cv$index[which.max(cv$cv - cv$cv.error <= min(cv$cv))]
  obj <- lars::lars(x, y)
  scaled_coefs <- scale(obj$beta, FALSE, 1 / obj$normx)
  l1 <- apply(X = scaled_coefs, MARGIN = 1, FUN = function(x) sum(abs(x)))
  coef(obj)[which.max(l1 / tail(l1, 1) > ideal_l1_ratio),]

Estou voltando a esta pergunta há algum tempo, pois acho que resolvi a solução correta.

Aqui está uma réplica usando o conjunto de dados mtcars:

library(glmnet)
`%ni%`<-Negate(`%in%')
data(mtcars)

x<-model.matrix(mpg~.,data=mtcars)
x=x[,-1]

glmnet1<-cv.glmnet(x=x,y=mtcars$mpg,type.measure='mse',nfolds=5,alpha=.5)

c<-coef(glmnet1,s='lambda.min',exact=TRUE)
inds<-which(c!=0)
variables<-row.names(c)[inds]
variables<-variables[variables %ni% '(Intercept)']

'variable' fornece a lista das variáveis ​​que resolvem a melhor solução.

Talvez a comparação com a regressão gradual da seleção direta ajude (consulte o seguinte link para um site de um dos autores http://www-stat.stanford.edu/~tibs/lasso/simple.html) Essa é a abordagem usada no capítulo 3.4.4 de Os elementos do aprendizado estatístico (disponível on-line gratuitamente). Pensei que o capítulo 3.6 desse livro ajudasse a entender a relação entre mínimos quadrados, melhor subconjunto e laço (além de alguns outros procedimentos). Também acho útil fazer a transposição do coeficiente t (coef (modelo)) e write.csv, para que eu possa abri-lo no Excel junto com uma cópia do gráfico (modelo) ao lado. Você pode classificar pela última coluna, que contém a estimativa de mínimos quadrados. Então você pode ver claramente como cada variável é adicionada em cada etapa por partes e como os coeficientes mudam como resultado. Claro que essa não é a história toda, mas espero que seja um começo.

larse glmnetoperar em matrizes brutas. Para incluir termos de interação, você precisará construir as matrizes. Isso significa uma coluna por interação (que é por nível por fator, se você tiver fatores). Veja lm()como funciona (aviso: existem dragões).

Para fazer isso agora, faça algo como: Para criar um termo de interação manualmente, você pode (mas talvez não deva , porque é lento):

int = D["x1"]*D["x2"]
names(int) = c("x1*x2")
D = cbind(D, int)

Em seguida, para usar isso no lars (supondo que você tenha uma ychance):

lars(as.matrix(D), as.matrix(y))

Eu gostaria de poder ajudá-lo mais com as outras perguntas. Encontrei este porque o lars está me dando dor e a documentação nele e na web é muito pequena.

LARS resolve o caminho da solução INTEIRA. O caminho da solução é linear por partes - há um número finito de pontos de "entalhe" (ou seja, valores do parâmetro de regularização) nos quais a solução é alterada.

Portanto, a matriz de soluções que você está obtendo é todas as soluções possíveis. Na lista que ele retorna, também deve fornecer os valores do parâmetro de regularização.