R - Lasso-Regression - unterschiedliche Lambda pro Regressor

11

Ich möchte Folgendes tun:

1) OLS-Regression (kein Bestrafungsterm), um Beta-Koeffizienten ; steht für die zur Regression verwendeten Variablen. Ich mache das durch jbjj

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

2) Lasso-Regression mit einem Bestrafungsbegriff, die Auswahlkriterien sind die Bayesian Information Criteria (BIC), angegeben durch

λj=log(T)T|bj|

Dabei steht für die Variable / Regressornummer, für die Anzahl der Beobachtungen und für die in Schritt 1) ​​erhaltenen Anfangsbetas. Ich möchte Regressionsergebnisse für diesen bestimmten Wert haben, der für jeden verwendeten Regressor unterschiedlich ist. Wenn es also drei Variablen gibt, gibt es drei verschiedene Werte .T b j λ j λ jjTbjλjλj

Das OLS-Lasso-Optimierungsproblem ist dann gegeben durch

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

Wie kann ich dies in R entweder mit dem Lars- oder dem glmnet-Paket tun? Ich kann kein Lambda angeben und bin mir nicht 100% sicher, ob ich beim Ausführen die richtigen Ergebnisse erhalte

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

Ich freue mich über jede Hilfe hier.


Aktualisieren:

Ich habe jetzt den folgenden Code verwendet:

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)

In Zeile 1 verwende ich eine Kreuzvalidierung mit meinem angegebenen Straffaktor ( ), die für jeden Regressor unterschiedlich ist . Zeile 2 wählt das "Lambda.min" von fit.cv aus. Dies ist das Lambda, das den minimalen mittleren Kreuzvalidierungsfehler ergibt. Zeile 3 führt eine Lasso-Anpassung ( ) für die Daten durch. Wieder habe ich den Straffaktor . Zeile 4 extrahiert die Koeffizienten aus Anpassungen, die zu dem in Zeile 2 gewählten "optimalen" gehören .λj=log(T)T|bj|alpha=1λλ

Jetzt habe ich die Beta-Koeffizienten für die Regressoren, die die optimale Lösung des Minimierungsproblems darstellen

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

mit einem Straffaktor . Der optimale Satz von Koeffizienten ist höchstwahrscheinlich eine Teilmenge der Regressoren, die ich ursprünglich verwendet habe. Dies ist eine Folge der Lasso-Methode, die die Anzahl der verwendeten Regressoren verringert.λj=log(T)T|bj|

Ist mein Verständnis und der Code korrekt?

Dom
quelle
2
Sie können LATEX-Markup in Ihrem Beitrag verwenden, der in Dollarzeichen eingeschlossen ist. $\alpha$wird . Bitte machen Sie dies, da die Leute dadurch Ihre Frage leichter verstehen und daher beantworten können. α
Sycorax sagt Reinstate Monica

Antworten:

15

Aus der glmnetDokumentation ( ?glmnet) geht hervor, dass es möglich ist, eine Differenzialschrumpfung durchzuführen. Dies bringt uns zumindest teilweise dazu, die Frage von OP zu beantworten.

penalty.factor: Für jeden Koeffizienten können separate Straffaktoren angewendet werden. Dies ist eine Zahl, die multipliziert wird lambda, um eine unterschiedliche Schrumpfung zu ermöglichen. Kann für einige Variablen 0 sein, was keine Schrumpfung impliziert, und diese Variable ist immer im Modell enthalten. Der Standardwert ist 1 für alle Variablen (und implizit unendlich für die in aufgeführten Variablen exclude). Hinweis: Die Straffaktoren werden intern neu skaliert, um die Summe zu ergeben nvars, und die lambdaReihenfolge spiegelt diese Änderung wider.

