Como interpretar glmnet?

Estou tentando ajustar um modelo de regressão linear multivariada com aproximadamente 60 variáveis ​​preditivas e 30 observações, por isso estou usando o pacote glmnet para regressão regularizada porque p> n.

Passei por documentação e outras perguntas, mas ainda não consigo interpretar os resultados, aqui está um código de exemplo (com 20 preditores e 10 observações para simplificar):

Eu crio uma matriz x com num linhas = num observações e num cols = num preditores e um vetor y que representa a variável de resposta

> x=matrix(rnorm(10*20),10,20)
> y=rnorm(10)

Eu ajustei um modelo glmnet deixando alfa como padrão (= 1 para penalidade de laço)

> fit1=glmnet(x,y)
> print(fit1)

Entendo que recebo previsões diferentes com valores decrescentes de lambda (ou seja, penalidade)

Call:  glmnet(x = x, y = y) 

        Df    %Dev   Lambda
  [1,]  0 0.00000 0.890700
  [2,]  1 0.06159 0.850200
  [3,]  1 0.11770 0.811500
  [4,]  1 0.16880 0.774600
   .
   .
   .
  [96,] 10 0.99740 0.010730
  [97,] 10 0.99760 0.010240
  [98,] 10 0.99780 0.009775
  [99,] 10 0.99800 0.009331
 [100,] 10 0.99820 0.008907

Agora, prevejo meus valores Beta, escolhendo, por exemplo, o menor valor lambda fornecido por glmnet

> predict(fit1,type="coef", s = 0.008907)

21 x 1 sparse Matrix of class "dgCMatrix"
                  1
(Intercept) -0.08872364
V1           0.23734885
V2          -0.35472137
V3          -0.08088463
V4           .         
V5           .         
V6           .         
V7           0.31127123
V8           .         
V9           .         
V10          .         
V11          0.10636867
V12          .         
V13         -0.20328200
V14         -0.77717745
V15          .         
V16         -0.25924281
V17          .         
V18          .         
V19         -0.57989929
V20         -0.22522859

Se eu escolher lambda com

cv <- cv.glmnet(x,y)
model=glmnet(x,y,lambda=cv$lambda.min)

Todas as variáveis ​​seriam (.).

Dúvidas e perguntas:

1. Não tenho certeza sobre como escolher lambda.
2. Devo usar as variáveis ​​não (.) Para ajustar outro modelo? No meu caso, eu gostaria de manter o máximo de variáveis ​​possível.
3. Como sei o valor p, ou seja, quais variáveis ​​predizem significativamente a resposta?

Peço desculpas pelo meu pobre conhecimento estatístico! E obrigado por qualquer ajuda.

Aqui está um fato não intuitivo - você não deve fornecer ao glmnet um único valor de lambda. A partir da documentação aqui :

Não forneça um valor único para lambda (para previsões após o uso de CV, prever ()). Forneça uma sequência decrescente de valores lambda. O glmnet confia em seus aquecimentos e começa a velocidade, e geralmente é mais rápido ajustar um caminho inteiro do que calcular um único ajuste.

cv.glmnetajudará você a escolher o lambda, como você mencionou em seus exemplos. Os autores do pacote glmnet sugerem, em cv$lambda.1sevez de cv$lambda.min, mas, na prática, tive sucesso com o último.

Depois de executar o cv.glmnet, você não precisa executar novamente o glmnet! Todo lambda na grade ( cv$lambda) já foi executado. Essa técnica é chamada "Warm Start" e você pode ler mais sobre isso aqui . Parafraseando da introdução, a técnica Warm Start reduz o tempo de execução de métodos iterativos usando a solução de um problema de otimização diferente (por exemplo, glmnet com um lambda maior) como o valor inicial para um problema de otimização posterior (por exemplo, glmnet com um lambda menor )

Para extrair a execução desejada cv.glmnet.fit, tente o seguinte:

small.lambda.index <- which(cv$lambda == cv$lambda.min)
small.lambda.betas <- cv$glmnet.fit$beta[, small.lambda.index]

Revisão (28/1/2017)

Não há necessidade de invadir o objeto glmnet como fiz acima; tomar @ conselho de alex23lemm abaixo e passar o s = "lambda.min", s = "lambda.1se"ou algum outro número (por exemplo, s = .007) para ambos coefe predict. Observe que seus coeficientes e previsões dependem desse valor que é definido por validação cruzada. Use uma semente para a reprodutibilidade! E não esqueça que, se você não fornecer um "s"in coefe predict, estará usando o padrão de s = "lambda.1se". Eu terminei esse padrão depois de vê-lo funcionar melhor em uma pequena situação de dados.s = "lambda.1se"também tende a fornecer mais regularização; portanto, se você estiver trabalhando com alfa> 0, também tenderá a um modelo mais parcimonioso. Você também pode escolher um valor numérico de s com a ajuda de plot.glmnet para chegar a um ponto intermediário (apenas não esqueça de exponenciar os valores do eixo x!).

Q1) Não tenho certeza sobre como escolher lambda. Q2) Devo usar as variáveis ​​não (.) Para ajustar outro modelo? No meu caso, eu gostaria de manter o máximo de variáveis ​​possível.

De acordo com a ótima resposta do @ BenOgorek, normalmente você permite que o ajuste use uma sequência lambda inteira e, ao extrair os coeficientes ideais, use o valor lambda.1se (diferente do que você fez).

Desde que você siga as três advertências abaixo, não lute contra a regularização nem ajuste o modelo: se uma variável foi omitida, foi porque deu uma penalidade geral mais baixa. As advertências são:

1. Para que os coeficientes regularizados sejam significativos, normalize explicitamente a média e o stdev da variável com scale(); não confie glmnet(standardize=T). Para justificação, consulte A padronização antes do Lasso é realmente necessária? ; basicamente uma variável com grandes valores pode ser injustamente punida na regularização.

2. Para ser reproduzível, corra com set.seedvárias sementes aleatórias e verifique a estabilidade dos coeficientes regularizados.

3. Se você deseja uma regularização menos severa, ou seja, inclui mais variáveis, use alfa <1 (ou seja, rede elástica adequada) em vez de crista simples. Eu sugiro que você varra alfa de 0 a 1. Se você fizer isso, para evitar o ajuste excessivo do alfa hiperparâmetro e do erro de regressão, você deve usar a validação cruzada, ou seja, use em cv.glmnet()vez de simples glmnet():

.

for (alpha in c(0,.1,.3,.5,.7,.9,1)) {
  fit <- cv.glmnet(..., alpha=alpha, nfolds=...)
  # Look at the CVE at lambda.1se to find the minimum for this alpha value...
}

Se você deseja automatizar essa pesquisa de grade com o CV, você pode codificá-lo ou usar o pacote caret no glmnet; O sinal de intercalação faz isso bem. Para o cv.glmnet nfoldsvalor do parâmetro, escolha 3 (mínimo) se seu conjunto de dados for pequeno ou 5 ou 10 se for grande.

Q3) Como sei o valor p, ou seja, quais variáveis ​​predizem significativamente a resposta?

Não, eles não são significativos . Conforme explicado em detalhes em Por que é desaconselhável obter informações estatísticas resumidas para os coeficientes de regressão do modelo glmnet?

Apenas deixe cv.glmnet()fazer a seleção da variável automaticamente. Com as advertências acima. E, claro, a distribuição da variável resposta deve ser normal (supondo que você está usando family='gaussian').