コラム「統計学におけるマイ未解決問題」では，χ2分布を用いたピアソンの適合度検定の話を取り上げた．ピアソンの適合度検定では，

　　χ2 ＝Σ（Ｏ−Ｅ）^2 ／Ｅ　　　（Ｏ：観測度数　Ｅ：期待度数）

を作って，このχ2 統計量が自由度＝階級数−１のカイ２乗分布をすることを用いて理論分布と実測値の適合度を検定するわけだが，χ2 統計量がｎ次元楕円：

　　ｘ1^2／ａ^2＋ｘ2^2／ｂ^2＋・・・＋ｘn^2／ｎ^2＝１

と同じ形に表されることを読者諸賢はお気づきであろうか？

　

　χ2 検定法はよく用いられる検定法であるが，期待度数Ｅが小さいとき不正確であるとか，理論的にいろいろな問題点のあることが指摘されている．そのため，当該コラムでは適合度検定の「exact probability」の話を掲載したのだが，その後，京大医学部・薬剤疫学教室の浜田知久馬先生（保健学）より，次のようなメールを頂戴した．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　佐藤先生

　

　京都大学の浜田です．ホームページ拝見しました．

　観測されたパターンの確率より小さくなる全てのパターンの確率を足し合わせて正確なｐ値を計算した場合，ほぼカイ２乗検定の結果に一致するはずです．ただこの計算はかなり大変です．多項分布の周辺和を固定して得られる全ての可能なパターンについて確率を計算しなければなりません．おそらく計算アルゴリズムにミスがあるのではないでしょうか？

　

　　　　　　　観測数　　期待数

　　Ａ型　　　62　　　　200*4/10=80

　　Ｏ型　　　84　　　　200*3/10=60

　　Ｂ型　　　30　　　　200*2/10=40

　　ＡＢ型　　24　　　　200*1/10=20

　　　計　　　200　　　 200

　

(1) DATA SUMMARY

Category　　　No.　　　　　 Prob.

1　　　　　　 62　　　　　　.4

2　　　　　　 84　　　　　　.3

3　　　　　　 30　　　　　　.2

4　　　　　　 24　　　　　　.1

(2) PEARSON CHI-SQUARE TEST

CHI^2= 16.95　　( P<0.01 )　 ( P= .0010516 )

(3) modified CHI-SQUARE TEST

modified CHI^2= 16.0829　　( P<0.01 )　 ( P= .0014559 )

(4) LIKELIHOOD RATIO TEST

LIKELIHOOD RATIO= 16.4112　　( P<0.01 )　 ( P= 1.28573E-03 )

　この例題についての私の計算では，PEARSON CHI-SQUARE TEST(DF=3)のｐ値は，0.0007237，正確なp値は0.0007778とほぼ一致しました．この計算はSASのFREQプロシジャで行いましたが，正確なp値の計算はSAS　Ver8の最新機能です．SASで計算したので間違いないと思います．

　

　また，EXCELのCHIDIST関数の計算結果も0.0007237に一致したので，計算されたピアソン・カイ２乗値 P= .0010516 と乖離がみられるのは，おそらく近似計算の精度の悪さのためだと思います．

　

　ＳＡＳのマニュアルの記述と文献です．

For one-way tables, PROC FREQ computes the exact chi-square goodness-of-fit test by the method of Radlow and Alf (1975).

PROC FREQ generates all possible one-way tables with the observed total sample size and number of categories. For each possible table, PROC FREQ compares its chi-square value with the value for the observed table. If the table's chi-square value is greater than or equal to the observed chi-square, PROC FREQ increments the exact p-value by the probability of that table, which is calculated under the null hypothesis using the multinomial frequency distribution. By default, the null hypothesis states that all categories have equal proportions. If you specify null hypothesis proportions or frequencies using the TESTP= or TESTF= option in the TABLES statement, then PROC FREQ calculates the exact chi-square test based on that null hypothesis.

　

Radlow, R. and Alf, E.F. (1975), "An Alternate Multinomial Assessment of the Accuracy of the Chi-Square Test of Goodness of Fit," Journal of the American Statistical Association, 70, 811 -813.

　

　　　　　　　　　　　　　　　　　　　　浜田

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　浜田先生には，ご多忙中にもかかわらずご回答を頂きましたことに，この場を借りて，御礼申し上げます．

　

　さて，小生の計算アルゴリズムにミスがあることは確実なので，再考してみることにした．早速，前述の文献を取り寄せてみたのだが，論文にある「exact χ2 test」とは，多項分布に基づいて正確に確率を計算したものではなく，ピアソンのχ2 乗検定の欠点を改良したものであることがわかった．中身は今一つ期待はずれだった．

　

