Erros padrão para previsão de laço usando R

60

Estou tentando usar um modelo LASSO para previsão e preciso estimar erros padrão. Certamente alguém já escreveu um pacote para fazer isso. Mas, até onde posso ver, nenhum dos pacotes no CRAN que fazem previsões usando um LASSO retornará erros padrão para essas previsões.

Portanto, minha pergunta é: Existe um pacote ou algum código R disponível para calcular erros padrão para previsões do LASSO?

Rob Hyndman
fonte
3
Para esclarecer a natureza subjacente dessa pergunta (uma vez que ela está indo e voltando b / t CV & SO), pergunto-me se poderíamos editar o título, Rob. Que tal 'Por que não parece haver um pacote para erros padrão do LASSO, eles são difíceis de calcular?', Ou algo assim, talvez acoplado a algumas edições menores no corpo para torná-lo consistente. Eu acho que isso tornaria mais claramente o tópico no CV, para que essa ambiguidade não surja e não tenhamos que ir e voltar.
gung - Reintegrar Monica
3
Eu poderia fazer a pergunta mais sobre a metodologia estatística, mas não era exatamente isso que eu queria saber. Deve haver um lugar para perguntas no CV sobre qual software implementa um determinado método. Discussão adicional em meta.stats.stackexchange.com/q/2007/159
Rob Hyndman
11
Você pode fazer isso facilmente em uma estrutura bayesiana usando o pacote monomvn, veja minha resposta abaixo.
Fabian9 /

Respostas:

20

Em uma nota relacionada, que pode ser útil, Tibshirani e colegas propuseram um teste de significância para o laço. O artigo está disponível e intitulado "Um teste de significância para o laço". Uma versão gratuita do artigo pode ser encontrada aqui

julieth
fonte
Link sem paywall para o artigo que você menciona: statweb.stanford.edu/~tibs/ftp/covtest.pdf
mvherweg
13

A resposta Sandipan Karmakar diz o que você deve fazer, isso deve ajudá-lo no "como":

> library(monomvn)
>
> ## following the lars diabetes example
> data(diabetes)
> str(diabetes)
'data.frame':   442 obs. of  3 variables:
 $ x : AsIs [1:442, 1:10] 0.038075.... -0.00188.... 0.085298.... -0.08906.... 0.005383.... ...
      ..- attr(*, "dimnames")=List of 2
      .. ..$ : NULL
  .. ..$ : chr  "age" "sex" "bmi" "map" ...

 $ y : num  151 75 141 206 135 97 138 63 110 310 ...

[...]

> ## Bayesian Lasso regression
> reg_blas <- with(diabetes, blasso(x, y))
t=100, m=8
t=200, m=5
t=300, m=8
t=400, m=8
t=500, m=7
t=600, m=8
t=700, m=8
t=800, m=8
t=900, m=5
> 
> ## posterior mean beta (setting those with >50% mass at zero to exactly zero)
> (beta <- colMeans(reg_blas$beta) * (colMeans(reg_blas$beta != 0)  > 0.5))
      b.1       b.2       b.3       b.4       b.5       b.6       b.7       b.8 
   0.0000 -195.9795  532.7136  309.1673 -101.1288    0.0000 -196.4315    0.0000 
      b.9      b.10 
 505.4726    0.0000 
> 
> ## n x nsims matrix of realizations from the posterior predictive:
> post_pred_y <- with(reg_blas, X %*% t(beta))
> 
> ## predictions:
> y_pred <- rowMeans(post_pred_y)
> head(y_pred)
[1]  52.772443 -78.690610  24.234753   9.717777 -23.360369 -45.477199
> 
> ## sd of y:
> sd_y <- apply(post_pred_y, 1, sd)
> head(sd_y)
[1] 6.331673 6.756569 6.031290 5.236101 5.657265 6.150473
> 
> ## 90% credible intervals
> ci_y <- t(apply(post_pred_y, 1, quantile, probs=c(0.05, 0.95)))
> head(ci_y)
             5%       95%
[1,]  42.842535  62.56743
[2,] -88.877760 -68.47159
[3,]  14.933617  33.85679
[4,]   1.297094  18.01523
[5,] -32.709132 -14.13260
[6,] -55.533807 -35.77809
fabianos
fonte
13

O LASSO bayesiano é a única alternativa ao problema de cálculo de erros padrão. Erros padrão são calculados automaticamente no LASSO Bayesiano ... Você pode implementar o LASSO Bayesiano com muita facilidade usando o esquema de amostragem Gibbs ...

O LASSO Bayesiano precisa que distribuições anteriores sejam atribuídas aos parâmetros do modelo. No modelo LASSO, temos a função objetivo with como o parâmetro de regularização. Aqui, como temos -norm for , é necessário um tipo especial de distribuição anterior para isso, a distribuição LAPLACE é uma mistura escalável de distribuição normal com distribuição exponencial como densidade de mistura. Com base nas posteriores condicionais completas de cada um dos parâmetros devem ser deduzidos.||yXβ||22+λ||β||1λ1β

Então, pode-se usar Gibbs Sampling para simular a cadeia. Ver Park & ​​Cassella (2008), "The Bayesian Lasso", JASA , 103 , 482 .

Existem três desvantagens inerentes ao LASSO frequentista:

  1. É preciso escolher por validação cruzada ou por outros meios.λ

  2. Os erros padrão são difíceis de calcular, pois o LARS e outros algoritmos produzem estimativas de pontos para .β

  3. A estrutura hierárquica do problema em questão não pode ser codificada usando o modelo frequentista, o que é bastante fácil na estrutura bayesiana.

Sandipan Karmakar
fonte
11

Para adicionar as respostas acima, o problema parece ser que mesmo um bootstrap é provavelmente insuficiente, pois a estimativa do modelo penalizado é tendenciosa e o bootstrapping fala apenas da variação - ignorando o viés da estimativa. Isso está bem resumido na vinheta do pacote penalizado na página 18 .

No entanto, se estiver sendo usado para previsão, por que é necessário um erro padrão do modelo? Não é possível cruzar a validação ou a inicialização adequada e produzir um erro padrão em torno de uma métrica relacionada à previsão, como o MSE?

B_Miner
fonte
3
O bootstrapping pode estimar e corrigir o viés, embora as amostras precisem ser bastante grandes.
Glen_b
3

Existe o pacote selectionInference em R, https://cran.r-project.org/web/packages/selectiveInference/index.html , que fornece intervalos de confiança ep valores para seus coeficientes ajustados pelo LASSO, com base no artigo a seguir :

Stephen Reid, Jerome Friedman e Rob Tibshirani (2014). Um estudo de estimativa de variância de erro na regressão do laço. arXiv: 1311.5274

PS: apenas perceba que isso produz estimativas de erro para seus parâmetros, não tenho certeza sobre o erro em sua previsão final, se é isso que você procura ... Suponho que você possa usar "intervalos de previsão de população" para isso, se quiser (por parâmetros de reamostragem de acordo com o ajuste após uma distribuição normal multivariada).

Tom Wenseleers
fonte