Scikit-learn χ²（卡方）统计量和相应的列联表

Posted 2023-02-23

【中文标题】Scikit-learn χ²（卡方）统计量和相应的列联表

【英文标题】：Scikit-learn χ² (chi-squared) statistic and corresponding contingency table

【发布时间】：2014-02-12 09:59:41

【问题描述】：

在 scikit-learn http://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.chi2.html 的卡方单变量特征选择函数的文档中，它指出

该分数可用于从 X 中选择具有最高 χ²（卡方）统计值的 n_features 特征，该统计必须包含布尔值或频率（例如，文档分类中的术语计数)，相对于类。

我很难理解相应的列联表是什么样的，尤其是在频率特征的情况下。

例如，考虑以下具有布尔特征和目标的数据集：

import numpy as np

>>> X = np.random.randint(2, size=50).reshape(10, 5)
array([[1, 0, 0, 0, 1],
       [1, 1, 0, 1, 1],
       [1, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 1],
       [1, 0, 0, 0, 1],
       [1, 0, 1, 1, 1],
       [0, 1, 1, 0, 0],
       [1, 0, 1, 1, 1],
       [1, 1, 1, 1, 0]])

>>> y = np.random.randint(2, size=10)
array([1, 0, 0, 0, 1, 1, 1, 1, 0, 1])

要构造关于第一个特征的列联表，我们可以这样做（请原谅我的 PEP8 违规）

import scipy as sp

>>> contingency_table = sp.sparse.coo_matrix(
...    (np.ones_like(y), (X[:, 0], y)), 
...    shape=(np.unique(X[:, 0]).shape[0], np.unique(y).shape[0])).A
array([[1, 2],
       [3, 4]])

所以现在我可以计算卡方统计量及其 p 值

>>> sp.stats.chi2_contingency(contingency_table)
(0.17857142857142855,
 0.67260381744151676,
 1,
 array([[ 1.2,  1.8],
       [ 2.8,  4.2]]))

这应该与 scikit-learn 的 chi2 一致

from sklearn.feature_selection import chi2

>>> chi2_, pval = chi2(X, y)
>>> chi2_[0], pval[0]
(0.023809523809523787, 0.87737055606414338)

...不。我误解了什么吗？

此外，在频率的情况下，列联表是什么样的？我以为会是这样的

contingency_table = sp.sparse.coo_matrix(
    (np.ones_like(y), (X[:, 0], y)), 
    shape=(X[:, 0].max()+1, np.unique(y).shape[0])).A

但对应的预期频率表很可能有几个零元素。

编辑：

为了进一步澄清，请考虑第一个特征X[:, 0]，即gender和目标y，例如handedness。

由此我们得到交叉表

                Right-handed    Left-handed (!right-handed)
Male            1               2
Female (!male)  3               4

我们可以通过设置期望频率，使用卡方检验来评估两个比例之间差异的显着性

sklearn.feature_selection.chi2 直接执行此操作，无需显式计算表格，并使用等效于scipy.stats.chisquare 的更有效过程获得分数。

在明确列举了上面显示的表格后，我想在应用scipy.stats.chi2_contingency 时验证它是否与chi2 一致，但令我沮丧的是，事实并非如此。我想问一下为什么不是。

【问题讨论】：

有趣：最近一小时我一直在尝试解决一个非常相似的问题，但方向相反：从本教程youtube.com/watch?v=VskmMgXmkMQ中找到的第一个示例列联表中，验证对应的提供给sklearn.feature_selection.chi2 的 2d 布尔数据集（X=is_male, y=is_deposit）应该会产生相同的结果，并尝试通过研究源代码来理解为什么不是这样。

我发现您从陈述一个包含 5 个特征的数据集开始，然后仅使用第一个特征构建您的案例的其余部分，这让我感到困惑。我认为这就是让@larsmans 的答案在这种情况下令人困惑的原因，因为他正在回答一个不同的问题。

这很酷。看起来权变矩阵是一个 2x2 矩阵，并且您在生成目标时获得所有特征集的单个 p-val，而 feature_selection chi2 在生成目标时独立于其他特征显示每个特征的 p-val

【参考方案1】：

鉴于您的数据，

