Логистическая регрессия: Scikit Learn против glmnet

15

Я пытаюсь продублировать результаты из sklearnбиблиотеки логистической регрессии, используя glmnetпакет в R.

Из sklearnрегрессионной логистической документации , она пытается свести к минимуму функцию стоимости при l2 казни

min_{w, c} \frac{1}{2} w^{T} w + C \sum_{i = 1}^{N} \log (\exp (- y_{i} (X_{i}^{T} w + c)) + 1)

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

Из виньеток из glmnetего реализация минимизирует несколько иной стоимость функции

min_{β, β_{0}} - [\frac{1}{N} Σ_{я знак равно 1}^{N} Y_{я} (β_{0} + {Икс}_{я}^{T} β) - журнал (1 + е^{(β_{0} + {Икс}_{я}^{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]$

С некоторыми изменениями во втором уравнении и установкой $\alpha=0$ ,

λ \underset{β, β_{0}}{мин} \frac{1}{N λ} Σ_{я знак равно 1}^{N} [- Y_{я} (β_{0} + {Икс}_{я}^{T} β) + журнал (1 + е^{(β_{0} + {Икс}_{я}^{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$

которая отличается от sklearnфункции стоимости только на коэффициент $\lambda$ если установлено $\frac1{N\lambda}=C$ , поэтому я ожидал одинаковую оценку коэффициента от двух пакетов. Но они разные. Я использую набор данных из учебника UCLA idre , прогнозирование admitна основе gre, gpaи rank. Есть 400 наблюдений, так и с $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

Результат Rкак-то близок к логистической регрессии без регуляризации, как можно увидеть здесь . Я что-то упускаю или делаю что-то явно не так?

Обновление: я также попытался использовать LiblineaRпакет Rдля выполнения того же процесса, но получил еще один другой набор оценок ( liblinearтакже решающий в 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

Обновление 2: отключение стандартизации в glmnetдает:

> 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

r logistic python scikit-learn glmnet hurrikale
источник

Вы когда-нибудь понимали это?

Хьюи,

8

$L_2$

Тем более что ваш greтермин в таком большем масштабе, чем другие переменные, это изменит относительные затраты на использование различных переменных для весов.

Также обратите внимание, что, включив в функции явный термин «перехват», вы упорядочиваете перехват модели. Обычно этого не делается, поскольку это означает, что ваша модель больше не является ковариантной для смещения всех меток на константу.

Дугал
источник

glmnetпозволяет отключить стандартизацию входных данных, но расчетные коэффициенты еще более различны, см. выше. Кроме того, я явно включил термин перехват в sklearnпотому что glmnetвключает один автоматически, так что это для того, чтобы убедиться, что вход для обеих моделей одинаковы.

ураган

2

@hurrikale Я думаю, что glmnet, вероятно, не регулирует перехват, но sklearn есть. Удалите столбец перехвата из Xи перейдите fit_intercept=True(по умолчанию) в LogisticRegression. Хотя, вероятно, что-то еще происходит.

Дугал

Я попробовал то, что вы предложили, и все же получил разные наборы коэффициентов: [-1.873, -0.0606, -1.175, -1.378, 0.00182, 0.2435]для sklearnи [-2.8181, 0.0001269, -0.0002259, -0.00020288, 0.00344, 0.000188]для glmnetв порядке [Intercept, rank_2, rank_3, rank_4, gre, gpa]. Меня беспокоит то, что они различаются как по величине, так и по положительному / отрицательному влиянию на вероятность, поэтому, не зная, почему они различаются, трудно выбрать один для интерпретации. И если в какой-либо реализации есть какая-либо ошибка, особенно важно, чтобы я знал, на какую из них можно положиться.

ураган

7

Ответ Дугала правильный, вы пересекаете перехват в, sklearnно не в R. Убедитесь, что вы используете, solver='newton-cg'так как по умолчанию solver ( 'liblinear') всегда пересекает перехват.

ср https://github.com/scikit-learn/scikit-learn/issues/6595

TomDLT
источник

Установка solver='newton-cg'сделала результаты sklearnи statsmodelsсоответствует. Большое спасибо.

Ирэн

0

Вы должны также использовать L1_wt=0аргумент вместе с alphaв fit_regularized()вызове.

Этот код в 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)

эквивалентно следующему коду из sklearn:

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

Надеюсь, это поможет!

Прафул Гупта
источник

Логистическая регрессия: Scikit Learn против glmnet

Ответы: