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 glmnet
pour obtenir le même résultat qu'avec la crête?
r
ridge-regression
glmnet
John
la source
la source
Réponses:
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.
où nous utilisons au lieu de pour , 1 / ( n - 1 ) s y s y = ∑ i ( y i - ˉ y ) 21/n 1/(n−1) sy
En différenciant par rapport à la version bêta, en mettant l'équation à zéro,
Et en résolvant pour la bêta, on obtient l'estimation,
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 λ
Comparez cette solution avec la dérivation standard de la régression de crête.
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λ 1 / sy λ∗ λ = λ∗/ sy
predict()
coef()
Sur la base de ces observations, la peine utilisée dans GLMNET doit être mis à l' échelle par un facteur de .sy/ N
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).
Code ajouté pour montrer X normalisé sans interception:
la source
Selon https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html , lorsque la famille est
gaussian
,glmnet()
devrait minimiserLors de l'utilisationX λ
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 .glmnet_2.0-13
glmnet(x, y, alpha=0)
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 minimisery y0
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).λ η
la source