>>> X = array([[1, 0, 0, 0, 1],
...        [1, 1, 0, 1, 1],
...        [1, 0, 0, 0, 0],
...        [0, 0, 0, 0, 0],
...        [0, 0, 0, 0, 1],
...        [1, 0, 0, 0, 1],
...        [1, 0, 1, 1, 1],
...        [0, 1, 1, 0, 0],
...        [1, 0, 1, 1, 1],
...        [1, 1, 1, 1, 0]])
>>> y = array([1, 0, 0, 0, 1, 1, 1, 1, 0, 1])

这是feature_selection.chi2 计算的：

>>> Y = np.vstack([1 - y, y])
>>> observed = np.dot(Y, X)
>>> observed
array([[3, 1, 1, 2, 2],
       [4, 2, 3, 2, 4]])

这些是每个类观察到的特征频率，即列联表。那么期望值：

>>> feature_count = X.sum(axis=0)
>>> class_prob = Y.mean(axis=1)
>>> expected = np.dot(feature_count.reshape(-1, 1), class_prob.reshape(1, -1)).T
>>> expected
array([[ 2.8,  1.2,  1.6,  1.6,  2.4],
       [ 4.2,  1.8,  2.4,  2.4,  3.6]])

最后，它运行一个 χ² 检验：

>>> from scipy.stats import chisquare
>>> score, pval = chisquare(observed, expected)
>>> score
array([ 0.02380952,  0.05555556,  0.375     ,  0.16666667,  0.11111111])
>>> pval
array([ 0.87737056,  0.81366372,  0.54029137,  0.6830914 ,  0.73888268])

分数是相关的部分：它们用于按判别力对特征进行排序。请注意，您将获得一个分数和一个 p 值每个特征。

【讨论】：

感谢您的回答。我认为我在提出这个问题方面做得并不好，因此，我们可能不在同一个页面上。从源代码（一种更有效但可读性较差的变体）和文档中可以清楚地看到您所说的。我所追求的是“R x C 表” s.t. chi2_contingency 将产生与 sklearn.feature_selection.chi2 相同的结果。

您展示的示例可以解释为将此类表分解为要传递给scipy.stats.chisquare 的数组，这将产生相同的结果。但是，我有兴趣获得原始的 R x C 表。

@larsmans：对于标准列联表分析，您的 expected 数组不正确。预期频率数组必须与观察频率数组具有相同的边际总和（即行和列总和）。

@WarrenWeckesser chisquare docstring 建议不然，你是说错了吗？

对不起，我应该写一个更长的评论（或者根本没有评论！）。您的expected 对sklearn.feature_selection.chi2 进行的计算是正确的。正如您所说，它没有对observed 列联表进行单一的 chi2 测试（当我说它对于“标准列联表分析”不正确时，我正在考虑这个测试）。它也没有对每列进行 2x2 列联表分析。我认为后一点是@tiao问题的本质。

【参考方案2】：

考虑X 的列x。 sklearn.feature_selection.chi2 测试是否 x 为 1 的 y 值的频率与 y 的频率一致 全体人口。 （@larsman 的回答显示了如何使用 numpy 和 scipy 重现计算。）这与标准的 2x2 列联表不同 分析x 和y。在 2x2 列联表分析中，y 的频率 其中x 为 0 也有助于测试。

假设我们为x 和y 形成列联表：

    | y=0  y=1
----+---------
x=0 |  a    b
x=1 |  c    d

设 n = a + b + c + d。这是样本数（即与 len(x) 和 len(y) 相同）。

令 nx = c + d。这是1 在x 中出现的次数。

令 py1 = (b + d)/n。这是 y 为 1 的全部人口的比例。

sklearn.feature_selection.chi2 使用预期的对 [c, d] 执行 chi2 测试 值 [(1-py1)*nx, py1*nx]。这与标准列联表不同 分析 2x2 表。

这是一个极端的例子。假设 x 和 y 的 2x2 列联表是

    |  y=0  y=1
----+----------
x=0 |   8    8
x=1 |  20  188

sklearn 计算得出的 chi2 分数为 1.58，p 值为 0.208。

scipy.stats.chi2_contingency 的列联表分析得出 chi2 得分为 18.6，p 值为 1.60e-5。

【讨论】：

这与 2x2 表的标准列联表分析不同。 这是我问题的本质；这不是标准的列联表分析吗？如果不是，它到底应该是什么（链接到文章/论文？）我认为这应该在文档中明确说明。