2008-01-20

カイ自乗値とalpha値から相当する自由度を求める

今、カイ自乗値がyのとき、自由度dの棄却検定を行って、タイプ１エラーがaであるとなるとする。
Rの関数を使えば、

a<-pchisq(y,d,lower.tail=FALSE)

という関係である。
今、aとdが与えられているときに、yを求める関数もRにはあって、

y<-qchisq(a,d,lower.tail=FALSE)

である。では、yとaが与えられているときに、dはどのようにして求めるか？
aが0.5未満のとき(厳密には、0.5より『十分に(！：後述、参照)小さいとき』)pchisq関数を用いて、一定の精度で推定する関数は末尾の通り。なお、この推定は次のような考えに基づく。

サンプル数Nの２ｘ２分割表のカイ自乗値の最大値は、２群のサンプル数N1,N2(N1+N2=N)が、きれいに分かれた、{[N1,0],[0,N2]}なる２ｘ２表のときに得られて、その値はNである。
- $\frac{(N1-\frac{N1*N1}{N})^2}{\frac{N1*N1}{N}}+\frac{(0-\frac{N1*N2}{N})^2}{\frac{N1*N2}{N}}+\frac{(0-\frac{N2*N1}{N})^2}{\frac{N2*N1}{N}}+\frac{(N2-\frac{N2*N2}{N})^2}{\frac{N2*N2}{N}}$ を式変形して得られる。
今、このNサンプルが、すべて弁別可能であるとき、Nカテゴリに２群のうちのどちらか１群のサンプルが１つあてはまるような、２ｘN分割表であるとみなせる。この分割表では、第１群に１がN1個、０がN2個、第２群に１がN2個、０がN1個あり、周辺度数は、２群のそれぞれの和がN1、N2であって、Nカテゴリのそれぞれに各１サンプルがあるような、そんな分割表である。この分割表のカイ自乗値もNである。
- $N1*\frac{(1-\frac{N1*1}{N})^2}{\frac{N1*1}{N}}+N2*\frac{(0-\frac{N1*1}{N})^2}{\frac{N1*1}{N}}+N2*\frac{(1-\frac{N2*1}{N})^2}{\frac{N2*1}{N}}+N1*\frac{(1-\frac{N2*1}{N})^2}{\frac{N2*1}{N}}$ を式変形して得られる
これは、サンプル数がNであるとき、自由度１のカイ自乗検定ではカイ自乗値の最大値がサンプル数Nに等しいことを示すとともに、自由度N(全員を弁別できるだけのカテゴリ数があるとき)のときに、そのカイ自乗値もNであることを意味している。全員が弁別できているときというのは、何も統計的に有意なことはいえないときであり、検定aは0.5になる。このことは、Rで

n<-seq(from=1,to=100,by=0.5)
N<-n^2
plot(N,pchisq(N,N,lower.tail=FALSE))

のようにすることで、pchisq(N,N,lower.tail=FALSE)が0.5に収束することからも確認できる。
これを利用して、aが0.5未満のとき(pchisq(N,N,lower.tail=FALSE)の0.5への収束を考慮に入れて、aが、pchisq(N,N,lower.tail=FALSE)未満のとき)に、 $a=pchisq(y,d,lower.tail=FALSE)$ を満足するようなdを１からyの範囲で挟み撃ちにて、探索すれば、カイ自乗値とその検定aの値とから、それを与える自由度dが求まるものである。
引数は、探索の範囲がx1とx2(上述の通り、x1=1,x2=y)、yはカイ自乗値、aはタイプ１エラー、thresは、探索終了のための、dの精度を与えている

DF<-function(x1,x2,y,a,thres){
 low=x1
 mid=(x1+x2)/2
 high=x2
 L=pchisq(y,low,lower.tail=FALSE)
 M=pchisq(y,mid,lower.tail=FALSE)
 H=pchisq(y,high,lower.tail=FALSE)
 if(H-L<thres)tmp<-(low+high)/2
 else if((L-a)*(M-a)<0) tmp<-DF(low,mid,y,a,thres)
 else if((M-a)*(H-a)<0) tmp<-DF(mid,high,y,a,thres)
 return (tmp)
}