裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

Julia でフィッシャーの正確確率検定

2021年01月02日 | ブログラミング

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 # 一致した
コメント