Fisherの正確検定

はじめに

Fisher(フィッシャー)の正確検定(Fisher's exact test)は,分割表(クロス集計表)の各行(各列)が独立かどうかを調べる方法です。直接確率法とも呼ばれます。

この方法はFisherが1935年に著した The Design of Experiments という本の序章の次の第2章の最初に出てくる有名な lady tasting tea の問題を解くために使われています。Fisherの帰無仮説の考え方を最初に説明したものとしても有名です。

2×2分割表の検定

2008年12月8日のニュースによれば,麻生内閣の支持率が前回と比べて半減しました。ネットで調べられる限りの結果を私のブログに載せました。これを見ると,20.9%から25.5%と,かなりばらつきがあります。このばらつきは偶然と考えていいでしょうか。

回答数のわかっている調査について,人数に直すと,次のような分割表(contingency table)が得られます。問題は,この表の行ごと(あるいは列ごと)の独立性の問題に焼き直すことができます(発表されているパーセント値は概数ですので,この数値は正確でないかもしれません)。

支持それ以外
読売228863
毎日217814
朝日4561618
NHK284851

ここでは,一番違いの大きい読売とNHKに限定して考えましょう(残りは後で考えます)。

支持それ以外
読売228863
NHK284851

これを考える前に,もっと単純な問題を考えましょう。男女各5人に意見を聞いたところ,次のように賛否が分かれました。男女で差はあるといえるでしょうか。

32
14

さらに別の問題にしてみましょう。壷の中に赤玉4個,白玉6個が入っていました。ここからランダムに5個取り出したところ,赤玉が3個,白玉が2個である確率を求めてください。

赤玉白玉
取り出した32
まだ壷の中14

10個から5個取り出す組合せは $_{10}C_5$ で,そのうち赤3,白2であるのは $_4C_3 \times {}_6C_2$ 通りですから,確率は

\[ \frac{{}_4C_3 \cdot {}_6C_2}{{}_{10}C_5} \]

になります。$_n C_r$ を求めるRの関数 choose() を使えば

> choose(4,3)*choose(6,2)/choose(10,5)
[1] 0.2380952

のように求められます。

どの玉を取り出す確率も同じという帰無仮説のもとに,10個から5個取り出して,これ以上に珍しい事象が起こる確率($p$ 値)を求めましょう。まず,取り出したものの合計が5個という条件を満たす分割表をすべて挙げれば,

赤玉白玉赤玉白玉赤玉白玉赤玉白玉赤玉白玉
取り出した4132231405
まだ壷の中0514233241

になり,それぞれの確率をさきほどと同じ方法で求めると

\[ 0.02380952, 0.2380952, 0.4761905, 0.2380952, 0.02380952 \]

になります。実際に得られた結果はこのうちの2番目の,確率 0.2380952 のものですので,「これか,これより珍しい結果」を「確率が 0.2380952 以下の結果」と解釈すれば,この5通りのうち真ん中の 0.4761905 を除いた4通りが該当します。それらを合計すると $p = 0.5238$ になります。

縦計,横計のことを周辺度数(marginal frequencies,marginal totals)といいます。周辺度数を固定したときの上のような確率分布を超幾何分布(hypergeometric distribution)といいます。

この手順をもっと簡単にしてくれるRの関数 fisher.test() があります。この引数はこの場合2×2の行列 $\begin{pmatrix} 3 & 2 \\ 1 & 4 \end{pmatrix}$ ですが,これはRでは matrix(c(3,1,2,4), nrow=2) と表します(デフォルトでは縦に読みます)ので,次のように入力します。

> fisher.test(matrix(c(3,1,2,4),nrow=2))

	Fisher's Exact Test for Count Data

data:  matrix(c(3, 1, 2, 4), nrow = 2) 
p-value = 0.5238
alternative hypothesis: true odds ratio is not equal to 1 
95 percent confidence interval:
   0.2180460 390.5629165 
sample estimates:
odds ratio 
  4.918388 

この方法をフィッシャーの正確検定(Fisher's exact test)または直接確率法といいます。

上の例では,10個のうちちょうど半分を取り出したので,分布は左右対称でしたが,10個のうち4個を取り出すことにすれば

赤玉白玉赤玉白玉赤玉白玉赤玉白玉赤玉白玉
取り出した4031221304
まだ壷の中0615243342

の確率は

\[ 0.004761905, 0.1142857, 0.4285714, 0.3809524, 0.07142857 \]

となり,対称でなくなります。この場合,左から2番目(確率0.1142857)が実際の結果とすれば,「これか,より珍しい場合」の合計は

\[ 0.004761905 + 0.1142857 + 0.07142857 \]

になります。これは両側検定ですが,片側検定なら 0.004761905 + 0.1142857 です。片側確率の2倍が両側確率にならない例の一つです。

上では「これか,より珍しい場合」を,フィッシャーに従って「確率がこれ以下の場合」と解釈しましたが,「$\chi^2$ 値がこれ以下の場合」と解釈することもできます。両側検定の場合,両者の結果が異なることもあります。また,両側確率は単に片側確率を2倍するのがよいという意見(F. Yates, ``Test of Significance for 2×2 Contingency Tables'', Journal of the Royal Statistical Society, Series A (General), Vol. 147, No. 3 (1984), pp. 426-463)もありますが,2×2 分割表以外(一般に複数の自由度がある場合)では,そもそも片側・両側という概念が定義できません。

最初の問題に戻って,読売とNHKの麻生内閣支持率の違い

支持それ以外
読売228863
NHK284851

を検定してみましょう。

> fisher.test(matrix(c(228,284,863,851),nrow=2))

結果は $p = 0.02335$ で,こんな違いが偶然に生じるのは100回に2回しかないことがわかります。

Rの matrix() はオプション byrow=TRUE を与えないと列ごとに数値を読むので,上では行と列を入れ替えた入力をしていることになりますが,縦計 $a + c$,$b + d$,横計 $a + b$,$c + d$ を固定したときに分割表 $\bigl[\begin{smallmatrix} a & b \\ c & d \end{smallmatrix} \bigr]$ を得る確率は

\[ \frac{{}_{a+c}C_a \cdot {}_{b+d}C_b}{{}_nC_{a+b}} = \frac{\dfrac{(a+c)!}{a!c!} \cdot \dfrac{(b+d)!}{b!d!}}{\dfrac{n!}{(a+b)!(c+d)!}} = \frac{(a+b)!(c+d)!(a+c)!(b+d)!}{n!a!b!c!d!} \]

ですから,行と列を入れ替えても同じ結果になります。

任意の分割表の場合

上の式を拡張して,任意の $m \times n$ 分割表 $[a_{ij}]$ について

\[ \frac{(\sum_j a_{1j})! \ldots (\sum_j a_{mj})! (\sum_i a_{i1})! \ldots (\sum_i a_{in})! }{n! a_{11}! a_{12}! \ldots a_{mn}!} \]

を計算することができます。元の分割表と同じ行和・列和を与えるすべての分割表について上の確率を計算し,元の分割表以下の確率の和を求めれば,任意の分割表についてフィッシャーの方法での $p$ 値が得られます。今のコンピュータは速いので,たいていの場合に,古い $\chi^2$ 検定は不要になります。

このページの最初の表の全データを入れてやってみましょう。

> fisher.test(matrix(c(228,863,217,814,456,1618,284,851),nrow=4,byrow=TRUE))

	Fisher's Exact Test for Count Data

data:  
p-value = 0.0712
alternative hypothesis: two.sided 

あれあれ,読売とNHKの違いが $p = 0.02335$ で有意だったのに,今度は $p = 0.0712$ で有意ではなくなりました(というような 0.05 を境にして結果の解釈が質的に変わってしまうような言い方には問題があるのですが,いずれにしても $p$ 値が減少してしまっています)。これは多重比較(multiple comparison)の問題といって,4行のデータ全体としてのばらつきは有意でなくても,その中から2行を選んで行う $_4C_2 = 6$ 通りの検定のうちどれかは有意になってしまうということはよくあります。7行のデータなら $_7C_2 = 21$ 通りの比較ができるので,何も違わなくても期待値として21通りのうち一つは5%水準で有意になります。多重比較には注意が必要です。

2×2 より大きい分割表の場合,メモリが足りなくなることがあります。このときは,fisher.test()workspace というオプションにメモリの量(4バイト単位)を指定します。デフォルトは workspace = 200000 です。これ以外にもいろいろなオプションがありますので,うまくいかないときは help() で調べてください。

オッズ比

fisher.test() はオッズ比(odds ratio,OR)も出力します。オッズ(odds)とは「当たりの確率」を「外れの確率」で割ったもので,オッズ比はその比です。内閣支持率の例では (354/432)/(418/370) = 0.725 がオッズ比です。分割表の縦横を入れ替えてもオッズ比は同じです。詳しくは2×2の表,オッズ比,相対危険度をご覧ください。

fisher.test() の出力するオッズ比は上の定義から少し外れていますので,少し説明しておきます。$a + c$ 個の赤玉から $a$ 個,$b + d$ 個の白玉から $b$ 個を取り出す分割表 $\bigl[\begin{smallmatrix} a & b \\ c & d \end{smallmatrix} \bigr]$ で,赤・白を取り出す確率 $p$,$q$ が等しくないとし,$a$ の分布は ${}_{a+c}C_a p^a(1-p)^c$,$b$ の分布は ${}_{b+d}C_b q^b(1-q)^d$ という2項分布でモデル化します。$a$,$c$ が与えられたとき,確率 ${}_{a+c}C_a p^a(1-p)^c$ を最大にするパラメータ $p$ を求めてみましょう。このように確率を最大にするように選んだパラメータを最尤推定量(maximum likelihood estimator,MLE)といいます。確率を $p$ で微分したものを 0 と置くと,$p/(1-p) = a/c$ という期待通りの式が出ます。同じように白玉についても $q/(1-q) = b/d$ で,両者の間に何の条件も付けなければ,オッズ比は期待通りの $\dfrac{p/(1-p)}{q/(1-q)} = \dfrac{a/c}{b/d}$ という式で求められます。ここで,赤玉を $a$ 個,白玉を $b$ 個選ぶ確率は,積

\[ {}_{a+c}C_a p^a(1-p)^c \times {}_{b+d}C_b q^b(1-q)^d \]

で与えられます。これはすべての可能な $a$,$b$ について合計すれば1になります。これを,取り出す個数 $a + b$ を固定したときの条件付き確率にするには,$a + b$ を固定してすべての可能な $a$ について合計すると1になるように比例定数を付け替えなければなりません。見通しをよくするために上の式を

\[ {}_{a+c}C_a p^a(1-p)^{-a} \times {}_{b+d}C_b q^{-a}(1-q)^a \times (1-p)^{a+c} q^{a+b} (1-q)^{(b+d)-(a+b)} \]

と変形すると,最後の $\times$ 以下は定数ですので,条件付き確率は結局

\[ \frac{{}_{a+c}C_a \cdot {}_{b+d}C_b \cdot \omega^a}{\sum\limits_a {}_{a+c}C_a \cdot {}_{b+d}C_b \cdot \omega^a}, \qquad \omega = \frac{p/(1-p)}{q/(1-q)} \]

になります。これを最大にする $\omega$ が,Rが出力するオッズ比です。これは $a + b =$ 一定 という条件付きの最尤推定量(conditional MLE)です。通常はオッズ比として,条件なしの最尤推定量(unconditional MLE)である $\dfrac{a/b}{c/d}$ を報告すればいいでしょう。

Fisherの検定はオッズ比が1であることの検定と理屈上は同じことですが,fisher.test() の採用する一般的な計算方法では,$p$ 値とオッズ比の信頼区間とは別の方法で計算されるので,同じ結果にならないことがあります。例えば

ex1 = matrix(c(6,12,12,5), nrow=2)

という例で試してみましょう。fisher.test(ex1) では,$p$ 値は 0.044 ですので5%水準で有意ですが,オッズ比の95%信頼区間は $[0.039, 1.056]$ ですので 1 を含みます。

このあたりの事情は exact2x2 というパッケージの解説に書かれている通りです。この状況を改善するには,次のように exact2x2 パッケージをインストールし,fisher.test() の代わりに fisher.exact() を使います:

install.packages("exact2x2")
library(exact2x2)
fisher.exact(ex1)

$p$ 値は fisher.test() と同じ 0.044 ですが,95%信頼区間は $[0.0435, 0.9170]$ となり,5%水準で有意でないことと矛盾しない結果になります。

参考文献:Michael P. Fay, Two-sided Exact Tests and Matching Confidence Intervals for Discrete Data, R Journal 2, 53-58 (2010). [PDF]

Fisherの正確検定は正確か?

昔からいろいろ議論があるところです。昔の議論は誤解に基づくものもありましたが,今でもシミュレーションしてみてカイ2乗検定のほうがいいではないかという話がよく出ます(Fisherの正確検定かカイ2乗検定か参照)。

まずこのページの上のほうにある「壷の中に赤玉4個,白玉6個が入っていました。ここからランダムに5個取り出したところ」という問題と,その上の「男女各5人に意見を聞いたところ」という問題は,少し違うものです。壷の問題は周辺度数が固定されているので超幾何分布を使うことに異論はないでしょうけれども,男女各5人の賛否は 4:6 に固定されているわけではなく,潜在的に起こりうる事象はもっとたくさんあります。このあたりに,カイ2乗検定のほうが正確な理由があります。

同様なことをBayesianの立場から論じた論文がGelmanのブログ I hate the so-called Fisher exact test: a pointer to one of my favorite articles - Statistical Modeling, Causal Inference, and Social Science からリンクされています。