Regressão logística: Scikit Learn vs glmnet

15

Estou tentando duplicar os resultados da sklearnbiblioteca de regressão logística usando o glmnetpacote em R.

A partir da documentação dasklearn regressão logística , ele está tentando minimizar a função de custo sob pena de l2 min w , c 1

minW,c12WTW+CEu=1Nregistro(exp(-yEu(XEuTW+c))+1)

Das vinhetas de glmnet, sua implementação minimiza uma função de custo ligeiramente diferente

minβ,β0 0-[1NEu=1NyEu(β0 0+xEuTβ)-registro(1+e(β0 0+xEuTβ))]+λ[(α-1)||β||22/2+α||β||1]

Com alguns ajustes na segunda equação, e definindo , λ min β , β 0 1α=0 0

λminβ,β0 01NλEu=1N[-yEu(β0 0+xEuTβ)+registro(1+e(β0 0+xEuTβ))]+||β||22/2

sklearnλ1Nλ=CadmitgregparankC=1λ=0,0025

#python sklearn
df = pd.read_csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
y, X = dmatrices('admit ~ gre + gpa + C(rank)', df, return_type = 'dataframe')
X.head()
>  Intercept  C(rank)[T.2]  C(rank)[T.3]  C(rank)[T.4]  gre   gpa
0          1             0             1             0  380  3.61
1          1             0             1             0  660  3.67
2          1             0             0             0  800  4.00
3          1             0             0             1  640  3.19
4          1             0             0             1  520  2.93

model = LogisticRegression(fit_intercept = False, C = 1)
mdl = model.fit(X, y)
model.coef_
> array([[-1.35417783, -0.71628751, -1.26038726, -1.49762706,  0.00169198,
     0.13992661]]) 
# corresponding to predictors [Intercept, rank_2, rank_3, rank_4, gre, gpa]


> # R glmnet
> df = fread("https://stats.idre.ucla.edu/stat/data/binary.csv")
> X = as.matrix(model.matrix(admit~gre+gpa+as.factor(rank), data=df))[,2:6]
> y = df[, admit]
> mylogit <- glmnet(X, y, family = "binomial", alpha = 0)
> coef(mylogit, s = 0.0025)
6 x 1 sparse Matrix of class "dgCMatrix"
                    1
(Intercept)      -3.984226893
gre               0.002216795
gpa               0.772048342
as.factor(rank)2 -0.530731081
as.factor(rank)3 -1.164306231
as.factor(rank)4 -1.354160642

A Rsaída está de alguma forma próxima da regressão logística sem regularização, como pode ser visto aqui . Estou perdendo alguma coisa ou fazendo algo obviamente errado?

Atualização: Eu também tentei usar o LiblineaRpacote Rpara conduzir o mesmo processo, e ainda recebi outro conjunto de estimativas diferente ( liblineartambém é o solucionador sklearn):

> fit = LiblineaR(X, y, type = 0, cost = 1)
> print(fit)
$TypeDetail
[1] "L2-regularized logistic regression primal (L2R_LR)"
$Type
[1] 0
$W
            gre          gpa as.factor(rank)2 as.factor(rank)3 as.factor(rank)4         Bias
[1,] 0.00113215 7.321421e-06     5.354841e-07     1.353818e-06      9.59564e-07 2.395513e-06

Atualização 2: desativar a padronização no glmnetfornece:

> mylogit <- glmnet(X, y, family = "binomial", alpha = 0, standardize = F)
> coef(mylogit, s = 0.0025)
6 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept)      -2.8180677693
gre               0.0034434192
gpa               0.0001882333
as.factor(rank)2  0.0001268816
as.factor(rank)3 -0.0002259491
as.factor(rank)4 -0.0002028832
hurrikale
fonte
Você já descobriu isso?
Huey

Respostas:

8

eu2

Especialmente porque o seu gretermo está em uma escala tão maior que as outras variáveis, isso alterará os custos relativos do uso das diferentes variáveis ​​para pesos.

Observe também que, ao incluir um termo de interceptação explícito nos recursos, você está regularizando a interceptação do modelo. Isso geralmente não é feito, pois significa que seu modelo não é mais covariável para mudar todos os rótulos por uma constante.

Dougal
fonte
glmnetpermite desativar a padronização das entradas, mas os coeficientes estimados são ainda mais diferentes, veja acima. Além disso, incluí explicitamente o termo de interceptação, sklearnporque glmnetinclui um automaticamente, para garantir que as entradas para os dois modelos sejam as mesmas.
hurrikale
2
@hurrikale Acho que o glmnet provavelmente não está regularizando a interceptação, mas o sklearn está. Solte a coluna de interceptação Xe passe fit_intercept=True(o padrão) para LogisticRegression. Provavelmente há algo mais acontecendo também.
Dougal 27/03
Eu tentei o que você sugeriu e ainda tenho diferentes conjuntos de coeficientes: [-1.873, -0.0606, -1.175, -1.378, 0.00182, 0.2435]para sklearne [-2.8181, 0.0001269, -0.0002259, -0.00020288, 0.00344, 0.000188]para glmnetem ordem de [Intercept, rank_2, rank_3, rank_4, gre, gpa]. Minha preocupação é que eles diferem tanto em magnitude quanto em afetar positivamente / negativamente a probabilidade; portanto, sem saber por que eles diferem, é difícil escolher um para interpretar. E se houver algum erro em uma das implementações, é particularmente importante que eu saiba em qual delas confiar.
hurrikale
7

A resposta de Dougal está correta, você regulariza a interceptação em sklearnmas não em R. Certifique-se de usar, solver='newton-cg'pois o solver padrão ( 'liblinear') sempre regulariza a interceptação.

cf https://github.com/scikit-learn/scikit-learn/issues/6595

TomDLT
fonte
Cenário solver='newton-cg'fez os resultados sklearne statsmodelsconsistente. Muito obrigado.
irene
0

Você também deve usar o L1_wt=0argumento junto com alphana fit_regularized()chamada.

Este código em statsmodels:

import statsmodels.api as sm
res = sm.GLM(y, X, family=sm.families.Binomial()).fit_regularized(alpha=1/(y.shape[0]*C), L1_wt=0)

é equivalente ao seguinte código de sklearn:

from sklearn import linear_model
clf = linear_model.LogisticRegression(C = C)
clf.fit(X, y)

Espero que ajude!

Praful Gupta
fonte