scikit learn：如何检查系数的显着性答案

【问题标题】：scikit learn: how to check coefficients significancescikit learn：如何检查系数的显着性
【发布时间】：2014-09-27 04:21:43
【问题描述】：

我尝试使用 SKLearn 对一个相当大的数据集进行 LR，该数据集有大约 600 个虚拟变量，只有很少的区间变量（我的数据集中有 300 K 行），由此产生的混淆矩阵看起来很可疑。我想检查返回的系数和方差分析的重要性，但我找不到如何访问它。有可能吗？对于包含大量虚拟变量的数据，最佳策略是什么？非常感谢！

【问题讨论】：

如果您的逻辑回归对象称为lr，请尝试查看lr.coef_。这是你要找的吗？
不，好吧，coef_ 是系数值，我想要的是这个值的意义：z-score 和 p-value。当您假设系数为 0 的检验假设（零假设 H_0=0）和备择假设 H_1!=0，然后 p 值基本上告诉您是否可以拒绝 H_0（当 H_0 很小时）或不是（当 H_0->1 时）
对于逻辑回归，我感觉你只能在每个样本的coef_ 上使用重采样和建立经验分布来获得那些。
嗯，是的，但我想知道是否有 sklearn 的内置方法，例如 R 中“glm 类”对象的摘要......
如果这有帮助，您还可以检查稳定性选择和随机逻辑回归的 sklearn 实现。这些可以为您提供稳定的功能选择。

标签： scikit-learn logistic-regression anova dummy-data

【解决方案1】：

Scikit-learn 故意不支持统计推断。如果您想要开箱即用的系数显着性检验（以及更多），您可以使用 Statsmodels 中的 Logit estimator。这个包模仿了R中的接口glm模型，所以你会觉得它很熟悉。

如果您仍想坚持使用 scikit-learn LogisticRegression，您可以使用渐近近似来分布最大似然估计。准确地说，对于最大似然估计向量theta，其方差-协方差矩阵可以估计为inverse(H)，其中H 是theta 处的对数似然的Hessian 矩阵。这正是下面的函数所做的：

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 摘要。

【讨论】：

@Rocketq 1) 您能否定义一下“可靠的 p 值”是什么意思？它是 MLE p 值的一个特例。所以我建议寻找“最大似然估计的渐近性质”的理论，以全面了解其可靠性。
@Rocketq 2) 是的，Statsmodels 确实以相同的方式计算逻辑回归的 p 值。参数的协方差矩阵（statsmodels.base.model.LikelihoodModelResults.normalized_cov_params 属性）在 statsmodels.base.model.LikelihoodModel.fit 方法中计算为逆 Hessian，并进一步用于 p 值估计和其他目的。据我所知，SPSS 基本上是一样的。
是的，这个p值正是Wald检验的意义。如果假设为真，两者都基于值(estimate-hypothesis) / std.dev(estimate) 是渐近标准正态的假设。见en.wikipedia.org/wiki/Wald_test#Test_on_a_single_parameter
@kand 因为专注。他们只是无法支持一切（资源非常有限），他们选择更全面地涵盖不同的 ML 算法，而不是做其他事情。
C 对系数本身有很大的影响（当然，通过它们对 p 值也有影响）。 Scikit-learn 默认使用 C=1； Statsmodels 根本不进行正则化（相当于 C=infinity）。因此，如果我们希望 scikit-learn 和 statsmodels 具有相似的系数，我们需要在 scikit-learn 中将 C 设置得很高。