Dlaczego regresja grzbietu glmnet daje mi inną odpowiedź niż ręczne obliczanie?

Używam glmnet do obliczania oszacowań regresji grzbietu. Mam wyniki, które wzbudziły we mnie podejrzenia, ponieważ glmnet naprawdę robi to, co myślę. Aby to sprawdzić, napisałem prosty skrypt R, w którym porównuję wynik regresji grzbietu wykonanej przez rozwiązanie i ten w glmnet, różnica jest znacząca:

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

Norma różnicy wynosi zwykle około 20, co nie może wynikać z różnych algorytmów numerycznych, muszę robić coś złego. Jakie ustawienia muszę ustawić glmnet, aby uzyskać taki sam efekt jak w przypadku kalenicy?

Różnica, którą obserwujesz, wynika z dodatkowego podziału przez liczbę obserwacji N, które GLMNET wykorzystuje w swojej funkcji celu i domyślnej standaryzacji Y przez odchylenie standardowe próbki, jak pokazano poniżej.

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

gdzie używamy zamiast 1 / ( n - 1 ) dla s y , s y$1/n$$1/(n-1)$$s_y$

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

Różniczkując względem beta, ustawiając równanie na zero,

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

I rozwiązując wersję beta, otrzymujemy oszacowanie,

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

Aby odzyskać dane szacunkowe (i odpowiadające im kary) na oryginale metryczny Y, GLMNET mnoży zarówno szacunki i lambdas by i powroty tych wyników do użytkownika,$s_y$

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

Porównaj to rozwiązanie ze standardową pochodną regresji kalenicy.

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

Zauważ, że jest skalowany przez dodatkowy czynnik N. Dodatkowo, gdy używamy lub funkcji, kara będzie domyślnie skalowane przez 1 / s y . To znaczy, kiedy używamy tych funkcji do uzyskania oszacowań współczynników dla niektórych λ $\lambda$predict()coef()$1/s_y$$\lambda^*$ , faktycznie uzyskujemy oszacowania dla .$\lambda=\lambda^*/s_y$

Na podstawie tych obserwacji, kara stosowana w GLMNET musi być skalowane przez czynnik .$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

Wyniki uogólniają się na włączenie przechwytywania i znormalizowanych zmiennych X. Modyfikujemy znormalizowaną macierz X, tak aby zawierała kolumnę jedynek, a macierz diagonalną, aby miała dodatkowy wpis zerowy w pozycji [1,1] (tj. Nie penalizować przecięcia). Następnie możesz unormować szacunki według odpowiednich odchyleń standardowych próbki (ponownie upewnij się, że używasz 1 / n przy obliczaniu odchylenia standardowego).

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

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

Dodano kod pokazujący znormalizowany X bez przechwytywania:

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

Według https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html , kiedy rodzina jest gaussian, glmnet()powinna zminimalizować

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

glmnet(x, y, alpha=1)$x$$\lambda$

$$\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}{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$$y$$\lambda/s_y$

$y$$y_0$

$$\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}{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,$$ or equivalently, to minimize $$\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.$$

For the lasso ($\alpha=1$), scaling $\eta$ back to report the penalty as $\eta s_y$ makes sense. Then for all $\alpha$, $\eta s_y$ has to be reported as the penalty to maintain continuity of the results across $\alpha$. This probably is the cause of the problem above. This is partly due to using (2) to solve (1). Only when $\alpha=0$ or $\alpha=1$ istnieje pewna równoważność między problemami (1) i (2) (tj. zgodność między $\lambda$ w (1) i $\eta$w 2)). Dla każdego innego$\alpha\in(0,1)$, problemy (1) i (2) to dwa różne problemy z optymalizacją i nie ma żadnej relacji jeden do jednego między $\lambda$ w (1) i $\eta$ w 2).