はじめに
フィッシャーの正確検定 (Fisher's exact test) は、2つの集合 $A_1,A_2\subseteq U$ があったとき、それらの重複 $A_1\cap A_2$ が有意に大きい (または小さい) ことを示したいときに使われます。$U$ は $A_1$, $A_2$ が動ける範囲を表す全体集合です。以下では上側検定のみ考えます。
例えば、「袋の中に100個の玉が入っていて、まず10個の玉を取り出し、それらに印をつけてから袋に戻してよく混ぜ、次に20個の玉を取り出したとき、印が付いていた玉が4個あります。これは偶然に起こり得るといえるでしょうか？」のような問題設定です。100個というのが全体集合 $U$、最初に取り出した10個が集合 $A_1$、次に取り出した20個が集合 $A_2$、そのうち印が付いていた4個が $A_1\cap A_2$ です。
$|U|=N$, $|A_1|=n_1$, $A_2=n_2$, $|A\cap B|=x$ とおくと、$A_1$ と $A_2$ が独立の場合に $x$ は超幾何分布 (hypergeometric distribution) に従うので、p値はPythonの場合以下で計算されます。
from scipy import stats
stats.hypergeom.sf(x-1,N,n1,n2)
ここで sf は生存関数 (survival function) を表します。等号を含まないため、$x$ から1を引いています。
しかし、集合が3つ以上の場合や、$A_1$ と $A_2$ の動ける範囲が異なる場合 ($A_1\subseteq U_1$, $A_2\subseteq U_2$, $U_1\neq U_2$) への拡張は容易ではありません。そこで、フィッシャーの正確検定の代わりに二項検定を使う方法を考え、数値計算で比較してみました。
二項検定を使う方法
$A_1$, $A_2$ の大きさを固定する代わりに、集合 $U$ の各要素が $A_1$ に含まれる確率 $p_1=n_1/N$ と $A_2$ に含まれる確率 $p_2=n_2/N$ を固定します。期待値について、$E[|A_1|]=n_1$, $E[|A_2|]=n_2$ が成り立ちます。
$A_1$ と $A_2$ が独立の場合、集合 $U$ の各要素が $A_1\cap A_2$ に含まれる確率は $p_1p_2=n_1n_2/N^2$ となります。従って、重複数 $x$ は二項分布 $\mathrm{Bin}(N,n_1n_2/N^2)$ に従い、p値はPythonの場合以下で計算されます。
stats.binom.sf(x-1,N,n1*n2/N**2)
フィッシャーの正確検定の近似と考えるよりも、異なる考え方に基づく別の検定であると解釈した方が良いと思います。
数値計算の設定
以下の条件でフィッシャーの正確検定と二項検定を比べました。
| 項目 | 設定 |
|---|---|
| $N$ | 1000 |
| $n_1$ | 100から1000まで100間隔 |
| $n_2$ | 100から1000まで100間隔 |
| 出力 | p<0.05となるために必要な重複数 $x$ |
なお、可能な最大重複 $x=\min(n_1,n_2)$ でもp<0.05にならない場合は0と表記しました。
結果
全体的にかなり似た結果となりました。どちらを使っても実用上は問題なさそうです。二項検定の方が少しだけ値が大きい、つまり、有意になりにくい傾向がありました。
応用
二項検定の利点は拡張がしやすいことです。例えば、3つの集合 $A_1$, $A_2$, $A_3$ があるとき、3つ全ての共通部分が有意に大きいことを示すことが出来ます。$|U|=N$, $|A_i|=n_i (i=1,2,3)$, $|A_1\cap A_2\cap A_3|=x$ とすると、p値は以下で計算出来ます。
stats.binom.sf(x-1,N,(n1/N)*(n2/N)*(n3/N))
同様に、$A_1$ のみに含まれる部分なら、その数を改めて $x$ とおいて、
stats.binom.sf(x-1,N,(n1/N)*(1-n2/N)*(1-n3/N))
などと出来ます。
別の例として、全体集合が異なる場合に適用します。$A_1\subseteq U_1$, $A_2\subseteq U_2$, $U_1\neq U_2$ として、$|U_1|=N_1$, $|U_2|=N_2$, $|U_1\cap U_2|=M$, $|A_1|=n_1$, $|A_2|=n_2$, $|A_1\cap A_2|=x$ とおくと、p値は (多分) 以下で計算できます。
stats.binom.sf(x-1,M,(n1/N1)*(n2/N2))
追記: カイ2乗検定との比較
分割表に対するカイ2乗検定との比較も追加します。コメントを下さった某先生に感謝いたします。p値はPythonの場合以下で計算できます。
p = stats.chi2_contingency([[x,n1-x],[n2-x,N-n1-n2+x]])[1]
上と同じ条件で計算した結果を以下に示します。
フィッシャーの正確検定とほぼ同じでした。もちろん、カイ2乗検定は近似的な手法のため、分割表のいずれかのマスの値が小さいときは不正確になる可能性があります。
カイ2乗検定は行数や列数が2より大きい分割表にも使えます。二項検定でも可能です。行数を $r$、列数を $c$ としたとき、 行和を $n_1,\ldots,n_r$, 列和を $m_1,\ldots,m_c$、総数を $N=\sum n_i=\sum m_j$ とすると、$i$ 行 $j$ 列の観測数 $x$ に対する二項検定のp値は (多分) 以下で計算出来ます。
stats.binom.sf(x-1,N,ni*mj/N**2)
$N$ 個の各要素が、$r$ 通りの値を取るカテゴリカル変数 $X$ と $c$ 通りの値を取るカテゴリカル変数 $Y$ の2つを持つと考えると分かりやすいと思います。$X$ と $Y$ が独立なら $p(X=i,Y=j)=p(X=i)p(Y=j)=(n_i/N)\times(m_j/N)$ が成り立ちます。