Warum gibt mir die glmnet-Grat-Regression eine andere Antwort als die manuelle Berechnung?

28

Ich verwende glmnet, um Schätzungen der Gratregression zu berechnen. Ich habe einige Ergebnisse erhalten, die mich misstrauisch gemacht haben, dass glmnet wirklich das tut, was ich denke, dass es tut. Um dies zu überprüfen, habe ich ein einfaches R-Skript geschrieben, in dem ich das Ergebnis der von solve durchgeführten Ridge-Regression mit dem in glmnet vergleiche. Der Unterschied ist signifikant:

n    <- 1000
p.   <-  100
X.   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

beta1 <- solve(t(X)%*%X+5*diag(p),t(X)%*%Y)
beta2 <- glmnet(X,Y, alpha=0, lambda=10, intercept=FALSE, standardize=FALSE, 
                family="gaussian")$beta@x
beta1-beta2

Die Norm der Differenz liegt normalerweise bei 20, was nicht an numerisch unterschiedlichen Algorithmen liegen kann. Ich muss etwas falsch machen. Welche Einstellungen muss ich vornehmen glmnet, um das gleiche Ergebnis wie bei Ridge zu erzielen?

John
quelle
1
Hast du diese Frage gesehen ?
cdeterman
1
Ja, aber mit der Normalisierung erhalte ich immer noch nicht das gleiche Ergebnis.
John
Könnten Sie dann Ihren Code posten?
Shadowtalker
Ich hatte gerade das gleiche Problem! a = Datenrahmen (a = Jitter (1:10), b = Jitter (1:10), c = Jitter (1:10), d = Jitter (1:10), e = Jitter (1:10) , f = Jitter (1:10), g = Probe (Jitter (1:10)), y = Seq (10, 100, 10)); coef (1m.ridge (y ~ a + b + c + d + e + f + g, a, Lambda = 2,57)); coef (glmnet (as.matrix (a [, 1: 7]), a $ y, family = "gaussian", alpha = 0, lambda = 2,57 / 10)) Die Ergebnisse weichen stark voneinander ab und ähneln sich erheblich, wenn Ich benutze viel höhere Lambdas für glmnet.
a11msp
Faszinierend. Die Koeffizienten scheinen sich ungefähr um den Faktor 10 zu unterscheiden.
tomka

Antworten:

27

Der Unterschied, den Sie beobachten, beruht auf der zusätzlichen Division durch die Anzahl der Beobachtungen, N, die GLMNET in ihrer objektiven Funktion verwendet, und der impliziten Standardisierung von Y durch seine Standardabweichung, wie unten gezeigt.

12NysyXβ22+λβ22/2

wo wir anstelle von für , 1 / ( n - 1 ) s y s y = i ( y i - ˉ y ) 21/n1/(n1)sy

sy=i(yiy¯)2n

Indem Sie in Bezug auf Beta unterscheiden und die Gleichung auf Null setzen,

XTXβXTysy+Nλβ=0

Wenn wir nach Beta suchen, erhalten wir die Schätzung,

β~GLMNET=(XTX+NλIp)1XTysy

Um die Schätzungen (und die entsprechenden Strafen) für die ursprüngliche Metrik von Y wiederherzustellen, multipliziert GLMNET sowohl die Schätzungen als auch die mit und gibt diese Ergebnisse an den Benutzer zurück.sy

& lgr;unstd. =syλ

β^GLMNET=syβ~GLMNET=(XTX+NλIp)1XTy
λunstd.=syλ

Vergleichen Sie diese Lösung mit der Standardableitung der Gratregression.

β^=(XTX+λIp)1XTy

Beachten Sie, dass um einen zusätzlichen Faktor von N skaliert wird. Wenn wir die Funktion oder verwenden, wird die Strafe implizit um skaliert . Das heißt, wenn wir diese Funktionen verwenden, um die Koeffizientenschätzungen für ein , erhalten wir effektiv Schätzungen für .1 / s y λ λ = λ / s yλpredict()coef()1/syλλ=λ/sy

