Quais são as diferenças entre a regressão de Ridge usando o glmnet de R e o scikit-learn do Python?

11

Estou examinando a seção LAB §6.6 sobre Regressão de Ridge / Lasso no livro 'Uma Introdução à Aprendizagem Estatística com Aplicações em R' de James, Witten, Hastie, Tibshirani (2013).

Mais especificamente, estou tentando aplicar o Ridgemodelo scikit-learn ao conjunto de dados 'Hitters' do pacote R 'ISLR'. Eu criei o mesmo conjunto de recursos, como mostrado no código R. No entanto, não consigo me aproximar dos resultados do glmnet()modelo. Eu selecionei um parâmetro de ajuste L2 para comparar. (argumento 'alpha' no scikit-learn).

Pitão:

regr = Ridge(alpha=11498)
regr.fit(X, y)

http://nbviewer.ipython.org/github/JWarmenhoven/ISL-python/blob/master/Notebooks/Chapter%206.ipynb

R:

Observe que o argumento alpha=0in glmnet()significa que uma penalidade de L2 deve ser aplicada (regressão de Ridge). A documentação avisa para não inserir um único valor para lambda, mas o resultado é o mesmo que no ISL, onde um vetor é usado.

ridge.mod <- glmnet(x,y,alpha=0,lambda=11498)

O que causa as diferenças?

Edit:
Ao usar a penalized()partir do pacote penalizado em R, os coeficientes são os mesmos do scikit-learn.

ridge.mod2 <- penalized(y,x,lambda2=11498)

Talvez a pergunta também possa ser: 'Qual é a diferença entre glmnet()e penalized()ao fazer a regressão de Ridge?

Novo wrapper python para o código Fortran real usado no pacote R glmnet
https://github.com/civisanalytics/python-glmnet

Jordi
fonte
5
Totalmente familiarizado com a regressão da crista glmnet. Mas, por padrão, sklearn.linear_model.Ridgefaz estimativa de interceptação não-penalizada (padrão) e a penalidade é tal que ||Xb - y - intercept||^2 + alpha ||b||^2é minimizada b. Pode haver fatores 1/2ou 1/n_samplesambos na frente da penalidade, tornando os resultados diferentes imediatamente. Para fatorar o problema de escala da penalidade, defina a penalidade como 0 nos dois casos, resolva quaisquer discrepâncias e verifique o que a adição da penalidade faz. E entre IMHO aqui é o lugar certo para fazer esta pergunta.

Respostas:

9

Minha resposta está faltando um fator de . Consulte @visitors answer abaixo para obter a comparação correta.1N


Aqui estão duas referências que devem esclarecer o relacionamento.

A documentação do sklearn diz que linear_model.Ridgeotimiza a seguinte função objetivo

|Xβy|22+α|β|22

O artigo glmnet diz que a rede elástica otimiza a seguinte função objetivo

|Xβy|22+λ(12(1α)|β|22+α|β|1)

Observe que as duas implementações usam maneiras totalmente diferentes, o sklearn usa para o nível geral de regularização, enquanto o glmnet usa para esse fim, reservando para negociação entre a regularização do cume e do laço. ααλα

Comparando as fórmulas, parece que a configuração e no glmnet deve recuperar a solução .α=0λ=2αsklearnlinear_model.Ridge

Matthew Drury
fonte
E eu também senti muita falta disso no comentário de @eickenberg. Eu tenho que usar standardize = FALSEem glmnet()obter os resultados idênticos.
Jordi
@ Jordi Você definitivamente deve padronizar se estiver usando linear_model.Ridgepara qualquer análise do mundo real.
Matthew Drury
Entendo que o linear_model.Ridgemodelo sklearn padroniza os recursos automaticamente. A normalização é opcional. Gostaria de saber por que preciso desativar a padronização glmnet()para que os modelos produzam resultados idênticos.
Jordi
10

A resposta de Matthew Drury deve ter um fator de 1 / N. Mais precisamente...

A documentação da glmnet afirma que a rede elástica minimiza a função de perda

1NXβy22+λ(12(1α)β22+αβ1)

A documentação do sklearn diz que linear_model.Ridgeminimiza a função de perda

Xβy22+αβ22

o que equivale a minimizar

1NXβy22+αNβ22

Para obter a mesma solução do glmnet e do sklearn, as duas funções de perda devem ser iguais. Isso significa configurar e no glmnet.α=0λ=2Nαsklearn

library(glmnet)
X = matrix(c(1, 1, 2, 3, 4, 2, 6, 5, 2, 5, 5, 3), byrow = TRUE, ncol = 3)
y = c(1, 0, 0, 1)
reg = glmnet(X, y, alpha = 0, lambda = 2 / nrow(X))
coef(reg)

Saída glmnet: –0.03862100, –0.03997036, –0.07276511, 0.42727955

import numpy as np
from sklearn.linear_model import Ridge
X = np.array([[1, 1, 2], [3, 4, 2], [6, 5, 2], [5, 5, 3]])
y = np.array([1, 0, 0, 1])
reg = Ridge(alpha = 1, fit_intercept = True, normalize = True)
reg.fit(X, y)
np.hstack((reg.intercept_, reg.coef_))

Saída do sklearn: –0.03862178, –0.0399697, –0.07276535, 0.42727921

Visitante
fonte
4
As diferentes definições de parâmetros e seu dimensionamento usados ​​em diferentes bibliotecas são uma fonte comum de confusão.
AaronDefazio 30/03
1
Eu não esperaria que tanto Gung quanto eu entendêssemos errado.
Michael R. Chernick
2
Sim, vocês dois entenderam errado. Seus motivos para rejeitar minha edição deixam claro que vocês dois não viram meu comentário "Fator ausente de 1 / N" em stats.stackexchange.com/review/suggested-edits/139985
visitante
Sua edição provavelmente foi rejeitada porque mudou muito mais do que apenas o que você reivindica. Se você deseja editar minha postagem e alterar apenas o fator ausente, faça isso, mas alterar meus links, texto e código também é um exagero. Os comentários sobre o seu tratamento injusto na sua resposta não são adequados e, não relacionados ao conteúdo da pergunta, remova-os. Sua redação também plagiou minha resposta. Esse não é o caminho certo para responder a uma edição rejeitada. Adoraríamos suas valiosas contribuições à nossa comunidade, mas familiarize-se com nossas normas antes de nos eviscerar.
Matthew Drury
1
@visitor Desculpe se saí um pouco áspero. Eu realmente só deveria estar tentando comunicar que você parece ser um bom colaborador em potencial para o site e quero que você tenha uma boa experiência. Temos algumas normas sociais, como qualquer outro grupo, e você terá uma experiência melhor se ficar ciente delas. Ainda acho que "a resposta de Matthew Drury está errada" é bastante dura, certamente existem maneiras melhores de comunicar que minha resposta está erroneamente faltando um fator de . "A resposta de X está errada" é um ataque pessoal. 1N
Matthew Drury