Pourquoi la régression de la crête glmnet me donne-t-elle une réponse différente du calcul manuel?

28

J'utilise glmnet pour calculer les estimations de régression de crête. J'ai obtenu des résultats qui m'ont rendu suspect dans la mesure où glmnet fait vraiment ce que je pense qu'il fait. Pour vérifier cela, j'ai écrit un script R simple où je compare le résultat de la régression de crête effectuée par résoudre et celui de glmnet, la différence est significative:

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

La norme de la différence est généralement d'environ 20, ce qui ne peut pas être dû à des algorithmes numériquement différents, je dois faire quelque chose de mal. Quels sont les paramètres que je dois définir glmnetpour obtenir le même résultat qu'avec la crête?

John
la source
1
Avez-vous vu cette question ?
cdeterman du
1
Oui, mais je n'obtiens toujours pas le même résultat en utilisant la normalisation.
John
Pourriez-vous poster votre code alors?
shadowtalker
Je viens d'avoir le même problème! a = data.frame (a = gigue (1:10), b = gigue (1:10), c = gigue (1:10), d = gigue (1:10), e = gigue (1:10) , f = gigue (1:10), g = échantillon (gigue (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)) Les résultats diffèrent un peu et deviennent beaucoup plus similaires lorsque J'utilise des lambdas beaucoup plus élevés pour glmnet.
a11msp
Intrigant. Les coefficients semblent différer approximativement du facteur 10.
tomka

Réponses:

27

La différence que vous observez est due à la division supplémentaire par le nombre d'observations, N, que GLMNET utilise dans leur fonction objective et à la standardisation implicite de Y par son écart-type d'échantillon comme indiqué ci-dessous.

12NysyXβ22+λβ22/2

où nous utilisons au lieu de pour , 1 / ( n - 1 ) s y s y = i ( y i - ˉ y ) 21/n1/(n1)sy

sy=i(yiy¯)2n

En différenciant par rapport à la version bêta, en mettant l'équation à zéro,

XTXβXTysy+Nλβ=0

Et en résolvant pour la bêta, on obtient l'estimation,

β~gLMNET=(XTX+Nλjep)-1XTysy

Pour récupérer les estimations (et leurs pénalités correspondantes) sur la métrique d'origine de Y, GLMNET multiplie les estimations et les lambdas par et renvoie ces résultats à l'utilisateur,sy

X u n s t d . = s y λ

β^gLMNET=syβ~gLMNET=(XTX+Nλjep)-1XTy
λunst.=syλ

Comparez cette solution avec la dérivation standard de la régression de crête.

β^=(XTX+λjep)-1XTy

Notez que est mis à l'échelle par un facteur supplémentaire de N. De plus, lorsque nous utilisons la fonction ou , la pénalité sera implicitement mise à l'échelle par . Autrement dit, lorsque nous utilisons ces fonctions pour obtenir les estimations de coefficient pour certains , nous obtenons effectivement des estimations pour .1 / s y λ λ = λ / s yλpredict()coef()1/syλλ=λ/sy

Sur la base de ces observations, la peine utilisée dans GLMNET doit être mis à l' échelle par un facteur de .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

Les résultats se généralisent à l'inclusion d'une variable d'interception et de variables X standardisées. Nous modifions une matrice X standardisée pour inclure une colonne d'unités et la matrice diagonale pour avoir une entrée de zéro supplémentaire en position [1,1] (c'est-à-dire ne pas pénaliser l'interception). Vous pouvez ensuite standardiser les estimations par leurs écarts-types d'échantillon respectifs (assurez-vous à nouveau d'utiliser 1 / n lors du calcul de l'écart-type).

β^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 ajouté pour montrer X normalisé sans interception:

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
la source
3
+6. Bienvenue sur CV et merci d'avoir répondu à cette vieille question d'une manière aussi claire.
amoeba dit Reinstate Monica
1
Ce devrait être la matrice d'identité au lieu de dans la solution de , correct? ˜ βββ~
user1769197
Je remarque également que pour la deuxième partie où vous avez dit "Les résultats se généralisent à l'inclusion d'une variable d'interception et de variables X normalisées"; pour cette partie, si vous excluez l'interception, puis en suivant les mêmes calculs, les résultats de glmnet deviennent différents du calcul manuel.
user1769197
Correct, j'ai mis à jour la solution avec la matrice d'identité à la place de au besoin. J'ai vérifié la solution pour X normalisé sans interception et j'obtiens toujours des résultats identiques (voir le code supplémentaire ci-dessus). β
skijunkie
3

Selon https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html , lorsque la famille est gaussian, glmnet()devrait minimiser

(1)12nje=1n(yje-β0-XjeTβ)2+λj=1p(α|βj|+(1-α)βj2/2).

Lors de l'utilisation glmnet(x, y, alpha=1)pour ajuster le lasso avec les colonnes en normalisées, la solution pour la pénalité rapportée est la solution pour minimiser Cependant, au moins dans , lors de l'utilisation pour ajuster la régression de crête, la solution pour une pénalité rapportée est la solution pour minimiser où est l'écart-type de . Ici, la pénalité aurait dû être signalée comme .Xλ

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

Ce qui pourrait arriver, c'est que la fonction standardise d'abord en puis minimise qui est effectivement de minimiser ou de manière équivalente, pour minimiser yy0

(2)12nje=1n(y0je-XjeTγ)2+ηj=1p(α|γj|+(1-α)γj2/2),
12nsy2je=1n(yje-β0-XjeTβ)2+ηαsyj=1p|βj|+η1-α2sy2j=1pβj2,
12nje=1n(yje-β0-XjeTβ)2+ηsyαj=1p|βj|+η(1-α)j=1pβj2/2.

Pour le lasso ( ), redimensionner pour signaler la pénalité car est logique. Ensuite, pour tous les , doit être signalé comme pénalité pour maintenir la continuité des résultats sur . C'est probablement la cause du problème ci-dessus. Cela est dû en partie à l'utilisation de (2) pour résoudre (1). Ce n'est que lorsque ou qu'il existe une certaine équivalence entre les problèmes (1) et (2) (c'est-à-dire une correspondance entre le dans (1) et le dans (2)). Pour tout autreα=1ηηsyαηsyαα=0α=1ληα(0,1), les problèmes (1) et (2) sont deux problèmes d'optimisation différents, et il n'y a pas de correspondance biunivoque entre le dans (1) et le dans (2).λη

Chun Li
la source
1
Je ne vois pas en quoi votre réponse diffère de la précédente. Pourriez-vous expliquer, s'il vous plaît?
Firebug
1
@Firebug Je voulais expliquer pourquoi la fonction signale le lambda de cette façon, qui ne semble pas naturel lorsqu'elle est vue uniquement du point de vue de la régression des crêtes, mais qui a du sens (ou doit être de cette façon) lorsqu'elle est vue du point de vue de l'ensemble du spectre y compris la crête et le lasso.
Chun Li