я попытался выполнить LR с помощью SKLearn для довольно большого набора данных с ~ 600 фиктивными и всего несколькими интервальными переменными (и 300 тыс. строк в моем наборе данных), и полученная матрица путаницы выглядит подозрительно. Я хотел проверить значимость возвращаемых коэффициентов и ANOVA, но не могу найти, как получить к ним доступ. Это вообще возможно? И какова наилучшая стратегия для данных, содержащих множество фиктивных переменных? Большое спасибо!
scikit Learn: как проверить значимость коэффициентов
Ответы (1)
Scikit-learn намеренно не поддерживает статистические выводы. Если вам нужны готовые тесты значимости коэффициентов (и многое другое), вы можете использовать Logit оценщик из Statsmodels. Этот пакет имитирует модели интерфейса glm
в R, поэтому он может показаться вам знакомым.
Если вы все еще хотите придерживаться логистической регрессии scikit-learn, вы можете использовать асимтотическое приближение к распределению оценок максимального правдоподобия. А именно, для вектора оценок максимального правдоподобия theta
его матрица дисперсии-ковариации может быть оценена как inverse(H)
, где H
— матрица Гессе логарифмического правдоподобия при theta
. Это именно то, что делает функция ниже:
import numpy as np
from scipy.stats import norm
from sklearn.linear_model import LogisticRegression
def logit_pvalue(model, x):
""" Calculate z-scores for scikit-learn LogisticRegression.
parameters:
model: fitted sklearn.linear_model.LogisticRegression with intercept and large C
x: matrix on which the model was fit
This function uses asymtptics for maximum likelihood estimates.
"""
p = model.predict_proba(x)
n = len(p)
m = len(model.coef_[0]) + 1
coefs = np.concatenate([model.intercept_, model.coef_[0]])
x_full = np.matrix(np.insert(np.array(x), 0, 1, axis = 1))
ans = np.zeros((m, m))
for i in range(n):
ans = ans + np.dot(np.transpose(x_full[i, :]), x_full[i, :]) * p[i,1] * p[i, 0]
vcov = np.linalg.inv(np.matrix(ans))
se = np.sqrt(np.diag(vcov))
t = coefs/se
p = (1 - norm.cdf(abs(t))) * 2
return p
# test p-values
x = np.arange(10)[:, np.newaxis]
y = np.array([0,0,0,1,0,0,1,1,1,1])
model = LogisticRegression(C=1e30).fit(x, y)
print(logit_pvalue(model, x))
# compare with statsmodels
import statsmodels.api as sm
sm_model = sm.Logit(y, sm.add_constant(x)).fit(disp=0)
print(sm_model.pvalues)
sm_model.summary()
Выходы print()
идентичны, и они являются p-значениями коэффициентов.
[ 0.11413093 0.08779978]
[ 0.11413093 0.08779979]
sm_model.summary()
также печатает хорошо отформатированную сводку в формате HTML.
(estimate-hypothesis) / std.dev(estimate)
является асимтотически стандартной нормой, если гипотеза верна. См. en.wikipedia.org/wiki/Wald_test#Test_on_a_single_parameter.
- person David Dale; 11.01.2018
lr
, попробуйте взглянуть наlr.coef_
. Это то, что вы ищите? - person eickenberg   schedule 04.08.2014coef_
каждой выборки. - person eickenberg   schedule 05.08.2014n_samples <= n_features
. - person eickenberg   schedule 05.08.2014