フィッシャーの検定とモンテカルロ法①
Rでフィッシャーの正確検定（fisher.test()）を行おうと思ったら、
Try increasing the size of the workspace
なるエラーメッセージがでてきた。
よくわからんが、とにかくworkspace=1e50 とか試して、workspaceなるものを増やしてみたが結局うまくいかなかった。
help(fisher.test)で、その説明を確認すると、どうやら2×2より大きいマトリックスに対してフィッシャーの正確検定を行う場合、simulate.p.valueをTRUEに設定してあげると、モンテカルロ法でp値が計算され、うまくいく。
友人もこのような出来事に出会っていた ＆ TAを担当している授業でSPSSでフィッシャーの正確検定ができないことが話題となっていたので、フィッシャーの正確検定およびモンテカルロ法についてまとめる。
記事が長くなってしまうので、今回はフィッシャーの正確検定を使う必要性がなぜあるのかについてのみ書く。
フィッシャーの正確検定はなぜ必要か
実用的には、フィッシャーの正確検定は、セル度数が少ない時にカイ二乗検定の代わりに使用される。（明確な基準はないだろうが、セル度数が5以下の場合はフィッシャーの正確検定をやっておけと教わった）
セル度数が少ない場合、カイ二乗検定を行う際に、以下の計算で使用される統計量（カイ二乗統計量と呼ばれたりするが、ややこしいので単に統計量と呼ぶ）がカイ二乗分布に従わなくなってしまうことが背景にある。
なお、ここでEは期待度数、Oは観測度数である。カイ二乗分布については以前に記事を書いたのでそちらを参照されたい。
abcxyzonetwothree.hatenablog.com
今回の記事では、本当に上記の統計量が、セル度数が少ない場合にカイ二乗分布に従わなくなってしまうのかについて実験してみる。
方法は以下の通り。
①異なる二項分布（0 or 1）から乱数を発生させる。
②0かつ0である場合はマトリックスの1行1列目、1かつ0である場合はマトリックスの1行2列目、0かつ1である場合はマトリックスの2行1列目、1かつ1である場合はマトリックスの2行2列目に収納する
これは、全く独立な二項分布からの乱数なので、互いに関連していないと言える。
以上①〜③を1万回繰り返す。互いに独立なものなので、この時の統計量はカイ二乗分布に従うはずである。
なお、①の時に何個の乱数を発生させるかを変化させることで、セル度数の数を少なくしたり多くしたりする。例えば、10個の乱数しか発生させない場合、セルの数は4つあるので、平均2.5にしかそのセル度数はならない。よって、この場合はカイ二乗分布に従わないことが予想される。一方で、例えば、1000個の乱数を発生させる時、セルの数は4つあるので、平均250が各セルの度数となる。よって、この場合はカイ二乗分布に従うと予想される。
今回は、乱数の数を10、100、1000の3通りで実験する。
結果は以下の通り。なお、赤色のヒストグラムが10個の乱数を発生させた場合（平均セル度数=2.5）、青色が100個の乱数を発生させた場合（平均セル度数=25）、水色のヒストグラムが1000個の乱数を発生させた場合（平均セル度数=250）、緑色が比較のために自由度1（2×2のマトリックスなので）のカイ二乗分布から1万個の乱数を発生させて、それをヒストグラムにしたものである。
これをみると、赤色のヒストグラムは緑色のヒストグラムと比べて、いびつな形をしていることがわかる。なだらかなスロープになっていないし、0が多い。一方で、青色のヒストグラムは、緑色のヒストグラムと近い形になっている。ただし、0の値が以上に大きいなどの違いはあるので完全ではない。しかし、少なくとも赤色のヒストグラムよりは、緑色のヒストグラムに近い形になっている。次に、水色のヒストグラムをみるとだいぶ緑色のヒストグラムに近い形になっている。少なくとも、赤色、青色のヒストグラムよりも緑色のヒストグラムに近い形になっている。
以上から、セル度数が少なくなればなるほど、カイ二乗分布に従わなくなると言える。
なお、今回の実験のコードは以下の通り。
N<-10000 #試行回数 rec10<-numeric(N) #カイ二乗値の記録場所 rec100<-numeric(N) #カイ二乗値の記録場所 rec1000<-numeric(N) #カイ二乗値の記録場所 for(i in 1:N){ x<-rbinom(10,1,0.5) #二項分布から乱数発生 y<-rbinom(10,1,0.5) #二項分布から乱数発生 data<-data.frame(cbind(x,y)) #x=0かつy=0、x=1かつy=0、x=0かつy=1、x=1かつy=1の4つのセルを持つマトリックスを作成 mat<-matrix(numeric(),nrow=2,ncol=2) a<-ifelse(data$x==0&data$y==0,1,0) b<-ifelse(data$x==1&data$y==0,1,0) c<-ifelse(data$x==0&data$y==1,1,0) d<-ifelse(data$x==1&data$y==1,1,0) mat[1,1]<-sum(a) mat[1,2]<-sum(b) mat[2,1]<-sum(c) mat[2,2]<-sum(d) #カイ二乗テスト chi_test<-chisq.test(mat) #統計量取り出し rec10[i]<-chi_test$statistic ############################### x<-rbinom(100,1,0.5) #二項分布から乱数発生 y<-rbinom(100,1,0.5) #二項分布から乱数発生 data<-data.frame(cbind(x,y)) #x=0かつy=0、x=1かつy=0、x=0かつy=1、x=1かつy=1の4つのセルを持つマトリックスを作成 mat<-matrix(numeric(),nrow=2,ncol=2) a<-ifelse(data$x==0&data$y==0,1,0) b<-ifelse(data$x==1&data$y==0,1,0) c<-ifelse(data$x==0&data$y==1,1,0) d<-ifelse(data$x==1&data$y==1,1,0) mat[1,1]<-sum(a) mat[1,2]<-sum(b) mat[2,1]<-sum(c) mat[2,2]<-sum(d) #カイ二乗テスト chi_test<-chisq.test(mat) #統計量取り出し rec100[i]<-chi_test$statistic ############################### x<-rbinom(1000,1,0.5) #二項分布から乱数発生 y<-rbinom(1000,1,0.5) #二項分布から乱数発生 data<-data.frame(cbind(x,y)) #x=0かつy=0、x=1かつy=0、x=0かつy=1、x=1かつy=1の4つのセルを持つマトリックスを作成 mat<-matrix(numeric(),nrow=2,ncol=2) a<-ifelse(data$x==0&data$y==0,1,0) b<-ifelse(data$x==1&data$y==0,1,0) c<-ifelse(data$x==0&data$y==1,1,0) d<-ifelse(data$x==1&data$y==1,1,0) mat[1,1]<-sum(a) mat[1,2]<-sum(b) mat[2,1]<-sum(c) mat[2,2]<-sum(d) #カイ二乗テスト chi_test<-chisq.test(mat) #統計量取り出し rec1000[i]<-chi_test$statistic } #結果の図示 x<-rchisq(N,df=1) par(mfrow=c(2,2)) hist(rec10,breaks=seq(0,17,0.5),ylim=c(0,10000),col="red") hist(rec100,breaks=seq(0,17,0.5),ylim=c(0,10000),col="blue") hist(rec1000,breaks=seq(0,17,0.5),ylim=c(0,10000),col="lightskyblue") hist(x,breaks=seq(0,17,0.5),ylim=c(0,10000),col="green")
この実験から、セル度数が少ない時に、その統計量がカイ二乗分布に従わなくなる可能性が示唆された。
今回は以上である。
次回は、セル度数が少ない時に、その使用が推奨されるフィッシャーの検定そのものについて書く。
p.s.
厳密に言えば、今回の実験は「セル度数が少ないから」とは言えない。本来は、サンプル数を固定した上でセル度数だけを変更して、実験するべきである。ただし、関連がないという仮定のもとで、そのような操作を行うのは、私の頭ではどのようにコードを書けば良いかわからなかったので、次善の策として今回のような実験を行なった。誰か教えてください。
