Python中的卡方检验

Question

我在R中使用了以下代码，以确定观察值（例如20、20、0和0）与预期值/比率（例如每种情况的25%）的匹配程度：

> chisq.test(c(20,20,0,0), p=c(0.25, 0.25, 0.25, 0.25))

    Chi-squared test for given probabilities

data:  c(20, 20, 0, 0)

X-squared = 40, df = 3, p-value = 1.066e-08

我该如何在Python中复制这个过程？我尝试使用scipy的chisquare函数，但我得到的结果非常不同；我不确定这是否是正确的函数。我已经搜索了scipy文档，但由于文档超过1000页，所以感觉很令人生畏；而numpy文档甚至比它多50％。

- SabreWolfy

Answer 1

scipy.stats.chisquare函数需要输入实际观测到的频数和期望的绝对频数，而不是比率。你可以通过以下方式获得所需的输入：

>>> observed = np.array([20., 20., 0., 0.])
>>> expected = np.array([.25, .25, .25, .25]) * np.sum(observed)
>>> chisquare(observed, expected)
(40.0, 1.065509033425585e-08)

尽管在期望值均匀分布于类别的情况下，你可以省略对期望值的计算：

>>> chisquare(observed)
(40.0, 1.065509033425585e-08)

第一个返回值是 χ² 统计量，第二个是测试的 p 值。

- Fred Foo

---

谢谢 - 这正是我需要的 :) 我确实尝试将观察值转换为比率，但可能误读了文档，因为该函数需要频率而不是比率。一些预期频率将不相等，所以我想我应该选择一个不同的例子 :) - SabreWolfy

Answer 2

我想指出的是，尽管答案在语法上似乎是正确的，但你不应该在你的例子中使用卡方分布，因为你观察到的频率太小，无法进行准确的卡方检验。

当每个类别中的观察值或期望值过小时，此测试无效。一个典型的规则是所有的观察值和期望值都应至少为5。请参阅： http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html#scipy.stats.chisquare

- emaxwell

---

据我所知，规则仅基于预期频率，而不是观察到的频率，因此此示例（其中所有预期频率均相等于10）应该是可以的。如果预期频率太小，R会发出警告...例如，http://fds.oup.com/www.oup.com/pdf/13/9780199219995.pdf; http://www.stat.sfu.ca/~cschwarz/Stat-650/Notes%/PDFbigbook-JMP/JMP-part016.pdf; http://udel.edu/~mcdonald/statsmall.html（谷歌搜索“卡方期望”规则”）。我不会因为您正确引用了（错误的？？？）Scipy文档而对您进行投票。 - Ben Bolker

谢谢你指出来！我可能要补充一下，对于80％的类别，只需要预期频率达到至少5即可。 - emaxwell

这取决于您希望近似值有多精确。根据我的经验，“所有预期频率>=5”的经验法则更常被引用。您所引用的规则稍微宽松一些（出于好奇，可以提供引文/链接吗？）。 - Ben Bolker

在您提供的第一个链接中（http://fds.oup.com/www.oup.com/pdf/13/9780199219995.pdf），它说要允许五分之一，同时不允许期望频率为0。 - emaxwell

Answer 3

另一种方法是从Python中调用您的R代码。您可以通过以下方式实现：

• 将R脚本作为命令行工具运行。请参见this link，了解有关使用Rscript从命令行运行R脚本的更多信息。然后，您可以通过使用subprocess或os.system执行系统调用来从Python中运行R脚本。所有数据交换都通过文本或二进制文件完成。我喜欢这种方法，因为它非常简单，并且很容易从Python代码中单独调试R脚本。缺点是所有数据都通过硬盘传输，这可能会非常慢。
• 通过使用rpy或rpy2直接从Python中运行R代码。这样，集成更加紧密，但是此链接也引入了自己的小问题。例如，根据我的经验，通过rpy调用的R代码的调试要困难一些。

- Paul Hiemstra

---

谢谢您的建议。我之前使用过 rpy，但在这里决定不使用它，因为我可能需要传输相当大量的具有复杂结构的数据。 - SabreWolfy

只是想将这个选项添加到答案范围中，也许其他人会喜欢这种方法。 - Paul Hiemstra