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 ridge-regression glmnet Jan
źródło

Widziałeś to pytanie ?

cdeterman,

Tak, ale nadal nie otrzymuję tego samego wyniku przy użyciu normalizacji.

Jan

Czy mógłbyś wtedy opublikować swój kod?

shadowtalker,

Właśnie miałem ten sam problem! a = ramka danych (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 = próbka (jitter (1:10)), y = seq (10 100,10)); coef (lm.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)) Wyniki różnią się nieco i stają się znacznie bardziej podobne, gdy Do glmnet używam znacznie wyższych lambd.

a11msp,

Intrygancki. Współczynniki wydają się różnić z grubsza o współczynnik 10.

tomka

Odpowiedzi:

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

gdzie używamy zamiast dla , $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}$

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

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

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

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

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

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{β}}_{G L M N E T} = s_{y} {\tilde{β}}_{G L M N E T.} = (X^{T.} X + N. λ {ja}_{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$

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

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

$\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 . 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{β}}_{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

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

skijunkie
źródło

+6. Witamy w CV i dziękuję za udzielenie odpowiedzi na to stare pytanie w tak jasny sposób.

Ameba mówi Przywróć Monikę

β

$\beta$

\tilde{β}

$\tilde{\beta}$

Zauważam również, że w drugiej części, w której powiedzieliście: „Wyniki uogólniają się, dodając przecięcie i znormalizowane zmienne X”; dla tej części, jeśli wykluczysz przechwytywanie, a następnie po tych samych obliczeniach wyniki glmnet staną się inne niż obliczenia ręczne.

user1769197

β

$\beta$

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

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

glmnet(x, y, alpha=1) $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$

$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,$ or equivalently, to minimize

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

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

Chun Li
źródło

Nie widzę, czym różni się twoja odpowiedź od poprzedniej. Czy mógłbyś wyjaśnić?

Firebug

@Firebug Chciałem wyjaśnić, dlaczego funkcja zgłasza lambda w ten sposób, co wydaje się nienaturalne, gdy patrzy się wyłącznie z perspektywy regresji grzbietu, ale ma sens (lub musi być w ten sposób), gdy patrzy się z perspektywy całego spektrum w tym zarówno grzbiet, jak i lasso.

Chun Li