Um die Frage vollständig zu beantworten, stehen Ihnen meines Erachtens zwei Ansätze zur Verfügung, je nachdem, was Sie erreichen möchten.

  1. Ihre Frage ist, wie Sie die Differentialschrumpfung anwenden glmnetund die Koeffizienten für einen bestimmten Wert abrufen können . Wenn einige Werte nicht 1 sind, wird bei jedem Wert von λ eine unterschiedliche Schrumpfung erreicht . Um eine Schrumpfung zu erreichen, beträgt die Schrumpfung für jedes b j ϕ j = log T.λpenalty.factorλbjWir müssen nur etwas Algebra machen. Seiϕjder Straffaktor fürbj, was geliefert werden würde. Aus der Dokumentation ist ersichtlich, dass diese Werte um einen Faktor vonCϕj=ϕ ' j stm=C m j = 1 logTneu skaliert werdenϕj=logTT|bj|ϕjbjpenalty.factorCϕj=ϕj. Dies bedeutetdassφ ' j ersetztφjin der unten Optimierung Ausdruck. Lösen Sie also nachC, geben Sie die Werteϕj anund extrahieren Sie dann die Koeffizienten fürλ=1. Ich würde empfehlen, zu verwenden.m=Cj=1mlogTT|bj|ϕjϕjCϕjglmnetλ=1coef(model, s=1, exact=T)

  2. Die zweite ist die "Standard" -Verwendung glmnet: Man führt eine wiederholte fache Kreuzvalidierung durch, um λ so auszuwählen , dass Sie die MSE außerhalb der Stichprobe minimieren. Dies ist, was ich unten ausführlicher beschreibe. Der Grund, warum wir CV verwenden und MSE außerhalb der Stichprobe prüfen , ist, dass MSE innerhalb der Stichprobe immer für λ = 0 minimiert wird , dh b ist eine gewöhnliche MLE. Die Verwendung von CV unter Variation von λ ermöglicht es uns, die Leistung des Modells bei Daten außerhalb der Stichprobe abzuschätzen und ein λ auszuwählen , das (in einem bestimmten Sinne) optimal ist.kλλ=0bλλ

Dieser glmnetAufruf gibt kein (und sollte es auch nicht, da er aus Leistungsgründen standardmäßig die gesamte λ- Trajektorie berechnet ). gibt die Koeffizienten für den λ- Wert zurück . Unabhängig von der von Ihnen angegebenen Auswahl von λ spiegelt das Ergebnis die Differenzstrafe wider, die Sie im Aufruf angewendet haben, um das Modell anzupassen.λλcoef(fits,s=something)λsomethingλ

Die Standardmethode zum Auswählen eines optimalen Werts von ist die Verwendung von anstelle von . Die Kreuzvalidierung wird verwendet, um das Ausmaß der Schrumpfung auszuwählen, das den Fehler außerhalb der Stichprobe minimiert, während die Spezifikation von einige Merkmale gemäß Ihrem Gewichtungsschema stärker schrumpft als andere.λcv.glmnetglmnetpenalty.factor

Dieses Verfahren wird optimiert

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

ϕjjthpenalty.factorλϕjλϕλϕλbλ

Dies ist im Grunde die Motivation, glmnetwie ich es verstehe: die bestrafte Regression zu verwenden, um ein Regressionsmodell zu schätzen, das hinsichtlich seiner Leistung außerhalb der Stichprobe nicht allzu optimistisch ist. Wenn dies Ihr Ziel ist, ist dies vielleicht doch die richtige Methode für Sie.

Sycorax sagt Reinstate Monica
quelle
+1 Das ist richtig. Ich werde auch hinzufügen, dass die Regularisierung der Regression als Bayes-Prior angesehen werden kann, dh Maximum a posteriori (MAP) ist regulierte Maximum Liklihood (ML). Wenn Sie in diesem Rahmen arbeiten, erhalten Sie mehr Flexibilität bei der Regularisierung, falls dies erforderlich sein sollte.
TLJ
Wenn ich laufe, pnlty = log(24)/(24*betas); fits = glmnet(x,y, alpha=1, intercept=FALSE, penalty.factor = pnlty) wie extrahiere ich dann die Regressor-Betas, die dem von mir angegebenen Lambda entsprechen, da das Lambda für jeden Risikofaktor unterschiedlich ist?
Dom
1
@ Dom Es wurde mir etwas zu spät klar, dass es einen offensichtlichen Weg gibt, genau das zu bekommen, was Sie wollen glmnet. Siehe meine überarbeitete Antwort.
Sycorax sagt Reinstate Monica
2
Passen Sie die Strafe nicht für jeden Prädiktor separat an. Dies würde in einigen Fällen nichts anderes als eine schrittweise Variablenauswahl bedeuten. Die bestrafte Regression verringert den mittleren quadratischen Fehler, indem eine sehr begrenzte Anzahl von Strafparametern angenommen und Informationen über Prädiktoren hinweg ausgeliehen werden.
Frank Harrell
2
@FrankHarrell Danke für den Kommentar! Es scheint, dass die Verwendung unterschiedlicher Strafen für jeden Prädiktor ein Bayes'sches Modell darstellt, das für jeden Parameter einen anderen Prior annimmt. Das scheint mir keine einzigartige Gefahr für die Bayes'sche Folgerung im Allgemeinen zu sein. Könnten Sie auch erläutern, wie die bestrafte Regression Informationen über Prädiktoren hinweg ausleiht? Ich bin mir nicht sicher, ob ich genau verstehe, wie das in einem solchen Szenario der Fall ist.
Sycorax sagt Reinstate Monica