　とはいえ，ｎ個のデータをｋ個のカテゴリーに分割する場合の数は，ｋ−１個の仕切り線をつけてやることに等しいから，

　　n+k-1Ｃk-1=(n+k-1)(n+k-2)･･･(n+2)(n+1)/(k-1)!

通りの場合分けが必要になる．上の例題はｎ＝２００，ｋ＝４の場合であるが，ｎやｋが大きいとき，この数は急激に増加するから，周辺和を固定して得られる全ての多項分布の確率（exact probability）を計算することがいかに大変な労力を要するのか，おわかり頂けるであろう．

　

　そこで，exact probabilityを求める本格的なプログラムを作る前のステップとして，その簡便法を考案してみることにした．以下に，そのプロシージャを記す．

　

（１）カテゴリー　　Ｃ1　　　Ｃ2　　・・・　Ｃk　　　計

　　　観測度数　　　ｎ1　　　ｎ2　　・・・　ｎk　　　ｎ

　　　期待度数　　　ｎｐ1　　ｎｐ2　・・・　ｎｐk　　ｎ

のi=1,2,･･･,kに対して，ｎiとｎｐiを比較して，ｎ1＋ｎ2＋・・・＋ｎk＝ｎの条件のもとで，｜ｎi−ｎｐi｜以上離れる確率を多項分布から計算する．

（２）i=1,2,･･･,kの最小のｐ値＝ｐminについて，

　　ｐ＝１−（１−ｐmin）^(k-1)

を計算する．

（３）ｐ値を帰無仮説の有意水準αと比較して，それ以下であったら帰無仮説を棄却する．

　

　(1)について補足しよう．すべてのパターンの確率を足し合わせて正確なｐ値を求めることは，ｎが少し大きくなると大変なわけであるが，ｎiの分布は２項分布Ｂ（ｎ，ｐi）となることを利用して，現在の度数よりもさらに偏ったものが起こる確率は，ｎやｋが大きい場合であっても，コンピュータを使えば容易に計算できる．上の例題では，

　　 Category　　　p-value

　　 1　　　　　　 .0113523

　　2　　　　　　 2.85481E-4

　　 3　　　　　　 .923735

　　 4　　　　　　 .409188

と計算された．

　

　(2)については，ｋカテゴリーのうち，たまたまそのうち１つのｐ値が０．０５より小さいからといって，有意に判定するわけにはいかない．ｋ個のカテゴリーすべてで観測度数と期待度数の差がないとすると，ｋ個のｐ値のうち，ｋ−１個は互いに独立である．そこで，ｋ個のうち少なくとも１つがｐminより小さい値をとる確率ｐは

　　ｐ＝１−（１−ｐmin）^(k-1)

で与えられる．これは，ｐmin，ｋがともに小さければ，近似的に

　　ｐ＝（ｋ−１）ｐmin

となる．

　

　上の例で計算してみると，ｐmin= 2.85481E-4，ｋ＝４より，

　　ｐ＝１−（１−2.85481E-4）^3＝8.56221E-4

すなわち，浜田先生の計算された正確なp値0.0007778に，かなり接近した値が得られたことになる．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　なお，(2)の計算については，最小のｐ値：ｐminではなく，ｐ値の積：

　　Πｐi

を評価する方式も考えられた．これはｐiの幾何平均：

　　k√（Πｐi）

を考えることとと同値だし，さらに，

　　−２ｌｎ（Πｐi）＝−２Σｌｎｐi

を考えることとも同値であるが，ｋ＝２の場合を考慮すると幾何平均を考える方が適当であろう．

　

　この計算は，

　　ｐ＝(.0113523*2.85481E-4*.923735*.409188)^(1/4)＝.0187082

となって，少なくともこの計算例では最小値ｐminを考える方が合理的と見える．

　

　両者の違いは，最小のｐ値を示したカテゴリーを１つ選び出して議論したいのか，分布の偏りの全体にわたる平均を議論したいのかの区別ということができよう．この計算結果は，最小値ｐminを考える方が合理的に評価されたのであって，exact probabilityの簡便検定法としては，やはりｐminを考えるべきであろうと思われた．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【参考文献】

　

Tate, M.W. and Hyer, L.A. (1973), "Inaccuracy of the χ2 Test of Goodnes of Fit When Expected Frequencies Are Small," Journal of the American Statistical Association, 68, 836 -841.