Basierend auf diesen Beobachtungen muss die in GLMNET verwendete Strafe mit einem Faktor von skaliert werden .sy/N

set.seed(123)

n    <- 1000
p   <-  100
X   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

beta1 <- solve(t(X)%*%X+10*diag(p),t(X)%*%(Y))[,1]

fit_glmnet <- glmnet(X,Y, alpha=0, standardize = F, intercept = FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

           [,1]        [,2]
[1,]  0.23793862  0.23793862
[2,]  1.81859695  1.81859695
[3,] -0.06000195 -0.06000195
[4,] -0.04958695 -0.04958695
[5,]  0.41870613  0.41870613
[6,]  1.30244151  1.30244151
[7,]  0.06566168  0.06566168
[8,]  0.44634038  0.44634038
[9,]  0.86477108  0.86477108
[10,] -2.47535340 -2.47535340

Die Ergebnisse verallgemeinern sich auf die Einbeziehung eines Abschnitts und standardisierter X-Variablen. Wir modifizieren eine standardisierte X-Matrix so, dass sie eine Spalte mit Einsen enthält, und die diagonale Matrix, dass sie an der Position [1,1] einen zusätzlichen Null-Eintrag aufweist (dh den Achsenabschnitt nicht benachteiligen). Sie können dann die Schätzungen anhand der jeweiligen Standardabweichungen der Stichprobe dekomprimieren (stellen Sie erneut sicher, dass Sie 1 / n verwenden, wenn Sie die Standardabweichung berechnen).

β^j=βj~sxj

β^0=β0~x¯Tβ^
mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)
X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}
X_scaled_ones <- cbind(rep(1,n), X_scaled)

beta3 <- solve(t(X_scaled_ones)%*%X_scaled_ones+1000*diag(x = c(0, rep(1,p))),t(X_scaled_ones)%*%(Y))[,1]
beta3 <- c(beta3[1] - crossprod(mean_x,beta3[-1]/sd_x), beta3[-1]/sd_x)

fit_glmnet2 <- glmnet(X,Y, alpha=0, thresh = 1e-20)
beta4 <- as.vector(coef(fit_glmnet2, s = sd_y*1000/n, exact = TRUE))

cbind(beta3[1:10], beta4[1:10])
             [,1]        [,2]
 [1,]  0.24534485  0.24534485
 [2,]  0.17661130  0.17661130
 [3,]  0.86993230  0.86993230
 [4,] -0.12449217 -0.12449217
 [5,] -0.06410361 -0.06410361
 [6,]  0.17568987  0.17568987
 [7,]  0.59773230  0.59773230
 [8,]  0.06594704  0.06594704
 [9,]  0.22860655  0.22860655
[10,]  0.33254206  0.33254206

Code hinzugefügt, um standardisiertes X ohne Intercept anzuzeigen:

set.seed(123)

n <- 1000
p <-  100
X <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)

X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}

beta1 <- solve(t(X_scaled)%*%X_scaled+10*diag(p),t(X_scaled)%*%(Y))[,1]

