Standardowe błędy do przewidywania lasso za pomocą R.

60

Próbuję użyć modelu LASSO do prognozowania i muszę oszacować standardowe błędy. Z pewnością ktoś już napisał paczkę, aby to zrobić. Ale o ile widzę, żaden z pakietów w CRAN, który wykonuje prognozy za pomocą LASSO, nie zwróci standardowych błędów dla tych prognoz.

Więc moje pytanie brzmi: czy jest dostępny pakiet lub jakiś kod R do obliczenia standardowych błędów dla prognoz LASSO?

Rob Hyndman
źródło
3
Aby wyjaśnić zasadniczy charakter tego pytania (ponieważ jest ono odbijane do tyłu i do przodu b / t CV i SO), zastanawiam się, czy moglibyśmy edytować tytuł, Rob. A może „Dlaczego wydaje się, że nie ma pakietu dla standardowych błędów LASSO, czy są trudne do obliczenia?” Lub coś w tym rodzaju, być może w połączeniu z kilkoma drobnymi zmianami w ciele, aby zapewnić spójność. Sądzę, że dzięki temu byłoby bardziej jednoznacznie na temat CV, aby ta dwuznaczność nie pojawiła się i nie musieliśmy iść do przodu i do tyłu.
gung - Przywróć Monikę
3
Mógłbym zadać pytanie o metodologię statystyczną, ale tak naprawdę nie chciałem wiedzieć. W CV powinno znaleźć się miejsce na pytania dotyczące tego, jakie oprogramowanie implementuje daną metodę. Dalsza dyskusja na stronie meta.stats.stackexchange.com/q/2007/159
Rob Hyndman
1
Możesz to łatwo zrobić w ramach Bayesian przy użyciu pakietu monomvn, zobacz moją odpowiedź poniżej.
fabians

Odpowiedzi:

20

W powiązanej notatce, która może być pomocna, Tibshirani i współpracownicy zaproponowali test istotności dla lasso. Papier jest dostępny, a zatytułowany „Test istotności dla lasso”. Darmowa wersja artykułu znajduje się tutaj

Julieth
źródło
Link bez paywall do artykułu, o którym wspominasz: statweb.stanford.edu/~tibs/ftp/covtest.pdf
mvherweg
13

Odpowiedź Sandipana Karmakara mówi, co robić, to powinno ci pomóc w „jak”:

> 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
fabianie
źródło
13

Bayesian LASSO jest jedyną alternatywą dla problemu obliczania standardowych błędów. Standardowe błędy są automatycznie obliczane w Bayesian LASSO ... Możesz bardzo łatwo wdrożyć Bayesian LASSO za pomocą schematu Gibbs Sampling ...

Bayesian LASSO potrzebuje wcześniejszych dystrybucji, które zostaną przypisane do parametrów modelu. W modelu LASSO mamy funkcję celu z jako parametr regularyzacji. Tutaj, ponieważ mamy -norm dla więc potrzebny jest do tego specjalny rodzaj wcześniejszej dystrybucji, rozkład LAPLACE skalowana mieszanina rozkładu normalnego z rozkładem wykładniczym jako gęstość mieszania. Na podstawie pełnego warunkowego posteriora każdego z parametrów należy wywnioskować.||yXβ||22+λ||β||1λ1β

Następnie można użyć Gibbs Sampling do symulacji łańcucha. Patrz Park & ​​Cassella (2008), „The Bayesian Lasso”, JASA , 103 , 482 .

Istnieją trzy nieodłączne wady LASSO:

  1. Trzeba wybrać metodą krzyżowej weryfikacji lub w inny sposób.λ

  2. Błędy standardowe są trudne do obliczenia, ponieważ LARS i inne algorytmy generują oszacowania punktowe dla .β

  3. Hierarchicznej struktury problemu nie da się zakodować za pomocą modelu częstościowego, co jest dość łatwe w ramach Bayesa.

Sandipan Karmakar
źródło
11

Aby dodać do powyższych odpowiedzi, wydaje się, że problem polega na tym, że nawet bootstrap jest prawdopodobnie niewystarczający, ponieważ oszacowanie z modelu ukaranego jest tendencyjne, a bootstrap będzie mówił tylko o wariancji - ignorując odchylenie oszacowania. Jest to ładnie podsumowane w winiecie dla ukaranego opakowania na stronie 18 .

Jeśli jednak jest używany do przewidywania, dlaczego wymagany jest standardowy błąd z modelu? Czy nie można odpowiednio sprawdzić poprawności lub uruchomić i wygenerować standardowy błąd wokół metryki związanej z prognozowaniem, takiej jak MSE?

B_Miner
źródło
3
Ładowanie początkowe może zarówno oszacować, jak i skorygować odchylenie, chociaż próbki muszą być dość duże.
Glen_b
3

W pakiecie R znajduje się pakiet selektywnych wniosków, https://cran.r-project.org/web/packages/selectiveInference/index.html , który zapewnia przedziały ufności i wartości p dla współczynników dopasowanych przez LASSO na podstawie poniższego dokumentu :

Stephen Reid, Jerome Friedman i Rob Tibshirani (2014). Badanie estymacji wariancji błędów w regresji lasso. arXiv: 1311.5274

PS: po prostu zdaj sobie sprawę, że to daje oszacowania błędów dla twoich parametrów, nie jestem pewien błędu w ostatecznej prognozie, jeśli to jest to, czego szukasz ... Przypuszczam, że możesz użyć „przedziałów prognoz populacji” , jeśli chcesz (przez parametry ponownego próbkowania zgodnie z dopasowaniem po wielowymiarowym rozkładzie normalnym).

Tom Wenseleers
źródło