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

28

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?

Jan
źródło
1
Widziałeś to pytanie ?
cdeterman,
1
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:

27

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.

12NysyXβ22+λβ22/2

gdzie używamy zamiast 1 / ( n - 1 ) dla s y , s y =1/n1/(n1)sy

sy=i(yiy¯)2n

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

XTXβXTysy+Nλβ=0

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

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

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,sy

λ

β^GLMNET=syβ~GLMNET.=(XT.X+N.λjap)1XTy
λunstd.=syλ

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

β^=(XTX+λIp)1XTy

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 λ λpredict()coef()1/syλ , faktycznie uzyskujemy oszacowania dla .λ=λ/sy

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

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

β^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

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
3
+6. Witamy w CV i dziękuję za udzielenie odpowiedzi na to stare pytanie w tak jasny sposób.
Ameba mówi Przywróć Monikę
1
ββ~
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
β
3

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

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

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

yy0

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

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

Chun Li
źródło
1
Nie widzę, czym różni się twoja odpowiedź od poprzedniej. Czy mógłbyś wyjaśnić?
Firebug
1
@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