裏 RjpWiki

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

Julia の ChisqTest がちょっと変

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

Julia の HypothesisTests で,ChisqTest の結果表示が不適切であることを示す。 実はこの記事は FisherExactTestについて書こうとしていたのだが, つまらないところで引っかかったということである。 まず,2 x 2 分割表の場合である。

Julia> x = [a b; c d]
   
2×2 Array{Int64,2}:
 10   5
  4  20

  4  20

R で実行してみると,以下のようになる。

chisq.test(matrix(c(10, 5, 4, 20), 2))
Pearson's Chi-squared test with Yates' continuity correction
data: matrix(c(10, 5, 4, 20), 2) X-squared = 7.9734, df = 1, p-value = 0.004747

chisq.test(matrix(c(10, 5, 4, 20), 2), correct=FALSE)
Pearson's Chi-squared test
data: matrix(c(10, 5, 4, 20), 2) X-squared = 10.029, df = 1, p-value = 0.001541

Julia でも,簡単にできるが,結果がちょっとおかしい。
   
Julia> obj1 = ChisqTest(x)
   
Pearson's Chi-square Test
   
-----------------
   
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.138067, 0.220907, 0.246548, 0.394477]
    point estimate:          [0.25641, 0.102564, 0.128205, 0.512821]
    95% confidence interval: [(0.1282, 0.4338), (0.0, 0.28), (0.0, 0.3056), (0.3846, 0.6902)]
Test summary:
    outcome with 95% confidence: reject h_0
    one-sided p-value:           0.0015
Details:
    Sample size:        39
    statistic:          10.028571428571428
    degrees of freedom: 1
    residuals:          [1.98898, -1.57243, -1.48842, 1.1767]
    std. residuals:     [3.16679, -3.16679, -3.16679, 3.16679]

まず,p-value のところであるが,one-sided p-value: 0.0015 となっている。
2 x 2 分割表の場合は,連続性の修正のあるなしで結果が違うのではあるが,上のRの結果では,p-value = 0.001541 となっている。
そもそも one-sided p-value と表示されているが,χ2検定には片側検定はない。
Julia では「P値はχ2分布の右側確率(上側確率)だから「片側」と書きました」というのか?
以下に示すように,自由度1のχ2分布で10.028571428571428以上の値を取る確率は確かに0.001541305308889962である。
しかし,これは両側検定のP値である。

Julia> using Distributions
Julia> obj2 = Chisq(1)
Julia> 1-cdf(obj2, 10.028571428571428)
   
0.001541305308889962   

Juliaは,検定によって表示されるP値が片側のものであったり,両側のものであったりするので,検定の結果得られるP値を自分の欲しいものに変えてくれるという便利な(本当か?)pvalue という関数がある。 Juliaが「one-sided p-value」と言うなら,これを使って両側検定のP値を求めようとした人は何を得るだろうか? 両側を指定するのはtail=:bothというから,以下のようにしたら,2倍された値が表示される。

Julia> pvalue(obj1, tail=:both)
   
0.003082610617779839

ちなみに,tail に :left と :right を指定してみると以下のようになる。ちょっと,驚きを禁じ得ない。

Julia> pvalue(obj1, tail=:left)
   
0.99845869469111
   
Julia> pvalue(obj1, tail=:right)
   
0.0015413053088899195

まあ,2 x 2 分割表の場合のχ2検定は,2群の比率の差の検定と p 値は同じになる。
2群の比率の差の検定の検定統計量Zの二乗が 2 x 2 分割表の場合のχ2検定の検定統計量に等しくなるからである。
両者が違うのは,2群の比率の差の検定は片側検定があり得るということである。
どうもChisqTestの作者は,このあたりのことがわかっていないようだ。
そこで,2 x 2 分割表より大きい分割表でのChisqTestの振る舞いをチェックしてみよう。 以下のようなデータで試す。

Julia> x = [1,2,1,2,3,2,1,2,2,2,1,2,3,3,2,1,2,2,3,2,1,2];
Julia> y = [2,1,2,3,2,2,1,2,2,3,2,2,1,2,2,1,2,2,3,2,1,2];

2変数の観測ベクトルがあるとき,分割表は StatsBaseのcounts 関数で集計する。

Julia> using StatsBase
Julia> tbl = counts(x, y) # StatsBase.counts
   
3×3 Array{Int64,2}:
 3  3  0
 1  9  2
 1  2  1
 1  2  1
   
Julia> obj2 = ChisqTest(tbl) # HypothesisTests.ChisqTest
   
Pearson's Chi-square Test
   
-----------------
   
Population details:
    parameter of interest:   Multinomial Probabilities
    value under h_0:         [0.0619835, 0.123967, 0.0413223, 0.173554, 0.347107, 0.115702, 0.0371901, 0.0743802, 0.0247934]
    point estimate:          [0.136364, 0.0454545, 0.0454545, 0.136364, 0.409091, 0.0909091, 0.0, 0.0909091, 0.0454545]
    95% confidence interval: [(0.0, 0.345), (0.0, 0.2541), (0.0, 0.2541), (0.0, 0.345), (0.2273, 0.6178), (0.0, 0.2996), (0.0, 0.2087), (0.0, 0.2996), (0.0, 0.2541)]
Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.2998
Details:
    Sample size:        22
    statistic:          4.8801587301587315
    degrees of freedom: 4
    residuals:          [1.4013, -1.04592, 0.0953463, -0.418718, 0.493464, -0.341882, -0.904534, 0.284268, 0.615457]
    std. residuals:     [1.86926, -1.7648, 0.119913, -0.814215, 1.21376, -0.626783, -1.14133, 0.453705, 0.732163]

片側検定手法がない場合にも,「one-sided p-value」を表示してしまった。
R では,

> chisq.test(matrix(c(3, 3, 0, 1, 9, 2, 1, 2, 1), ncol=3, byrow=TRUE)) Pearson's Chi-squared test data: matrix(c(3, 3, 0, 1, 9, 2, 1, 2, 1), ncol = 3, byrow = TRUE) X-squared = 4.8802, df = 4, p-value = 0.2998

となる。

単なる表示上の誤りということであろうが,統計学をよく知らない人が Julia の表示結果を見て誤解しないよう祈るばかりである。

コメント