フィッシャーの正確確率検定　カイ二乗検定　3群以上の比較　Rで実施する方法

フィッシャーの正確確率検定、カイ二乗検定の3群以上の比較をRで実施する方法の解説。

＞＞もう統計で悩むのを終わりにしませんか？

↑1万人以上の医療従事者が購読中

カイ二乗検定の3群以上の比較
フィッシャー正確確率検定の3群以上の比較
- フィッシャーの正確確率検定とカイ二乗検定の結果の比較
まとめ
参考書籍

カイ二乗検定の3群以上の比較

三群以上の割合の比較はどうやればいいのか？

Bonferroni型のp値調整を使う方法がある。

toukeier.hatenablog.com

Rで行う場合、pairwise.prop.test()という関数を使う。

カイ二乗検定　3群以上の比較で多重調整しない場合

例として下記の2x3の分割表のデータでやってみる。

	Group1	Group2	Group3
No	49 (100%)	97 (96%)	63 (86%)
Yes	0 (0%)	4 (4%)	10 (14%)

パーセントは列パーセントを記している。

分割表作成と検定のためのRスクリプトは以下の通り。

tab <- matrix(c(49,0,97,4,63,10),nr=2)

dimnames(tab) <- list(c("No","Yes"),c(paste("Group",1:3,sep="")))

tab

round(prop.table(tab, 2),2)

> tab
    Group1 Group2 Group3
No      49     97     63
Yes      0      4     10
> 
> round(prop.table(tab, 2),2)
    Group1 Group2 Group3
No       1   0.96   0.86
Yes      0   0.04   0.14

多重比較の結果、Group1と3、Group2と3がそれぞれ0.05を下回って、統計学的有意に異なるように見える。

> pairwise.prop.test(t(tab), p.adj="none")

        Pairwise comparisons using Pairwise comparison of proportions 

data:  t(tab) 

       Group1 Group2
Group2 0.383  -     
Group3 0.018  0.041 

P value adjustment method: none 
Warning messages:
1: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
  Chi-squared approximation may be incorrect
2: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
  Chi-squared approximation may be incorrect

例題は期待値が5以下のマス目が20%以上なので、 $\chi^2$ 検定が不適切であるメッセージが出ている。

$\chi^2$ 検定が不適切な場合は、フィッシャーの正確確率検定（後述）を行う。

以下は、期待値を計算した分割表。

	Group1	Group2	Group3	Total
No	45.9	97.4	68.4	209
Yes	3.1	6.3	4.6	14
Total	49	101	73	223

赤で示した二つのセルが期待値5未満だ。

どのグループ間比較でも、期待値5未満のセルが一つまたは二つ含まれる。

一つのセルだけでも $1/4=25 \%$ 、二つのセルなら $2/4=50 \%$ なので、いずれの場合も $20 \%$ を超えている。

ゆえに、 $\chi^2$ 検定を使うのは適切ではない。

＞＞もう統計で悩むのを終わりにしませんか？

↑1万人以上の医療従事者が購読中

カイ二乗検定　ホルムの多重調整

ホルムの方法は、もっとも小さいp値に最大比較ペア数をかける。

例はペア数は3なので、3をかける。

2番目のp値に2をかける。一番大きいp値はそのまま。と計算される。

1と3： $0.018 \times 3 = 0.054$

2と3： $0.041 \times 2 = 0.081$

もっとも小さいp値の3倍が0.05を下回らなかった時点で、計算は終了だ。

多重比較調整するとどの群間も統計学的有意に異ならない。

> pairwise.prop.test(t(tab))

        Pairwise comparisons using Pairwise comparison of proportions 

data:  t(tab) 

       Group1 Group2
Group2 0.383  -     
Group3 0.054  0.081 

P value adjustment method: holm 
Warning messages:
1: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
  Chi-squared approximation may be incorrect
2: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
  Chi-squared approximation may be incorrect

カイ二乗検定　ホックバーグの多重調整

もう一つの方法はホックバーグの方法。

こちらのほうが検出力が高く、おすすめ。

ホルムと同様に、大きいp値はそのまま、2番目に大きいp値を2倍して、一番小さいp値を3倍する。

チェックの仕方がホルムの時と逆で、大きいp値から始める。

もしもっとも大きいp値が0.05を下回っていたら、それ以外のp値がどうであれ、すべて統計学的有意に異なると判断する。

もっとも大きいp値が0.05を下回っていなかったら、次に2番目に大きいp値の2倍と0.05を比較する。

ここでもし0.05より小さければ、最も小さいp値は不問で、ともに統計学的有意に異なると判断する。

今回の例は、どのレベルでも0.05を下回らなかったので、どの群間にも統計学的有意差はなかった。

> pairwise.prop.test(t(tab),p.adjust="hochberg")

        Pairwise comparisons using Pairwise comparison of proportions 

data:  t(tab) 

       Group1 Group2
