Julia の HypothesisTests で，FisherExactTest のデフォルトの結果表示が不適切であることを示す。
Julia が取り扱えるのは 2 x 2 分割表の場合ダケである。残念！！
a = 13; b = 8; c = 34; d = 5;
R で実行してみると，
> options(digits=15)
> fisher.test(matrix(c(13, 8, 34, 5), 2))
Fisher's Exact Test for Count Data
data: matrix(c(13, 8, 34, 5), 2)
p-value = 0.045498106171
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.0525509475953605 1.0314128734283439
sample estimates:
odds ratio
0.245559700550878
Julia でやってみる。
using HypothesisTests
obj = FisherExactTest(13, 8, 34, 5)
Fisher's exact test
-------------------
Population details:
parameter of interest: Odds ratio
value under h_0: 1.0
point estimate: 0.245543
95% confidence interval: (0.0525, 1.0314)
Test summary:
outcome with 95% confidence: fail to reject h_0
two-sided p-value: 0.0561
Details:
contingency table:
13 8
34 5
おや，違いがでました。
R では p-value が 0.045498106171 でしたが，，
Julia の ExactFisherTest では 0.05606237951881565 と出ました。
ずいぶん違いますねえ，Python ではどうなんでしょうか，川崎さん？（川崎さんって誰？）
はい，こちら川崎です。Python では ...
やけに簡単に答えが表示されました。
>>> from scipy import stats
>>> stats.fisher_exact([[13, 8], [34, 5]])
(0.23897058823529413, 0.04549810617099971)
ということです。返される2つのあたいは，オッズ比とP値ということですから...， R と同じようですね...
どーなんでしょ？？？？
この説明のために，フィッシャーの正確確率検定のプログラムを書いてみる。
最終結果だけではなく，P値がどのように計算されるか，途中経過を記録する。
出力中の Prob の列が，四分表の a, b, c, d がそれぞれの値を取ったときの，そのような四分表の生起確率である。
using SpecialFunctions, Printf
function Fisher(a, b, c, d)
function lchoose(n, k)
return lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1)
end
function Stats(a, e, f, g, n)
return exp(lchoose(e, a) + lchoose(f, g - a) - lchoose(n, g))
end
e, f, g, h = a + b, c + d, a + c, b + d
n = e + f
@printf("%5s %5s %5s %5s %15s\n", "a", "b", "c", "d", "Prob")
for a0= max(0, e + g - n):min(e, g)
@printf("%5d %5d %5d %5d %.15f\n", a0, e - a0, g - a0, f - g + a0, Stats(a0, e, f, g, n))
end
end
Fisher(13, 8, 34, 5)
a b c d Prob
8 13 39 0 0.000000039383661
9 12 38 1 0.000002218612928
10 11 37 2 0.000050584374769
11 10 36 3 0.000623873955480
12 9 35 4 0.004679054666097
13 8 34 5 0.022675418766472
14 7 33 6 0.073425165529527
15 6 32 7 0.161535364164963
16 5 31 8 0.242303046247443
17 4 30 9 0.245470406329106
18 3 29 10 0.163646937552738
19 2 28 11 0.068120974005207
20 1 27 12 0.015894893934548
21 0 26 13 0.001572022477043
Julia の FisherExactTest が出力する P 値の下側確率は，0.028031189759407826 である。
pvalue(obj, tail=:left)
0.028031189759407826
この値は，a が 13の場合（観察された四分表） と，13 より小さい場合（12 〜 8 までの四分表）の生起確率の和である。
p0 = +(0.022675418766472, 0.004679054666097, 0.000623873955480, 0.000050584374769, 0.000002218612928, 0.000000039383661)
0.028031189759406997
この結果は片側検定なので，その結果を2倍して両側検定のP値として出力しているわけである。
2 * p0
0.056062379518813994
pvalue(obj)
0.05606237951881565
片側検定のP値を2倍して両側検定のP値にするというのは一部の教科書にも書かれているやりかたではある。
しかし，本来 Fisher が提案した方法は，観察された四分表の生起確率と，それより小さい生起確率の和をもって両側検定のP値とするというやりかたである。
この方法では，観察された四分表とは逆の方向にある四分表の生起確率を加えるということになる。
対象となる四分表の生起確率の和は以下の2つである。
p1 = +(0.015894893934548, 0.001572022477043)
0.017466916411591003
そして，得られるのが以下の値で，R の fisher.test() や Python の scipy.stats.fisher_exact() が提示する結果なのである。
p = p0 + p1 # 0.045498106171
0.045498106170998004
FisherExactTest() でこの計算法での P 値を表示させるには，以下のように pvalue() で method=:minlike を指定する。
この指定が，デフォルトになっていないのはちょっと解せない。
pvalue(obj, method=:minlike)
0.04549810617100024 # 一致した
