R - Regresja Lasso - inna Lambda na regresor

11

Chcę wykonać następujące czynności:

1) Regresja OLS (bez kary), aby uzyskać współczynniki beta ; oznacza zmienne użyte do regresji. Robię to przezbjj

lm.model = lm(y~ 0 + x)
betas    = coefficients(lm.model)

2) Regresja Lasso z terminem karnym, kryteriami wyboru są Bayesowskie Kryteria Informacyjne (BIC), podane przez

λj=log(T)T|bj|

gdzie oznacza numer zmiennej / regresora, oznacza liczbę obserwacji, a dla początkowych bet uzyskanych w kroku 1). Chcę uzyskać wyniki regresji dla tej konkretnej wartości , która jest inna dla każdego zastosowanego regresora. Dlatego jeśli są trzy zmienne, będą trzy różne wartości .jTbjλjλj

Problem optymalizacji OLS-Lasso jest następnie podawany przez

minbϵRn={t=1T(ytbXt)2+Tj=1m(λt|bj|)}

Jak mogę to zrobić w R z pakietem lars lub glmnet? Nie mogę znaleźć sposobu na określenie lambda i nie jestem w 100% pewien, czy otrzymam prawidłowe wyniki po uruchomieniu

lars.model <- lars(x,y,type = "lasso", intercept = FALSE)
predict.lars(lars.model, type="coefficients", mode="lambda")

Doceniam każdą pomoc tutaj.


Aktualizacja:

Użyłem teraz następującego kodu:

fits.cv = cv.glmnet(x,y,type="mse",penalty.factor = pnlty)
lmin    = as.numeric(fits.cv[9]) #lambda.min
fits    = glmnet(x,y, alpha=1, intercept=FALSE, penalty.factor = pnlty)
coef    = coef(fits, s = lmin)

W wierszu 1 używam weryfikacji krzyżowej z moim określonym współczynnikiem kary ( ), który jest inny dla każdego regresora . Wiersz 2 wybiera „lambda.min” z fits.cv, który jest lambda dającym minimalny średni błąd walidacji krzyżowej. Linia 3 wykonuje dopasowanie lasso ( ) na danych. Ponownie użyłem współczynnika kary . Wiersz 4 wyodrębnia współczynniki z dopasowań, które należą do „optymalnego” wybranego w wierszu 2.λλλj=log(T)T|bj|alpha=1λλ

Teraz mam współczynniki beta dla regresorów, które przedstawiają optymalne rozwiązanie problemu minimalizacji

minbϵRn={t=1T(ytbXt)2+Tj=1m(λt|bj|)}

z czynnikiem karnym . Optymalny zestaw współczynników jest najprawdopodobniej podzbiorem regresorów, które początkowo stosowałem, jest to konsekwencja metody Lasso, która zmniejsza liczbę używanych regresorów.λj=log(T)T|bj|

Czy moje rozumienie i kod są poprawne?

Dom
źródło
2
Możesz użyć znaczników LATEX we wpisie, ujętym w znaki dolara. $\alpha$staje się . Proszę, zrób to, ponieważ dzięki temu ludzie łatwiej zrozumieją twoje pytanie, a zatem odpowiedzą na nie. α
Sycorax mówi Przywróć Monikę

Odpowiedzi:

15

Z glmnetdokumentacji ( ?glmnet) wynika, że ​​możliwe jest wykonanie skurczu różnicowego. To pozwala nam przynajmniej częściowo odpowiedzieć na pytanie OP.

penalty.factor: Do każdego współczynnika można zastosować osobne współczynniki kary. Jest to liczba mnożąca się, lambdaaby umożliwić skurcz różnicowy. Może wynosić 0 dla niektórych zmiennych, co oznacza brak skurczu, a zmienna ta jest zawsze uwzględniana w modelu. Wartość domyślna to 1 dla wszystkich zmiennych (i domyślnie nieskończoność dla zmiennych wymienionych w exclude). Uwaga: współczynniki kar są wewnętrznie przeskalowywane do sumy nvars, a lambdasekwencja odzwierciedla tę zmianę.