fit_glmnet <- glmnet(X_scaled,Y, alpha=0, standardize = F, intercept = 
FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

             [,1]        [,2]
 [1,]  0.23560948  0.23560948
 [2,]  1.83469846  1.83469846
 [3,] -0.05827086 -0.05827086
 [4,] -0.04927314 -0.04927314
 [5,]  0.41871870  0.41871870
 [6,]  1.28969361  1.28969361
 [7,]  0.06552927  0.06552927
 [8,]  0.44576008  0.44576008
 [9,]  0.90156795  0.90156795
[10,] -2.43163420 -2.43163420
Skijunkie
quelle
3
+6. Willkommen bei CV und vielen Dank für die übersichtliche Beantwortung dieser alten Frage.
Amöbe sagt Reinstate Monica
1
Es sollte die Identitätsmatrix anstelle von in der Lösung von , richtig? ˜ βββ~
user1769197
Ich stelle auch fest, dass Sie im zweiten Teil gesagt haben: "Die Ergebnisse verallgemeinern sich auf die Einbeziehung eines Abschnitts und standardisierter X-Variablen." Wenn Sie für diesen Teil den Achsenabschnitt ausschließen, weichen die Ergebnisse von glmnet nach denselben Berechnungen von der manuellen Berechnung ab.
user1769197
Richtig, ich habe die Lösung nach Bedarf mit der Identitätsmatrix anstelle von aktualisiert . Ich habe die Lösung für standardisiertes X ohne Unterbrechung überprüft und erhalte dennoch identische Ergebnisse (siehe zusätzlichen Code oben). β
Skijunkie
3

Nach https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html , wenn die Familie ist gaussian, glmnet()sollte minimieren

(1)12ni=1n(yiβ0xiTβ)2+λj=1p(α|βj|+(1α)βj2/2).

Wenn glmnet(x, y, alpha=1)das Lasso mit den Spalten in standardisiert angepasst wird, ist die Lösung für die gemeldete Strafe die Lösung zum Minimieren von Zumindest bei der Anpassung an die Gratregression ist die Lösung für eine gemeldete Strafe jedoch die Lösung zum Minimieren von Dabei ist die Standardabweichung von . Hier sollte die Strafe als .xλ

12ni=1n(yiβ0xiTβ)2+λj=1p|βj|.
glmnet_2.0-13glmnet(x, y, alpha=0)λ
12ni=1n(yiβ0xiTβ)2+λ12syj=1pβj2.
syyλ/sy

Was passieren könnte, ist, dass die Funktion zuerst zu standardisiert und dann minimiert. was effektiv minimieren soll oder äquivalent, um zu minimieren yy0

(2)12ni=1n(y0ixiTγ)2+ηj=1p(α|γj|+(1α)γj2/2),
12nsy2i=1n(yiβ0xiTβ)2+ηαsyj=1p|βj|+η1α2sy2j=1pβj2,
12ni=1n(yiβ0xiTβ)2+ηsyαj=1p|βj|+η(1α)j=1pβj2/2.

Für das Lasso ( ) ist es sinnvoll , zurück zu skalieren , um die Strafe als zu melden . Dann gilt für alle , hat als Strafe gemeldet werden Kontinuität der Ergebnisse über halten . Dies ist wahrscheinlich die Ursache des obigen Problems. Dies ist teilweise auf die Verwendung von (2) zur Lösung von (1) zurückzuführen. Nur wenn oder gibt es eine gewisse Äquivalenz zwischen den Problemen (1) und (2) (dh eine Entsprechung zwischen dem in (1) und dem in (2)). Für jedes andere& eegr; & eegr; s y α & eegr; s y α α = 0 α = 1 & lgr; & eegr; α ( 0 , 1 ) & lgr; & eegr;α=1ηηsyαηsyαα=0α=1ληα(0,1)sind die Probleme (1) und (2) zwei verschiedene Optimierungsprobleme, und es gibt keine Eins-zu-Eins-Entsprechung zwischen dem in (1) und dem in (2).λη

Chun Li
quelle
1
Ich kann nicht sehen, wo sich Ihre Antwort von der vorherigen unterscheidet. Könnten Sie das bitte erklären?
Firebug
1
@Firebug Ich wollte Aufschluss darüber geben, warum die Funktion das Lambda auf diese Weise meldet. Dies erscheint unnatürlich, wenn man es nur aus der Perspektive der Gratregression betrachtet, ist aber sinnvoll (oder muss dies auch sein), wenn man es aus der Perspektive des gesamten Spektrums betrachtet einschließlich Grat und Lasso.
Chun Li