Jak interpretować glmnet?

Próbuję dopasować wielowymiarowy model regresji liniowej z około 60 zmiennymi predykcyjnymi i 30 obserwacjami, więc używam pakietu glmnet do regresji regularnej, ponieważ p> n.

Przeglądałem dokumentację i inne pytania, ale nadal nie mogę zinterpretować wyników, oto przykładowy kod (z 20 predyktorami i 10 obserwacjami w celu uproszczenia):

Tworzę macierz x z num rzędów = num obserwacji i num cols = num predyktorów i wektor y reprezentujący zmienną odpowiedzi

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

Dopasowuję model glmnet, pozostawiając alpha jako domyślną (= 1 dla kary Lasso)

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

Rozumiem, że otrzymuję różne prognozy ze spadającymi wartościami lambda (tj. Kara)

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

Teraz przewiduję moje wartości Beta, wybierając na przykład najmniejszą podaną wartość lambda 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

Jeśli zamiast tego wybieram lambda z

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

Wszystkie zmienne to (.).

Wątpliwości i pytania:

1. Nie jestem pewien, jak wybrać lambda.
2. Czy należy użyć zmiennych innych niż (.), Aby dopasować do innego modelu? W moim przypadku chciałbym zachować jak najwięcej zmiennych.
3. Jak poznać wartość p, tj. Które zmienne znacznie przewidują odpowiedź?

Przepraszam za słabą wiedzę statystyczną! I dziękuję za wszelką pomoc.

Oto nieintuicyjny fakt - tak naprawdę nie powinieneś dawać glmnetowi jednej wartości lambda. Z dokumentacji tutaj :

Nie podawaj pojedynczej wartości lambda (zamiast prognoz dla CV użyj predykcji ()). Zamiast tego podaj malejącą sekwencję wartości lambda. glmnet polega na swoich początkowych prędkościach na szybkość, a często jest szybszy, aby znaleźć całą ścieżkę niż obliczyć pojedynczy atak.

cv.glmnetpomoże ci wybrać lambda, jak wspomniałeś w swoich przykładach. Autorzy pakietu glmnet sugerują cv$lambda.1sezamiast tego cv$lambda.min, ale w praktyce odnoszę sukces z tym ostatnim.

Po uruchomieniu cv.glmnet nie musisz ponownie uruchamiać glmnet! Każda lambda w grid ( cv$lambda) została już uruchomiona. Ta technika nazywa się „Warm Start” i możesz przeczytać więcej na jej temat tutaj . Parafrazując od wprowadzenia, technika Warm Start skraca czas działania metod iteracyjnych poprzez zastosowanie rozwiązania innego problemu optymalizacji (np. Glmnet z większą lambda) jako wartości początkowej dla późniejszego problemu optymalizacji (np. Glmnet z mniejszą lambda ).

Aby wyodrębnić pożądany przebieg cv.glmnet.fit, spróbuj tego:

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

Wersja (28.01.2017)

Nie trzeba hakować obiektu glmnet, tak jak ja powyżej; skorzystaj z porady @ alex23lemm poniżej i przekaż s = "lambda.min", s = "lambda.1se"lub inny numer (np. s = .007) do obu coefi predict. Zauważ, że twoje współczynniki i prognozy zależą od tej wartości, która jest ustalana przez walidację krzyżową. Użyj materiału siewnego, aby uzyskać powtarzalność! I nie zapominajcie, że jeśli nie dostarczyć "s"w coefa predict, będziesz przy użyciu domyślnego s = "lambda.1se". Rozgrzałem się do tego domyślnego ustawienia, gdy zobaczyłem, że działa lepiej w sytuacji małych danych.s = "lambda.1se"zapewnia również większą regularyzację, więc jeśli pracujesz z alfa> 0, będzie również dążyć do bardziej oszczędnego modelu. Możesz także wybrać wartość liczbową s za pomocą plot.glmnet, aby dostać się gdzieś pomiędzy (nie zapomnij jednak potęgować wartości z osi x!).

P1) Nie jestem pewien, jak wybrać lambda. Q2) Czy powinienem używać zmiennych innych niż (.) W celu dopasowania do innego modelu? W moim przypadku chciałbym zachować jak najwięcej zmiennych.

Zgodnie ze świetną odpowiedzią @ BenOgorek, zazwyczaj pozwalasz dopasowaniu na użycie całej sekwencji lambda, a następnie przy wydobywaniu optymalnych współczynników użyj wartości lambda.1se (w przeciwieństwie do tego, co zrobiłeś).

Dopóki przestrzegasz trzech ostrzeżeń poniżej, nie walcz z regularyzacją ani nie poprawiaj modelu: jeśli zmienna została pominięta, to dlatego, że dawała niższą ogólną karę. Ostrzeżenia są następujące:

1. Aby znormalizowane współczynniki były znaczące, upewnij się, że wcześniej wyraźnie znormalizowałeś średnią i stdev zmiennej scale(); nie polegaj na glmnet(standardize=T). W celu uzasadnienia patrz: Czy normalizacja przed Lasso jest naprawdę konieczna? ; w zasadzie zmienna o dużych wartościach może zostać niesprawiedliwie ukarana w regularyzacji.

2. Aby zapewnić powtarzalność, uruchom z set.seedkilkoma losowymi nasionami i sprawdź stabilizowane współczynniki.

3. Jeśli chcesz mniej surowej regularyzacji, tj. Uwzględniono więcej zmiennych, użyj alfa <1 (tj. Odpowiedniej elastycznej siatki) zamiast zwykłego grzbietu. Sugeruję zamiatanie alfa od 0 do 1. Jeśli masz zamiar to zrobić, to aby uniknąć nadmiernego dopasowania hiperparametru alfa i błędu regresji, musisz użyć walidacji krzyżowej, tj. Użyj cv.glmnet()raczej niż prostej 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...
}

Jeśli chcesz zautomatyzować takie wyszukiwanie siatki za pomocą CV, możesz albo sam go zakodować, albo użyć pakietu caret na glmnet; Caret robi to dobrze. Dla cv.glmnet nfoldswartości parametru, pick 3 (minimum), jeśli zbiór danych jest mała, lub 5 lub 10, jeśli jest duża.

P3) Jak poznać wartość p, tj. Które zmienne znacząco przewidują odpowiedź?

Nie, nie mają znaczenia . Jak wyjaśniono szczegółowo w Dlaczego nie jest wskazane uzyskiwanie statystycznych informacji podsumowujących dla współczynników regresji z modelu glmnet?

Po prostu pozwól cv.glmnet()dokonać wyboru zmiennej automatycznie. Z zastrzeżeniami powyżej. I oczywiście rozkład zmiennej odpowiedzi powinien być normalny (zakładając, że używasz family='gaussian').