Aby jednak w pełni odpowiedzieć na pytanie, uważam, że dostępne są dwa podejścia, w zależności od tego, co chcesz osiągnąć.

  1. Twoje pytanie brzmi: jak zastosować kurczenie różnicowe glmneti odzyskać współczynniki dla określonej wartości . Podanie st niektórych wartości nie 1 powoduje różnicowy skurcz przy dowolnej wartości . Aby osiągnąć skurcz, skurcz dla każdego wynosi , po prostu musimy wykonać algebrę. Niech będzie czynnikiem karnym dla , co zostanie dostarczone . Z dokumentacji wynika, że ​​te wartości są ponownie skalowane o współczynnik st . Oznacza to, żeλpenalty.factorλbjϕj=logTT|bj|ϕjbjpenalty.factorCϕj=ϕjm=Cj=1mlogTT|bj|ϕjzastępuje w poniższym wyrażeniu optymalizacyjnym. więc dla , podaj wartości do , a następnie wyodrębnij współczynniki dla . Poleciłbym użyć .ϕjCϕjglmnetλ=1coef(model, s=1, exact=T)

  2. Drugi to „standardowy” sposób użycia glmnet: jeden przeprowadza wielokrotne sprawdzanie poprawności -krotnie, aby wybrać tak aby zminimalizować MSE poza próbą. To właśnie opisuję bardziej szczegółowo. Powodem, dla którego używamy CV i sprawdzamy MSE poza próbą, jest to, że MSE w próbce zawsze będzie zminimalizowane dla , tj. jest zwykłym MLE. Używanie CV podczas zmieniania pozwala nam oszacować wydajność modelu na danych poza próbą i wybrać optymalną (w pewnym sensie) .kλλ=0bλλ

To glmnetwywołanie nie określa (również nie powinno, ponieważ domyślnie oblicza całą trajektorię ze względu na wydajność). powróci współczynniki dla wartości . Ale bez względu na wybór który podasz, wynik będzie odzwierciedlał karę różnicową zastosowaną w wezwaniu do dopasowania modelu.λ λ λλλcoef(fits,s=something)λsomethingλ

Standardowym sposobem wyboru optymalnej wartości jest użycie zamiast . Walidacja krzyżowa służy do wyboru stopnia skurczu, który minimalizuje błąd poza próbą, podczas gdy specyfikacja skurczy niektóre funkcje bardziej niż inne, zgodnie z twoim schematem ważenia.λcv.glmnetglmnetpenalty.factor

Ta procedura jest optymalizowana

minbRmt=1T(ytbXt)2+λj=1m(ϕj|bj|)

gdzie jest czynnikiem karnym dla funkcji (co podajesz w argumencie). (Różni się to nieco od wyrażenia optymalizacyjnego; zwróć uwagę, że niektóre indeksy dolne są różne.) Zauważ, że termin jest taki sam we wszystkich funkcjach, więc jedynym sposobem, aby niektóre funkcje były zmniejszone bardziej niż inne, jest . Co ważne, i nie są takie same; to skalar, a to wektor! W tym wyrażeniu jest stałe / zakłada się, że jest znane; to znaczy optymalizacja wybierze optymalne , a nie optymalneϕjjthpenalty.factorλϕjλϕλϕλbλ.

Jest to w zasadzie motywacja, glmnetjak rozumiem: stosowanie regresji karnej w celu oszacowania modelu regresji, który nie jest nadmiernie optymistyczny w odniesieniu do jego wydajności poza próbą. Jeśli to jest twój cel, być może jest to w końcu właściwa metoda.

Sycorax mówi Przywróć Monikę
źródło
+1 To jest poprawne. Dodam również, że regularyzacja regresji może być postrzegana jako wcześniejsze bayesowskie, tj. Maksimum a posteriori (MAP) jest uregulowane maksymalnym prawdopodobieństwem (ML). Praca w tych ramach zapewnia sobie większą elastyczność regularyzacji, jeśli zajdzie taka potrzeba.
TLJ
Jeśli uruchomię, pnlty = log(24)/(24*betas); fits = glmnet(x,y, alpha=1, intercept=FALSE, penalty.factor = pnlty) jak wyodrębnić beta regresora, które odpowiadają określonej przeze mnie lambdzie, ponieważ lambda jest inna dla każdego czynnika ryzyka?
Dom
1
@Dom Przyszło mi do głowy trochę za późno, że istnieje oczywisty sposób na uzyskanie dokładnie tego, czego chcesz glmnet. Zobacz moją poprawioną odpowiedź.
Sycorax mówi Przywróć Monikę
2
Uważaj na dostosowanie kary osobno dla każdego predyktora. W niektórych przypadkach oznaczałoby to nic innego jak stopniowy wybór zmiennych. Regresja karana zmniejsza średni błąd kwadratu, zakładając bardzo ograniczoną liczbę parametrów kary i informacje o pożyczkach między predyktorami.
Frank Harrell,
2
@FrankHarrell Dzięki za komentarz! Wydaje się, że zastosowanie różnych kar dla każdego predyktora jest równoznaczne z modelem bayesowskim, który zakłada inny uprzedni dla każdego parametru. Nie wydaje mi się to, że stanowi wyjątkowe zagrożenie w porównaniu z wnioskami bayesowskimi. Czy mógłbyś również wyjaśnić, w jaki sposób regresja karna pożycza informacje między predyktorami? Nie jestem pewien, czy w pełni rozumiem, jak to jest w takim scenariuszu.
Sycorax mówi Przywróć Monikę