フィッシャーの検定とモンテカルロ法②

前回からの続き。前回の記事は以下。

abcxyzonetwothree.hatenablog.com

今回は、フィッシャーの検定について。

例えば、以下のようなデータを得たとしよう。

この時、縦と横の周辺度数が固定されたとして、あり得るマトリックスは何だろうか。考えられるのは以下の4つのマトリックスである。

この時、「幸せか不幸せか」「金持ちか貧乏か」という2つの事象が独立であるという帰無仮説のもとで、それぞれのマトリックスが出現する確率は何であろうか。

フィッシャー大先生は、そのマトリックスが出現する回数が帰無仮説のもとで以下のようになることを示したようである。

なお、セル（1,1）をa、セル（1,2）をb、セル（2,1）をc、セル（2,2）をdとしている。

 これに従うと、mat1〜mat4に出る確率はそれぞれ以下のようになる。

なお、確率の合計値はちゃんと1になる。

ここで、もともとのデータはmat3に該当する。mat3が出る確率は、0.07142857である。これ以下の確率で出現するマトリックスは他にmat4（その確率は、0.07142857）が存在する。従って、mat3が出るか、mat3が出る確率以下の確率で出現するマトリクッス（今回の場合は、mat4）が出てくる可能性は0.07142857+0.07142857で、0.1428571である。これがフィッシャーの正確検定のp値である。これが、0.05未満だと当該マトリクッスか、それ以下の確率で出てくるマトリクッスが無視できるほど小さい（有意）ので、帰無仮説を棄却することになる。

以上がフィッシャーの正確検定の説明である。

次に、フィッシャーの正確検定とカイ二乗検定がどれくらい検定結果が異なるのかを実験してみた。実用的に困るのは、一方が有意になり、片方が有意にならない場合である。今回は、10%水準で実験してみた。手順は以下の通りである。

①異なる2つの二項分布（0 or 1）から乱数を発生させる。

②0かつ0である場合はマトリックスの1行1列目、1かつ0である場合はマトリックスの1行2列目、0かつ1である場合はマトリックスの2行1列目、1かつ1である場合はマトリックスの2行2列目に収納する

③このマトリックスでカイ二乗検定およびフィッシャーの正確検定を行う

以上①〜③を1万回繰り返す。これで、片方が有意（p<0.1）で片方が有意でない場合を問題あり、両方とも有意もしくは両方とも有意ではない場合は問題なしと判断し、1万回中問題ありとどれくらい判断されるのかを確認する。

なお、①の時に何個の乱数を発生させるかを変化させることで、セル度数の数を少なくしたり多くしたりする。例えば、10個の乱数しか発生させない場合、セルの数は4つあるので、平均2.5にしかそのセル度数はならない。一方で、例えば、1000個の乱数を発生させる時、セルの数は4つあるので、平均250が各セルの度数となる。

今回は乱数の数を10、100、1000の3通りで実験する。

まず、乱数の数を10とした場合、問題なしと判断されたのは8923回、問題ありと判断されたのは1032回である（なお、セル度数が小さすぎるなどの問題で45回検定が行えなかった）。すなわち、おおよそ10%の確率で問題が発生してしまう。

次に、乱数の数を100とした場合、問題なしと判断されたのは9717回、問題ありと判断されたのは283回である。すなわち、2.83%の確率で問題が発生した。

最後に、乱数の数を1000とした場合、問題なしと判断されたのは9948回、問題ありと判断されたのは52回である。すなわち、0.52%の確率で問題が発生した。

ここから言えるのは、やはりセル度数が少ない時は、フィッシャーの正確検定とカイ二乗検定の結果が異なる傾向にあるということである。フィッシャーの正確検定が「良い」検定かはわからないが、少なくともセル度数が少ない時に両者の結果が異なってしまう。この結論は、前回の記事の、セル度数が少ない時に（カイ二乗）統計量がカイ二乗分布に従わなくなってしまう結果とも整合的である。

以上である。

次は、いよいよモンテカルロ法について。

なお、今回の実験のコードは以下の通り。

#p値の違いの判定をする関数（0=問題なし、1=問題あり、2=どちらかがp=0.1もしくは欠損値の場合）
test<-function(L,R){
  if(is.na(L)==T){L<-0.1}
  if(is.na(R)==T){R<-0.1}
  if(L < 0.1 && R < 0.1){result<-0}
  if(L > 0.1 && R > 0.1){result<-0}
  if(L < 0.1 && R > 0.1){result<-1}
  if(L > 0.1 && R < 0.1){result<-1}
  if(L==0.1 | R==0.1){result<-2}
  result
}

################################

N<-10000 #試行回数
rec_chi_p<-numeric(N)
rec_fisher_p<-numeric(N)
rec<-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<-chisq.test(mat,correct = F) #カイ二乗検定
  fisher<-fisher.test(mat) #フィッシャーの正確検定
  
  rec_chi_p[i]<-chi$p.value
  rec_fisher_p[i]<-fisher$p.value
  
  rec[i]<-test(rec_chi_p[i],rec_fisher_p[i])
}
table(rec)

################################

N<-10000 #試行回数
rec_chi_p<-numeric(N)
rec_fisher_p<-numeric(N)
rec<-numeric(N) #記録場所

for(i in 1:N){
  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<-chisq.test(mat,correct = F) #カイ二乗検定
  fisher<-fisher.test(mat) #フィッシャーの正確検定
  
  rec_chi_p[i]<-chi$p.value
  rec_fisher_p[i]<-fisher$p.value
  
  rec[i]<-test(rec_chi_p[i],rec_fisher_p[i])
}
table(rec)

################################

N<-10000 #試行回数
rec_chi_p<-numeric(N)
rec_fisher_p<-numeric(N)
rec<-numeric(N) #記録場所

for(i in 1:N){
  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<-chisq.test(mat,correct = F) #カイ二乗検定
  fisher<-fisher.test(mat) #フィッシャーの正確検定
  
  rec_chi_p[i]<-chi$p.value
  rec_fisher_p[i]<-fisher$p.value
  
  rec[i]<-test(rec_chi_p[i],rec_fisher_p[i])
}
table(rec)