Group2 0.383  -     
Group3 0.054  0.081 

P value adjustment method: hochberg 
Warning messages:
1: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
  Chi-squared approximation may be incorrect
2: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
  Chi-squared approximation may be incorrect

ホルムは小さいp値から0.05を下回り続けなければならない点で厳しく、ホックバーグは、大きいp値が一回0.05を下回った段階でそれより小さいp値は不問で全部有意というところが、検出力が高いゆえんだ。

フィッシャー正確確率検定の3群以上の比較

フィッシャーの正確確率検定を3群以上の比較で使う方法。

RVAideMemoireパッケージのfisher.multcomp()を使う。

RVAideMemoireパッケージをインストールする。

インストールは一回だけ。

install.packages("RVAideMemoire")

p.methodで多重比較調整の方法を指定する。

library(RVAideMemoire)

fisher.multcomp(tab, p.method="none")

fisher.multcomp(tab, p.method="holm")

fisher.multcomp(tab, p.method="hochberg")

フィッシャーの正確確率検定とカイ二乗検定の結果の比較

このデータはフィッシャーの正確確率検定を使うことが適切である典型例と思う。

	$\chi^{2}$	$\chi^{2}$ with Holm/Hochberg	Fisher	Fisher with Holm/Hochberg
Group1:Group2	0.383	0.383	0.3038	0.3038
Group1:Group3	0.018	0.054	0.005656	0.01697
Group2:Group3	0.041	0.081	0.02481	0.04962

Group1と3、Group2と3のp値が、 $\chi^2$ に比べるとFisherのほうが断然低い。

ホルムなら、Group1と3が有意、Group2と3が有意、Group1と2は有意でなく計算終了。

ホックバーグなら、Group1と2は有意でないが、Group2と3が有意で、Group1と3を確認する前に終了。

結果として、Group1と3、Group2と3の群間に統計学的有意な差を認めた（赤で示した有意確率）。

つまり、Group3がダントツでYesの割合が高いという結論だ。

以下は、上記の結果を計算するためのRスクリプト。

> library(RVAideMemoire)
> 
> fisher.multcomp(tab, p.method="none")

        Pairwise comparisons using Fisher's exact test for count data

data:  tab

       Group1:Group2 Group1:Group3 Group2:Group3
No:Yes        0.3038      0.005656       0.02481

P value adjustment method: none
> 
> fisher.multcomp(tab, p.method="holm")

        Pairwise comparisons using Fisher's exact test for count data

data:  tab

       Group1:Group2 Group1:Group3 Group2:Group3
No:Yes        0.3038       0.01697       0.04962

P value adjustment method: holm
> 
> fisher.multcomp(tab, p.method="hochberg")

        Pairwise comparisons using Fisher's exact test for count data

data:  tab

       Group1:Group2 Group1:Group3 Group2:Group3
No:Yes        0.3038       0.01697       0.04962

P value adjustment method: hochberg

ちなみに、fisher.multcomp()は $2 \times n$ を超え $m \times n$ の分割表に適用できる。

下記の結果から抜き出すと、検討に値する組み合わせは以下の5通り。

H:MでGp1:Gp3　0.848485
H:MでGp1:Gp4　0.535714
M:LでGp1:Gp2　0.05031
M:LでGp1:Gp4　0.008772
M:LでGp1:Gp3　0.005409

結論として、

MとL x Gp1とGp4
MとL x Gp1とGp3

この2つの組み合わせが統計学的有意に異なっていた。

ちなみに、調整したp値が1を超える場合、1と表現されている。

> tab1 <- matrix(c(0,6,2,1,3,19,3,2,27,2,0,13),nr=3)

> dimnames(tab1) <- list(c("H","M","L"),c(paste("Gp",1:4,sep="")))
> tab1
  Gp1 Gp2 Gp3 Gp4
H   0   1   3   2
M   6   3   2   0
L   2  19  27  13

> fisher.multcomp(tab1, p.method="none")

        Pairwise comparisons using Fisher's exact test for count data

data:  tab1

     Gp1:Gp2   Gp1:Gp3  Gp1:Gp4 Gp2:Gp3 Gp2:Gp4 Gp3:Gp4
H:M 0.400000 0.0606061 0.035714  0.5238  0.4000       1
H:L 1.000000 1.0000000 1.000000  0.6411  0.5646       1
M:L 0.003145 0.0003005 0.000516  0.6407  0.2790       1

P value adjustment method: none

> fisher.multcomp(tab1, p.method="hochberg")

        Pairwise comparisons using Fisher's exact test for count data

data:  tab1

    Gp1:Gp2  Gp1:Gp3  Gp1:Gp4 Gp2:Gp3 Gp2:Gp4 Gp3:Gp4
H:M 1.00000 0.848485 0.535714       1       1       1
H:L 1.00000 1.000000 1.000000       1       1       1
M:L 0.05031 0.005409 0.008772       1       1       1

P value adjustment method: hochberg