大于2x2列联表的Fisher精确检验

Question

你好，scipy库里有Fisher精确检验的实现，但它只适用于2x2的列联表。我想在大于2x2的表格上进行这个检验，比如5x2或5x3的表格。

我知道在R语言里有一个fisher.test的函数可以完成这个工作，但我想在我的Python代码中实现。

有没有人知道在Python中有没有可以处理更大表格的Fisher精确检验的实现？

另外，我也不确定在大于2x2的表格上做Fisher精确检验是否合适。

谢谢！

Answer 1

通过一些简单的计算，我们可以把2 x 2的超几何概率公式扩展到任何大小的r x c表格。例如，2 x 3的表格的概率可以这样计算：

(N1!N2!M1!M2!M3!) / (N!a!b!c!d!e!f!)

这里，N1和N2是每一行的总和，M1、M2、M3是每一列的总和，N是所有数据的总和，而a、b、c、d、e、f则是每个单元格里的数字。

你可以通过深度优先搜索（dfs）来实现费舍尔精确检验。代码如下：

import math

def _dfs(mat, pos, r_sum, c_sum, p_0, p):

    (xx, yy) = pos
    (r, c) = (len(r_sum), len(c_sum))

    mat_new = []

    for i in range(len(mat)):
        temp = []
        for j in range(len(mat[0])):
            temp.append(mat[i][j])
        mat_new.append(temp)

    if xx == -1 and yy == -1:
        for i in range(r-1):
            temp = r_sum[i]
            for j in range(c-1):
                temp -= mat_new[i][j]
            mat_new[i][c-1] = temp
        for j in range(c-1):
            temp = c_sum[j]
            for i in range(r-1):
                temp -= mat_new[i][j]
            mat_new[r-1][j] = temp
        temp = r_sum[r-1]
        for j in range(c-1):
            temp -= mat_new[r-1][j]
        if temp <0:
            return
        mat_new[r-1][c-1] = temp

        p_1 = 1
        for x in r_sum:
            p_1 *= math.factorial(x)
        for y in c_sum:
            p_1 *= math.factorial(y)

        n = 0
        for x in r_sum:
            n += x
        p_1 /= math.factorial(n)

        for i in range(len(mat_new)):
            for j in range(len(mat_new[0])):
                p_1 /= math.factorial(mat_new[i][j])
        if p_1 <= p_0 + 0.00000001:
            #print(mat_new)
            #print(p_1)
            p[0] += p_1
    else:
        max_1 = r_sum[xx]
        max_2 = c_sum[yy]
        for j in range(c):
            max_1 -= mat_new[xx][j]
        for i in range(r):
            max_2 -= mat_new[i][yy]
        for k in range(min(max_1,max_2)+1):
            mat_new[xx][yy] = k
            if xx == r-2 and yy == c-2:
                pos_new = (-1, -1)
            elif xx == r-2:
                pos_new = (0, yy+1)
            else:
                pos_new = (xx+1, yy)
            _dfs(mat_new, pos_new, r_sum, c_sum, p_0, p)


def fisher_exact(table):

    row_sum = []
    col_sum = []

    for i in range(len(table)):
        temp = 0
        for j in range(len(table[0])):
            temp += table[i][j]
        row_sum.append(temp)
    
    for j in range(len(table[0])):
        temp = 0
        for i in range(len(table)):
            temp += table[i][j]
        col_sum.append(temp)

    mat = [[0] * len(col_sum)] * len(row_sum)
    pos = (0, 0)

    p_0 = 1

    for x in row_sum:
        p_0 *= math.factorial(x)
    for y in col_sum:
        p_0 *= math.factorial(y)

    n = 0
    for x in row_sum:
        n += x
    p_0 /= math.factorial(n)

    for i in range(len(table)):
        for j in range(len(table[0])):
            p_0 /= math.factorial(table[i][j])

    p = [0]
    _dfs(mat, pos, row_sum, col_sum, p_0, p)

    return p[0]

你可以测试这段代码，比如：

print(fisher_exact([[1,24],[5,20],[14,11],[11,14]]))

这样会得到结果：

0.0001228337404686859

这个结果和R语言给出的结果是一样的。虽然可能有更优雅的方法，但这段代码确实能给出正确的结果。

Answer 2

是的，进行大于2x2的表格的Fisher精确检验是可以的。

目前在Python中没有特别干净、经过广泛测试的解决方案。一个解决办法是使用rpy2，然后从Python调用R的函数：

import numpy as np
import rpy2.robjects.numpy2ri
from rpy2.robjects.packages import importr
rpy2.robjects.numpy2ri.activate()

stats = importr('stats')
m = np.array([[4,4],[4,5],[10,6]])
res = stats.fisher_test(m)
print 'p-value: {}'.format(res[0][0])
>> p-value: 0.668165917041

另一个解决办法是深入研究R实现所使用的C代码，并直接调用那段代码。这里有一个链接，指向某人的GitHub项目，他们回到了最初的Fortran实现，并从Python调用了它。