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

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?

r ridge-regression glmnet John
quelle

Hast du diese Frage gesehen ?

cdeterman

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:

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.

\frac{1}{2 N} {‖ \frac{y}{s_{y}} - X β ‖}_{2}^{2} + λ ‖ β ‖_{2}^{2} / 2

$\frac{1}{2N}\left\|\frac{y}{s_y}-X\beta\right\|^2_{2}+\lambda\|\beta\|^2_{2}/2$

wo wir anstelle von für , $1/n$ $1/(n-1)$ $s_y$

s_{y} = \frac{\sum_{i} (y_{i} - \bar{y})^{2}}{n}

$s_y=\frac{\sum_i(y_i-\bar{y})^2}{n}$

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

X^{T} X β - \frac{X^{T} y}{s_{y}} + N λ β = 0

$X^TX\beta-\frac{X^Ty}{s_y}+N\lambda\beta =0$

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

{\tilde{β}}_{G L M N E T} = (X^{T} X + N λ I_{p})^{- 1} \frac{X^{T} y}{s_{y}}

$\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}\frac{X^Ty}{s_y}$

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. $s_y$

{\hat{β}}_{G L M N E T} = s_{y} {\tilde{β}}_{G L M N E T} = (X^{T} X + N λ I_{p})^{- 1} X^{T} y

$\hat{\beta}_{GLMNET}=s_y\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}X^Ty$

λ_{u n s t d .} = s_{y} λ

$\lambda_{unstd.}=s_y\lambda$

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

\hat{β} = (X^{T} X + λ I_{p})^{- 1} X^{T} y

$\hat{\beta}= (X^TX+\lambda I_p)^{-1}X^Ty$

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 . $\lambda$ predict()coef() $1/s_y$ $\lambda^*$ $\lambda=\lambda^*/s_y$

Basierend auf diesen Beobachtungen muss die in GLMNET verwendete Strafe mit einem Faktor von skaliert werden . $s_y/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).

{\hat{β}}_{j} = \frac{\tilde{β_{j}}}{s_{x_{j}}}

$\hat\beta_{j}=\frac{\tilde{\beta_j}}{s_{x_j}}$

{\hat{β}}_{0} = \tilde{β_{0}} - {\bar{x}}^{T} \hat{β}

$\hat\beta_{0}=\tilde{\beta_0}-\bar{x}^T\hat{\beta}$

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

+6. Willkommen bei CV und vielen Dank für die übersichtliche Beantwortung dieser alten Frage.

Amöbe sagt Reinstate Monica

Es sollte die Identitätsmatrix anstelle von in der Lösung von , richtig?

β

$\beta$

\tilde{β}

$\tilde{\beta}$

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).

β

$\beta$

Skijunkie

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

\begin{matrix} (1) & \frac{1}{2 n} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + λ \sum_{j = 1}^{p} (α | β_{j} | + (1 - α) β_{j}^{2} / 2) . \end{matrix}

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda\sum_{j=1}^p(\alpha|\beta_j| +(1-\alpha)\beta_j^2/2). \tag{1}$

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$ $\lambda$

\frac{1}{2 n} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + λ \sum_{j = 1}^{p} | β_{j} | .

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \sum_{j=1}^p |\beta_j|.$ glmnet_2.0-13glmnet(x, y, alpha=0)

λ

$\lambda$

\frac{1}{2 n} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + λ \frac{1}{2 s_{y}} \sum_{j = 1}^{p} β_{j}^{2} .

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \frac{1}{2s_y} \sum_{j=1}^p \beta_j^2.$

s_{y}

$s_y$

y

$y$

λ / s_{y}

$\lambda/s_y$

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

\begin{matrix} (2) & \frac{1}{2 n} \sum_{i = 1}^{n} (y_{0 i} - x_{i}^{T} γ)^{2} + η \sum_{j = 1}^{p} (α | γ_{j} | + (1 - α) γ_{j}^{2} / 2), \end{matrix}

$\frac{1}{2n} \sum_{i=1}^n (y_{0i}-x_i^T\gamma)^2 +\eta \sum_{j=1}^p(\alpha|\gamma_j| +(1-\alpha)\gamma_j^2/2), \tag{2}$

\frac{1}{2 n s_{y}^{2}} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + η \frac{α}{s_{y}} \sum_{j = 1}^{p} | β_{j} | + η \frac{1 - α}{2 s_{y}^{2}} \sum_{j = 1}^{p} β_{j}^{2},

$\frac{1}{2n s_y^2} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta \frac{\alpha}{s_y} \sum_{j=1}^p |\beta_j| +\eta \frac{1-\alpha}{2s_y^2} \sum_{j=1}^p \beta_j^2,$

\frac{1}{2 n} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + η s_{y} α \sum_{j = 1}^{p} | β_{j} | + η (1 - α) \sum_{j = 1}^{p} β_{j}^{2} / 2.

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta s_y \alpha \sum_{j=1}^p |\beta_j| +\eta (1-\alpha) \sum_{j=1}^p \beta_j^2/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 $\alpha=1$ $\eta$ $\eta s_y$ $\alpha$ $\eta s_y$ $\alpha$ $\alpha=0$ $\alpha=1$ $\lambda$ $\eta$ $\alpha\in(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). $\lambda$ $\eta$

Chun Li
quelle

Ich kann nicht sehen, wo sich Ihre Antwort von der vorherigen unterscheidet. Könnten Sie das bitte erklären?

Firebug

@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