组合行和列来创建用于 Fisher 精确检验的 2x2 表

Question

我需要使用 python 对以下交叉表

ct

执行独立性测试：

由于有些值小于5，我无法执行独立性的卡方检验。相反，我需要执行费舍尔精确检验。

由于 Fisher 在 Scipy 上的精确测试实现仅支持 2x2 表，因此我实现了以下解决方案：

from scipy.stats import fisher_exact

# Combine rows and columns to create a 2x2 table
table_2x2 = np.array([[ct[1][4] + ct[2][4] + ct[1][3] + ct[2][3], ct[3][4] + ct[4][4] + ct[3][3] + ct[3][3]],
                      [ct[1][2] + ct[2][2] + ct[1][1] + ct[2][1], ct[3][2] + ct[4][2] + ct[3][1] + ct[4][1]]])

# Perform Fisher's exact test on the 2x2 table
odds_ratio, p_value = fisher_exact(table_2x2)

# Display the results
print(f'Odds Ratio: {odds_ratio}')
print(f'P-value: {p_value}')

您认为这是一个有效的解决方案吗？如果没有，还有其他建议在Python中实现这个吗？ R 解决方案不是一个选项，因为这应该在 Python 中完成。

提前谢谢您。

Answer 1

如果没有，还有其他建议可以用Python实现这个吗？

如果您愿意接受随机排列测试，您可以使用

scipy.stats.permutation_test

创建自己的测试。我们将使用与

scipy.stats.chi2_contingency

相同的检验统计量，但原假设将类似于 Fisher 精确检验。

import numpy as np
from scipy import stats

table = np.asarray([[20, 49, 25, 4], [35, 54, 43, 12], [27, 44, 29, 8], [7, 20, 16, 4]])
ref = stats.chi2_contingency(table)

def untab(table):
    # convert 2d contingency table to two samples
    x = []
    y = []
    m, n = table.shape
    for i in range(m):
        for j in range(n):
            count = table[i, j]
            x += [i]*count
            y += [j]*count
    return np.asarray(x), np.asarray(y)

x, y = untab(table)

def statistic(x):
    table = stats.contingency.crosstab(x, y).count
    return stats.chi2_contingency(table).statistic

res = stats.permutation_test((x,), statistic, alternative='greater', 
                             permutation_type='pairings')

print(res.pvalue, ref.pvalue)  # 0.6592 0.6500840391351904

对于原始帖子中显示的列联表，与卡方检验相比，p 值几乎没有差异。尽管表中的某些计数很小，但零分布似乎与具有适当自由度数的卡方分布非常相似：

import matplotlib.pyplot as plt
plt.hist(res.null_distribution, bins=30, density=True, label='normalized histogram')

# see https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2_contingency.html
# for degrees of freedom
df = table.size - sum(table.shape) + table.ndim - 1
dist = stats.chi2(dof)
x = np.linspace(0, 40, 300)
plt.plot(x, dist.pdf(x), label='chi2')
plt.legend()

有关理论（和实践）的更多信息，请参阅有关 重采样和蒙特卡罗方法的 SciPy 教程，尤其是 2c、相关样本排列测试。