Regresja logistyczna: Scikit Learn vs glmnet

15

Próbuję powielić wyniki z sklearnbiblioteki regresji logistycznej przy użyciu glmnetpakietu w języku R.

Z dokumentacjisklearn regresji logistycznej próbuje zminimalizować funkcję kosztu w ramach kary l2

min_{w, do} \frac{1}{2)} w^{T.} w + do \sum_{ja = 1}^{N.} \log (\exp (- y_{ja} (X_{ja}^{T.} w + do)) + 1)

$\min_{w,c} \frac12 w^Tw + C\sum_{i=1}^N \log(\exp(-y_i(X_i^Tw+c)) + 1)$

Z winiet z glmnetjego realizacja minimalizuje nieco inny kosztów funkcji

min_{β, β_{0}} - [\frac{1}{N.} \sum_{ja = 1}^{N.} y_{ja} (β_{0} + x_{ja}^{T.} β) - \log (1 + {mi}^{(β_{0} + x_{ja}^{T.} β)})] + λ [(α - 1) | | β | |_{2)}^{2)} / 2) + α | | β | |_{1}]

$\min_{\beta, \beta_0} -\left[\frac1N \sum_{i=1}^N y_i(\beta_0+x_i^T\beta)-\log(1+e^{(\beta_0+x_i^T\beta)})\right] + \lambda[(\alpha-1)||\beta||_2^2/2+\alpha||\beta||_1]$

Z w drugim równaniu i ustawiając $\alpha=0$ ,

λ min_{β, β_{0}} \frac{1}{N. λ} \sum_{ja = 1}^{N.} [- y_{ja} (β_{0} + x_{ja}^{T.} β) + \log (1 + {mi}^{(β_{0} + x_{ja}^{T.} β)})] + | | β | |_{2)}^{2)} / 2)

$\lambda\min_{\beta, \beta_0} \frac1{N\lambda} \sum_{i=1}^N \left[-y_i(\beta_0+x_i^T\beta)+\log(1+e^{(\beta_0+x_i^T\beta)})\right] + ||\beta||_2^2/2$

który różni się od sklearnfunkcji kosztu tylko współczynnikiem $\lambda$ jeśli jest ustawiony $\frac1{N\lambda}=C$ , więc spodziewałem się tego samego oszacowania współczynnika z dwóch pakietów. Ale są różne. Używam zestawu danych z UCLA Idre poradnik , przewidywania admitna podstawie gre, gpai rank. Istnieje 400 obserwacji, więc przy $C=1$ , $\lambda = 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

Dane Rwyjściowe są w jakiś sposób zbliżone do regresji logistycznej bez regularyzacji, co można zobaczyć tutaj . Czy coś brakuje lub robię coś ewidentnie niewłaściwego?

Aktualizacja: Próbowałem również użyć LiblineaRpakietu w Rcelu przeprowadzenia tego samego procesu, a mimo to otrzymałem inny zestaw szacunków ( liblinearjest to również solver 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

Aktualizacja 2: wyłączenie standaryzacji w glmnetdaje:

> 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
źródło

Czy kiedykolwiek to rozgryzłeś?

— Huey,

8

regresja logistyczna sklearn domyślnie nie standaryzuje danych wejściowych, co zmienia znaczenie regularyzacji ; prawdopodobnie robi to Glmnet. $L_2$

Zwłaszcza, że twój gretermin jest na tak większą skalę niż inne zmienne, zmieni to względne koszty użycia różnych zmiennych dla wag.

Pamiętaj również, że dołączając do funkcji wyraźny termin przechwytywania, regulujesz przechwytywanie modelu. Zasadniczo tak się nie dzieje, ponieważ oznacza to, że model nie jest już współzmiennym z przesunięciem wszystkich etykiet o stałą.

— Dougal
źródło

glmnetpozwala wyłączyć standaryzację wejść, ale szacowane współczynniki są jeszcze bardziej różne, patrz wyżej. Ponadto wyraźnie uwzględniłem pojęcie przechwytywania w, sklearnponieważ glmnetobejmuje ono jeden automatycznie, aby upewnić się, że dane wejściowe do obu modeli są takie same.

— hurrikale

2

@ururikale Myślę, że glmnet prawdopodobnie nie reguluje przechwytywania, ale sklearn jest. Upuść kolumnę przechwytującą z Xi przekaż fit_intercept=True(domyślnie) do LogisticRegression. Prawdopodobnie dzieje się też coś innego.

— Dougal

Próbowałem tego, co zasugerowałeś, a jednak uzyskałem różne zestawy współczynników: [-1.873, -0.0606, -1.175, -1.378, 0.00182, 0.2435]dla sklearni [-2.8181, 0.0001269, -0.0002259, -0.00020288, 0.00344, 0.000188]dla glmnetkolejności [Intercept, rank_2, rank_3, rank_4, gre, gpa]. Obawiam się, że różnią się one zarówno pod względem wielkości, jak i pozytywnie / negatywnie wpływają na prawdopodobieństwo, więc nie wiedząc, dlaczego się różnią, trudno wybrać taką, która będzie interpretowana. A jeśli jest jakakolwiek szansa na błąd w jednej z implementacji, szczególnie ważne jest, aby wiedzieć, na której polegać.

— hurrikale

7

Odpowiedź Dougala jest poprawna, regulujesz przechwytywanie w, sklearnale nie w R. Upewnij się, że używasz, solver='newton-cg'ponieważ default solver ( 'liblinear') zawsze reguluje przechwycenie.

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

— TomDLT
źródło

Ustawienie solver='newton-cg'wykonane z wyników sklearni statsmodelsspójne. Wielkie dzięki.

— Irene

0

Należy również użyć L1_wt=0argumentu oraz alphaw fit_regularized()połączeniu.

Ten kod w 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)

jest równoważne z następującym kodem z sklearn:

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

Mam nadzieję, że to pomoże!

— Praful Gupta
źródło