Używanie LASSO z pakietu lars (lub glmnet) w R do wyboru zmiennych

39

Przepraszam, jeśli to pytanie jest trochę podstawowe.

Chciałbym użyć selekcji zmiennych LASSO dla modelu wielokrotnej regresji liniowej w R. Mam 15 predyktorów, z których jeden jest kategoryczny (czy to spowoduje problem?). Po ustawieniu mojego i Y używam następujące polecenia:xy

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

Mój problem polega na tym, kiedy używam coef(model). Zwraca macierz z 15 wierszami, z każdym dodatkowym predyktorem dodawanym za każdym razem. Jednak nie ma sugestii, który model wybrać. Czy coś przeoczyłem? Czy istnieje sposób, aby uzyskać pakiet Larsa, aby zwrócił tylko jeden „ najlepszy ” model?

Istnieją inne posty sugerujące użycie glmnetzamiast tego, ale wydaje się to bardziej skomplikowane. Próba jest następujący, przy użyciu tego samego i y . Czy coś tu przegapiłem ?: xy

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

Ostatnie polecenie zwraca listę moich zmiennych, większość ze współczynnikiem, chociaż niektóre mają wartość = 0. Czy to właściwy wybór „ najlepszego ” modelu wybranego przez LASSO? Jeśli następnie dopasuję model liniowy ze wszystkimi moimi zmiennymi, które miały współczynniki not=0, otrzymam bardzo podobne, ale nieco inne, oszacowania współczynników. Czy jest powód tej różnicy? Czy akceptowalny byłby remont modelu liniowego z tymi zmiennymi wybranymi przez LASSO i uznanie go za mój ostateczny model? W przeciwnym razie nie widzę żadnych wartości p dla istotności. Czy coś przeoczyłem?

Robi

type.gaussian="covariance" 

upewnić się, że glmnetużywa wielu regresji liniowej?

Czy automatyczna normalizacja zmiennych wpływa w ogóle na współczynniki? Czy jest jakiś sposób na uwzględnienie warunków interakcji w procedurze LASSO?

Chciałbym użyć tej procedury bardziej jako demonstracji tego, jak można użyć LASSO, niż w jakimkolwiek modelu, który faktycznie zostanie użyty do każdego ważnego wnioskowania / prognozy, jeśli to cokolwiek zmieni.

Dziękujemy za poświęcenie czasu na przeczytanie tego. Wszelkie ogólne uwagi na temat LASSO / lars / glmnet będą również mile widziane.

James
źródło
4
Na marginesie, jeśli chcesz interpretować wynik, zademonstruj, że zestaw zmiennych wybranych przez lasso jest stabilny. Można to zrobić za pomocą symulacji Monte Carlo lub przez załadowanie własnego zestawu danych.
Frank Harrell,

Odpowiedzi:

28

Korzystanie glmnetjest bardzo łatwe, gdy się go zrozumie, dzięki doskonałej winiecie w http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html (możesz również sprawdzić stronę pakietu CRAN). Jeśli chodzi o najlepszą lambda glmnet, zasadą jest użycie

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

zamiast lambda.min.

Aby zrobić to samo lars, musisz to zrobić ręcznie. Oto moje rozwiązanie

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,]

Pamiętaj, że nie jest to dokładnie to samo, ponieważ zatrzymuje się na węźle lasso (gdy wchodzi zmienna) zamiast w dowolnym punkcie.

Należy pamiętać, że glmnetjest to teraz preferowany pakiet, jest on aktywnie utrzymywany, bardziej niż lars, i że wcześniej pojawiły się pytania dotyczące odpowiedzi glmnetvs lars(stosowane algorytmy różnią się).

Jeśli chodzi o twoje pytanie o użycie lasso do wyboru zmiennych, a następnie dopasowania OLS, jest to ciągła debata. Google for OLS po Lasso i jest kilka artykułów na ten temat. Nawet autorzy Elements of Statistics Learning przyznają, że jest to możliwe.

Edycja : Oto kod do dokładniejszego odtwarzania tego, co glmnetdziałalars

  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),]
Juancentro
źródło
+1 Świetna odpowiedź! Czy ty lub ktokolwiek mógłby wyjaśnić, dlaczego lambda.1se jest regułą, zamiast lambda.min?
Erosennin
Po 4 latach pisania tego (i nie używając lasso przez chwilę) moja pamięć zniknęła. Przepraszam!
Juancentro
8

Powracam do tego pytania jakiś czas temu, ponieważ myślę, że rozwiązałem właściwe rozwiązanie.

Oto replika korzystająca z zestawu danych 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)']

„zmienne” daje listę zmiennych, które rozwiązują najlepsze rozwiązanie.

Jason
źródło
1
Szukałem kodu i stwierdziłem, że „testowanie” nie zostało jeszcze zdefiniowane, a zatem kod: „final.list <-testing [-removed] #removing zmienne” daje błąd: obiekt nie został znaleziony Spoglądając na kod I załóżmy, że zamiast używania „testowania” należy użyć „cp.list”, aby kod był: final.list <-cp.list [-removed] #usuning zmienne final.list <-c (final.list, duplikaty) #dodanie w tych zmiennych, które zostały usunięte, a następnie dodane później Daj mi znać, czy to jest poprawne Z poważaniem
3
`% ni%` <-Negate (`% ni%`); ## wygląda źle. Podczas gdy `% ni%` <-Negate (`% w%`); ## wygląda dobrze. Wydaje mi się, że formatyzator wymiany stosów popsuł to ...
Chris
Czy potrafisz opracować sposób wyboru parametrów nfolds=5i alpha=0.5?
colin,
7

Być może porównanie z regresją stopniową wyboru do przodu pomoże (patrz poniższy link do strony jednego z autorów http://www-stat.stanford.edu/~tibs/lasso/simple.html). Takie podejście zastosowano w rozdziale 3.4.4 elementów uczenia statystycznego (dostępne online za darmo). Pomyślałem, że rozdział 3.6 tej książki pomógł zrozumieć związek między najmniejszymi kwadratami, najlepszym podzbiorem i lasso (plus kilka innych procedur). Przydaje mi się również transpozycja współczynnika t (coef (model)) i write.csv, dzięki czemu mogę otworzyć go w programie Excel wraz z kopią wykresu (model) z boku. Możesz posortować według ostatniej kolumny, która zawiera oszacowanie najmniejszych kwadratów. Następnie możesz wyraźnie zobaczyć, w jaki sposób każda zmienna jest dodawana na każdym kroku i jak w rezultacie zmieniają się współczynniki. Oczywiście to nie jest cała historia, ale mam nadzieję, że to będzie początek.

Joel Cadwell
źródło
3

larsi glmnetdziałają na surowych matrycach. Aby uwzględnić warunki interakcji, musisz sam zbudować macierze. Oznacza to jedną kolumnę na interakcję (która jest na poziom na czynnik, jeśli masz czynniki). Sprawdź, lm()jak to robi (ostrzeżenie: są smoki).

Aby zrobić to teraz, zrób coś takiego: Aby ręcznie dokonać określenie interakcji, to może (ale może nie powinienem , bo to spowalnia) zrobić:

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

Następnie użyj tego w Larach (zakładając, że masz do czynienia z ykopaniem):

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

Chciałbym pomóc Ci w przypadku innych pytań. Znalazłem ten, ponieważ lars sprawia mi ból, a dokumentacja w nim i w Internecie jest bardzo cienka.

kousu
źródło
2
„Ostrzeżenie: są smoki” Z tym jest całkiem łatwo model.matrix().
Gregor
2

LARS rozwiązuje CAŁĄ ścieżkę rozwiązania. Ścieżka rozwiązania jest fragmentarycznie liniowa - istnieje skończona liczba punktów „wycięcia” (tj. Wartości parametru regularyzacji), w których zmienia się rozwiązanie.

Zatem macierz otrzymywanych rozwiązań to wszystkie możliwe rozwiązania. Na liście, którą zwraca, powinien również podać wartości parametru regularyzacji.

Adam
źródło
Dziękuję za Twoją odpowiedź. Czy istnieje sposób wyświetlania wartości parametru regularyzacji? Ponadto czy istnieje sposób, aby wybrać między rozwiązaniami opartymi na tym parametrze? (Również parametr lambda?)
James
Zauważ, że liniowość częściowa nie oznacza, że ​​linie są poziome, a zatem rozwiązanie zmienia się cały czas w przypadku lambda. Na przykład dla celów predykcyjnych można by uzyskać siatkę wartości lambda nie tylko w węzłach, ale także między nimi. Jest całkiem możliwe, że jakiś punkt między węzłami daje najlepszą lepszą wydajność predykcyjną.
Richard Hardy