クロス集計表の統計：
カイ二乗分布に関連したいくつかの検定
目次
・余談：クロス集計表に関する統計の個人的に不思議に思っている性質
集計したデータに関する検定
今まで、多様な形のデータに関して、様々な解析を行ってきた。今回は、さらに別の形で取りまとめられたデータに関して、解析する手法を紹介しよう。
コイントスを考える。例えば、1000回コインを投げて、表と裏が出た数を集計し、以下のようなデータを得たとしよう。
果たして、このコインは平等なコインと言えるだろうか。つまり、（表の出る比率）=（裏の出る比率）=1/2といえるだろうか。
あるいは、アンケート調査を行うことを考える。例えば、無作為に男性100人、女性100人を選び、「虫は好きか？」という質問をすることを考えよう。そして得られた回答を集計して、以下のようなデータを得たとしよう。
果たして、「虫が好きか否か」に関して、男女差がないと言えるだろうか。言い換えれば、「男性の虫好きの比率」=「女性の虫好きの比率」と言えるだろうか。
上記のような、各カテゴリごとに集計した表を分割表contingency table、あるいは2標本以上ならクロス集計表cross-tabulation tableと呼ぶ。上記の問題は、一見すればそれぞれ別問題であるが、χ^2（カイ二乗）と呼ばれる統計量を計算するカイ二乗検定によって、いずれも解決することができる。本項では、クロス集計表に関する検定として、カイ二乗検定を紹介し、また、それに関連するいくつかの統計を紹介する。
カイ二乗検定に関するいろんな呼び名
私は、クロス集計表のカイ二乗検定に対して、かなり苦手意識があった。というのも、非常にたくさんの呼び名があるのだ！ 全く同じ計算をしているにもかかわらず、検定をかけるデータの形によって、あたかも別の検定かのような名前がある。以下に、まずまとめておく。
それぞれ、どういうときに呼ばれているか、解説をしよう。
・母比率の検定Hypothesis testing for the population proportion
これは、標本が1つで、多くの場合は二者択一の分割表に関して、それぞれの比率が、帰無仮説と有意に異なるかを検定するときに使われる名前である。ちょうど、上記のコイントスの例がこれにあたる。この時の帰無仮説は、（表の出る比率）=（裏の出る比率）=1/2である。
他にも、サイコロを振って、1とそれ以外の出目を集計した場合も該当する。この時の帰無仮説は、（1が出る比率）= 1/6、あるいは、（1以外が出る比率）= 5/6である。
・適合度の検定Hypothesis testing for goodness of fit
これは、標本が1つで、3つ以上のカテゴリがある分割表に関して、それぞれの比率が、帰無仮説と有意に異なるかを検定するときに使われる名前である。先ほどのサイコロの例を取り上げ直すと、サイコロの各出目をすべて集計した場合にあたる。この時の帰無仮説は、（1が出る比率）= （2が出る比率）=（3が出る比率）=（4が出る比率）=（5が出る比率）=（6が出る比率）= 1/6である。
全ての確率が等しい場合だけでなくても、例えば、各事象が(1/6, 1/2, 1/3)の確率で生じると期待されるような現象に対して、帰無仮説を(1/6, 1/2, 1/3)として検定できる。母比率の検定とは、適合度の検定において、2つのカテゴリしかない場合を指していると言える。
・母比率の差の検定Hypothesis testing for the difference in the population proportions
これは、多くの場合、標本が2つで、二者択一のクロス集計表に関して、それぞれの比率が有意に異なるかを検定するときに使われる名前である。ちょうど、上記のアンケートの例がこれにあたる。この時の帰無仮説は、「男性の虫好きの比率」=「女性の虫好きの比率」である。
また、特に母比率の差の検定というときは、2つの標本のサンプルサイズが大きく異なっているときに使われる傾向があるように見える。だからこそ、「母比率」の差の検定なのだろう。
・独立性の検定Hypothesis testing of independence
これは、標本が2つ以上で、3つ以上のカテゴリがあるクロス集計表に関して、各標本の比率が有意に異なるかを検定するときに使われる名前である。先ほどのアンケートの例を取り上げ直すと、「好き/嫌い」の二者択一ではなく、「好き/どちらかといえば好き/どちらかといえば嫌い/嫌い」の4択にした時に、それぞれを集計した場合にあたる。この時の帰無仮説は、
「男性における好きの比率」=「女性における好きの比率」
かつ「男性におけるどちらかといえば好きの比率」=「女性におけるどちらかといえば好きの比率」
かつ「男性におけるどちらかといえば嫌いの比率」=「女性におけるどちらかといえば嫌いの比率」
かつ「男性における嫌いの比率」=「女性における嫌いの比率」
である。
母比率の差の検定とは、独立性の検定において、2つのカテゴリしかない場合を指していると言える。
もちろん、それぞれがさしている状況は微妙に違うから、呼び分け自体が完全に不要だとまでは思わない。母比率の検定は、母平均の検定（1標本のt検定）と同じような呼び方だし、母比率の差の検定は母平均の差の検定（2標本のt検定など）と同様だ。どっちかというと問題は、一見して何を検定しているのかわからない適合性の検定や独立性の検定だろう。
学べば、言いたいことはわかるのだ。適合性とは、「帰無仮説=期待比率」に対する「適合の度合い=類似性」のことを指している。独立性とは、2つの標本の比率が異なるか=独立なふるまいか、のことを指している。ただ、あまりに名称は多いし、一貫性が無いのだ。別に適合度も独立性も、母比率や母比率の差の検定でよかったのではないか、あるいは適合性、独立性の検定に組み込むか。そもそも、適合性は類似性で、独立性は相違性というような真逆の名前になっているから、適合性も独立性でよかったのではないか、とか。分散分析も、同じ状況っちゃ、そうだが。まあ、これは愚痴だが、せめて今後は、ちょっと長くても「母比率のカイ二乗検定」「適合性のカイ二乗検定」「母比率の差のカイ二乗検定」「独立性のカイ二乗検定」というように、内実が初学者にもわかりやすい形が良いんじゃないかと思っている。
カイ二乗検定の統計量χ^2
カイ二乗検定では、その名の通り、χ^2（カイ二乗）と呼ばれる統計量を計算する。χ（カイ）は、ギリシャ文字で、x（エックス）とは異なる文字であることに注意しよう。母比率および適合度のカイ二乗検定では、カイ二乗は以下のように計算をする。
また、母比率の差および独立性のカイ二乗検定では、カイ二乗は以下のように計算をする。
いずれも共通していることは、実測値（xiやyi）と期待値（XpiやYpi）の差の二乗を、期待値で除したものの総和をとっている点である。これらの総和χ^2が、母比率（適合度）の場合はm-1、母比率の差（独立性）の場合は(m-1)×(標本数-1)の自由度に従うと考えて、検定を行う。
母比率（適合度）のカイ二乗検定では、m = 2のとき、以下のようになる。
確率p1の下、n回中x1回の観測が起こるという、上記のような条件設定は、二項検定の時と同様である。つまり、標本（回数）は、二項分布に従ってデータが生成されると考えられる。統計量を以下のように式変形する。
確率p1の下の二項分布は、nが十分に大きければ、平均np1、分散np1(1 - p1)の正規分布に従う。ゆえに、上記の(x1 - np1)/√(np1(1 - p1))は標準正規分布に従うと考えられる。標準正規分布に従う変数の二乗の和がカイ二乗分布になるので、だからこそ上記の統計量は帰無分布がカイ二乗分布なのだ。そして、標準正規分布の2乗を足した回数が自由度になるので、上記は2乗が1つ、つまり、自由度は1のカイ二乗分布に従う。
そして、上記の式変形は、帰無仮説を明瞭にしている。帰無仮説は「x1の出現率が期待比率p1に一致する」こと、つまり「x1/n = p1」⇔「x1= np1」⇔「x1 - np1 = 0」である。この帰無仮説の形から、まさに、標本比率と母比率の比較の検定なのである。標本比率を標本平均、母比率を母平均と読み替えればこれは1標本のt検定と同じような発想になっていることがわかるだろう。ちなみに、上記の式からわかるように、χ^2の平方根をとって、標準正規分布を使って検定することも可能である。その場合は、z検定とも呼ばれる。
m = 2なら、p1が決まればp2が自動的に決定されるので帰無仮説は上記のとおりである。mが3以上の時、「x1/n = p1」かつ「x2/n = p2」かつ……かつ「xm/n = pm」である。ただし、p1 + p2 + …… + pm = 1であるから、m - 1個までの帰無仮説が成立した時点てm個目の帰無仮説も成立する。また、この時、分布は二項分布ではなく、多項分布に従うと考える。
母比率の差（独立性）のカイ二乗検定では、m = 2のとき、以下のようになる。
帰無仮説は、「x1/X = y1/Y = p1」である。上記の式変形と同様に、この帰無仮説の妥当性を見てみよう。
上記の結果は、2×2のクロス集計表におけるカイ二乗検定の統計量としてよく紹介される形である。さらに式変形を続ける。
さて、標本1については、確率p1の下、X = x1 + x2回中x1回の観測が起こり、標本2については、確率p1の下、Y = y1 + y2回中y1回の観測が起こっている。そして、確率p1の下の二項分布は、XやYが十分に大きければ、それぞれ平均Xp1、Yp1、分散Xp1(1 - p1)、Yp1(1 - p1)の正規分布に従う。これらをそれぞれの期待比率に基づく分布に書き換えると、それぞれ平均p1、p1、分散p1(1 - p1)/X、p1(1 - p1)/Y（定数倍の再生性）の正規分布に従う。つまり、標本1については確率p1の下、X回中x1回の観測が起こったことは、平均p1、分散p1(1 - p1)/Xの正規分布から、x1/Xの値が観測されたとみなせ、同様に標本2については平均p1、分散p1(1 - p1)/Yの正規分布から、y1/Yの値が観測されたとみなせる。
これらの観測された標本比率の差x1/X - y1/Yは、和の再生性から平均0、分散p1(1 - p1)/X + p1(1 - p1)/Yに従う。標本比率の差が従う正規分布の標準偏差で、標本比率の差を割ったものは標準正規分布に従うわけだから、上記の統計量χ^2は標準正規分布に従う変数の2乗の和で表されている。m = 2、標本が2のとき、上記のように和は1つだから、自由度1のカイ二乗分布に従う。
帰無仮説は「x1の出現率とy1の出現率が一致する」こと、つまり「x1/X = y1/Y」⇔「x1/X - y1/Y = 0」である。この帰無仮説の形から、まさに、母比率の差の検定なのである。母比率を母平均と読み替えればこれは2標本のt検定と同じような発想になっていることがわかるだろう。
カイ二乗検定の注意点
上記のように、統計量χ^2の性質を紹介してきたが、その中で、「n（あるいはXやY）が十分に大きいとき」という仮定をしていた。カイ二乗検定が有効に働くためには、観測数がそれなりに大きいときに限られる。サンプルサイズの大きさが、二項分布の正規分布への近似の観点から重要であることは、上記の議論から明白だろう。つまり、カイ二乗検定は、漸近検定の一種であるのだ。あまりに小標本の時には、カイ二乗検定は適切ではないことを覚えておこう。いくつかの基準や考え方がネットで紹介されているので、確認してみるとよいだろう（井口先生のページ、黒木玄先生のTwitter）。一応、小標本のときには、正確検定である二項検定、多項検定、そして、Fisherの正確確率検定（補足）を用いることが推奨されている。ただし、黒木玄先生は、小標本なら自動的に正確検定を行うことに対して警鐘を鳴らしている（検出力の観点から必ずしもいい選択とは言えないため、後述のシミュレーションも参照）。
カイ二乗分布のシミュレーション
では、いくつかの例を通して、カイ二乗分布をシミュレートしてみよう。例えば、サイコロを例にする。確率1/6で1の出目が観察されると期待される状況で、300回ふったとしよう。これをシミュレートしたデータは二項分布の乱数生成を使って以下のようにできる。
------------------------------------------------------
library(plotn)
library(rstatix)
n <- 300
p <- 1/6
x <- rbinom(n = n, size = 1, prob = p)
d <- data.frame(Trial = 1:n, One = x == 1)
table(d$One)
##
## FALSE TRUE
## 248 52
------------------------------------------------------
今回は、300回中、1が52回出たというような感じだ。期待値は50回だからおおむね正しいだろう。では、この「一連の流れ」を10000回行い、そのたびにχ^2を計算し、その分布を確認しよう。分割表の形から、自由度1のカイ二乗分布に従うことが期待される。
------------------------------------------------------
r <- 10000
x <- sapply(rep(n,r), function(x) rbinom(n = x, size = 1, prob = p))
chisq <- apply(t(x), 1, function(d){
tab <- table(d)
((tab[1] - sum(tab)*(1-p))^2)/(sum(tab)*(1-p)*p)
})
histn(chisq, freq = F)#図1の描画
xx <- seq(0, 20, length = 200)
overdraw(points(xx, dchisq(xx, df = 1), type = "l", col = "red"))
------------------------------------------------------
図1 確率1/6、総数300で、シミュレートしたχ^2の分布。赤線は自由度1のカイ二乗分布。
確かに分布はよく一致する。この分布は、観察総数が十分に大きいので、確率の影響をほとんど受けない。
------------------------------------------------------
n <- 300
p <- 1/2
x <- rbinom(n = n, size = 1, prob = p)
d <- data.frame(Trial = 1:n, One = x == 1)
table(d$One)
##
## FALSE TRUE
## 152 148
x <- sapply(rep(n,r), function(x) rbinom(n = x, size = 1, prob = p))
chisq <- apply(t(x), 1, function(d){
tab <- table(d)
((tab[1] - sum(tab)*(1-p))^2)/(sum(tab)*(1-p)*p)
})
histn(chisq, freq = F)#図2の描画
xx <- seq(0, 20, length = 200)
overdraw(points(xx, dchisq(xx, df = 1), type = "l", col = "red"))
------------------------------------------------------
図2 確率1/2、総数300で、シミュレートしたχ^2の分布。赤線は自由度1のカイ二乗分布。
やはりよく一致する。では今度は、すべての出目をそれぞれちゃんと数えることを考える。この状況は、多項分布の乱数生成を使って以下のようにシミュレートできる。
------------------------------------------------------
n <- 300
p <- rep(1/6, 6)
x <- rmultinom(n = n, size = 1, prob = p)
d <- data.frame(Trial = 1:n, Dice = apply(t(x), 1, FUN = function(x) which.max(x)))
table(d$Dice)
##
## 1 2 3 4 5 6
## 53 42 57 45 56 47
------------------------------------------------------
イカサマのないサイコロを仮定しているから、確かにすべての出目は期待値通り50前後である。では、ここからχ^2を計算しよう。この統計量は、自由度6 - 1 = 5のカイ二乗分布に従うことが期待される。
------------------------------------------------------
x <- sapply(rep(n,r), function(x) {
apply(t(rmultinom(n = n, size = 1, prob = p)), 1, function(x) which.max(x))
})
chisq <- apply(t(x), 1, function(d){
tab <- table(d)
sum(((tab - sum(tab*p))^2)/(sum(tab)*p))
})
histn(chisq, freq = F)#図3の描画
xx <- seq(0, 20, length = 200)
overdraw(points(xx, dchisq(xx, df = 5), type = "l", col = "red"))
------------------------------------------------------
図3 各確率1/6、総数300で、シミュレートしたχ^2の分布。赤線は自由度5のカイ二乗分布。
ちゃんと自由度に応じた変化がみられる。
次はクロス集計表の例を挙げよう。以下のように、2つの標本は同じ確率1/6に従う二項分布から生成されたことを考える。状況としては、2つのサイコロがあって、それぞれ150回ずつふり、1の出目とそれ以外を分けて集計し、2つのサイコロは1の出目の出現率が等しいかを検定している。2つのサイコロは、ともに25前後1が出現することが期待される。
------------------------------------------------------
n <- 300
p <- 1/6
x <- rbinom(n = n, size = 1, prob = p)
d <- data.frame(Trial = 1:n, Type = c(rep("x",n/2), rep("y",n/2)), Yes = x == 1)
table(d$Type, d$Yes)
##
## FALSE TRUE
## x 123 27
## y 132 18
------------------------------------------------------
このような状況の時、同様にχ^2を計算して、その分布を確認する。このときの自由度は、(2 - 1)(2 - 1) = 1である。
------------------------------------------------------
x <- sapply(rep(n,r), function(x) rbinom(n = x, size = 1, prob = p))
chisq <- apply(t(x), 1, function(d){
nt <- length(d)
dt <- data.frame(Type = c(rep("x",nt/2), rep("y",nt/2)), Yes = d == 1)
tab <- table(dt)
xt <- sum(tab[1,])
yt <- sum(tab[2,])
Ft <- sum(tab[,1])
Tt <- sum(tab[,2])
etab <- matrix(c(Ft*xt/n, Tt*xt/n, Ft*yt/n, Tt*yt/n), nrow = 2, byrow = T)
sum(((tab - etab)^2)/etab)
})
histn(chisq, freq = F)#図4の描画
xx <- seq(0, 20, length = 200)
overdraw(points(xx, dchisq(xx, df = 1), type = "l", col = "red"))
------------------------------------------------------
図4 2つの標本はともに確率1/6、総数150の二項分布から生成されたとき、シミュレートしたχ^2の分布。赤線は自由度1のカイ二乗分布。
確かに、よく一致する。
では、4つの標本を考え、すべてで1/6の確率で確認できる事象があるとする。この時のクロス集計表は以下の通り。今度は、4つのサイコロというわけだ。1の期待値は、12.5である。
------------------------------------------------------
n <- 300
p <- 1/6
x <- rbinom(n = n, size = 1, prob = p)
d <- data.frame(Trial = 1:n, Type = c(rep("x",n/4), rep("y",n/4),
rep("z",n/4), rep("w",n/4)), Yes = x == 1)
table(d$Type, d$Yes)
##
## FALSE TRUE
## w 62 13
## x 60 15
## y 64 11
## z 62 13
------------------------------------------------------
このような状況の時、同様にχ^2を計算して、その分布を確認する。このときの自由度は、(2 - 1)(4 - 1) = 3である。
------------------------------------------------------
x <- sapply(rep(n,r), function(x) rbinom(n = x, size = 1, prob = p))
chisq <- apply(t(x), 1, function(d){
nt <- length(d)
dt <- data.frame(Type = c(rep("x",nt/4), rep("y",nt/4),
rep("z",nt/4), rep("w",nt/4)), Yes = d == 1)
tab <- table(dt)
xt <- sum(tab[1,])
yt <- sum(tab[2,])
zt <- sum(tab[3,])
wt <- sum(tab[4,])
Ft <- sum(tab[,1])
Tt <- sum(tab[,2])
etab <- matrix(c(Ft*xt/n, Tt*xt/n, Ft*yt/n, Tt*yt/n,
Ft*zt/n, Tt*zt/n, Ft*wt/n, Tt*wt/n), nrow = 4, byrow = T)
sum(((tab - etab)^2)/etab)
})
histn(chisq, freq = F)#図5の描画
xx <- seq(0, 20, length = 200)
overdraw(points(xx, dchisq(xx, df = 3), type = "l", col = "red"))
------------------------------------------------------
図5 4つの標本はすべて確率1/6、総数75の二項分布から生成されたとき、シミュレートしたχ^2の分布。赤線は自由度3のカイ二乗分布。
やはり、よく一致する。
では、最後に小標本の時、いかにカイ二乗分布への近似が不適切かを確認しよう。サンプルサイズを10としよう。
------------------------------------------------------
n <- 10
p <- 1/2
x <- rbinom(n = n, size = 1, prob = p)
d <- data.frame(Trial = 1:n, One = x == 1)
table(d$One)
##
## FALSE TRUE
## 6 4
x <- sapply(rep(n,r), function(x) rbinom(n = x, size = 1, prob = p))
chisq <- apply(t(x), 1, function(d){
tab <- table(d)
((tab[1] - sum(tab)*(1-p))^2)/(sum(tab)*(1-p)*p)
})
histn(chisq, freq = F)#図6の描画
xx <- seq(0, 20, length = 200)
overdraw(points(xx, dchisq(xx, df = 1), type = "l", col = "red"))
------------------------------------------------------
図6 確率1/2、総数10で、シミュレートしたχ^2の分布。赤線は自由度1のカイ二乗分布。
明らかに分布が一致していない。また、予測されるカイ二乗分布よりも、実際の分布はより大きな値のほうへ裾が重い分布となっている。これが後述するような、小標本において正確検定よりも、カイ二乗検定の危険率が高い原因となっている。
別の状況も確認しておく。
------------------------------------------------------
n <- 20
p <- 1/6
x <- rbinom(n = n, size = 1, prob = p)
d <- data.frame(Trial = 1:n, One = x == 1)
table(d$One)
##
## FALSE TRUE
## 16 4
x <- sapply(rep(n,r), function(x) rbinom(n = x, size = 1, prob = p))
chisq <- apply(t(x), 1, function(d){
tab <- table(d)
((tab[1] - sum(tab)*(1-p))^2)/(sum(tab)*(1-p)*p)
})
histn(chisq, freq = F)#図7の描画
xx <- seq(0, 20, length = 200)
overdraw(points(xx, dchisq(xx, df = 1), type = "l", col = "red"))
------------------------------------------------------
図7 確率1/6、総数20で、シミュレートしたχ^2の分布。赤線は自由度1のカイ二乗分布。
やはり、やや裾が重い。
カイ二乗検定のシミュレーション
次に、様々なサンプルサイズ（5～1000）で、カイ二乗検定の危険率を、正確検定である二項検定の危険率の比較してみよう。上記のような、クロス集計表が与えられた時、chisq.test関数を使うことでカイ二乗検定を行える。二項検定はbinom.test関数で行える。なお、以下のコマンドを行うと、めちゃくちゃエラーが出てくる。これは、小標本の時、カイ二乗検定を行うと、不適切だよ～、とRが教えてくれるためである。
------------------------------------------------------
r <- 10000
s <- c(5,10,20,30,40,50,75,100,150,200,250,300,400,500,750,1000)
chip <- NULL
binp <- NULL
for(i in s){
x <- sapply(rep(i,r), function(x) rbinom(n = x, size = 1, prob = 1/6))
chip_t <- apply(t(x), 1, function(d){
y <- c(sum(d == 0), i - sum(d == 0))
try(chisq.test(y, p = c(5/6, 1/6))$p.value, silent = T)}
)
binp_t <- apply(t(x), 1, function(d) {
y <- c(sum(d == 0), i - sum(d == 0))
try(binom.test(y, p = 5/6)$p.value, silent = T)}
)
chip <- c(chip, sum(chip_t < 0.05))
binp <- c(binp, sum(binp_t < 0.05))
}
ps <- data.frame(size = s, chi = chip, bin = binp)
plotn(ps[,1], ps[,2:3]/r, mode = "m", xlab = "サンプルサイズ", ylab = "危険率",
legend = T, line = T, leg.sp = 4)#図8の描画
overdraw(abline(h = 0.05, col = "red"))
------------------------------------------------------
図8 確率1/6、総数5～1000で、カイ二乗検定と二項検定を行ったときの危険率
全体傾向として、カイ二乗検定は危険率がやや高く、二項検定は危険率がやや低い。カイ二乗検定の危険率が小標本の時高くなりうることは、上記の議論から予測がつくが、一方で、二項検定はかなり保守的であるということがわかる。検出力を高めたいならカイ二乗、危険率を抑えたいなら二項検定という感じで、自分の研究にあった検定を吟味したほうがよさそうだ。サンプルサイズが大きくなれば、本当にどちらがいいとも言えないので、悩むくらいならちゃんとサンプルサイズをそろえようということだろうか。あと、なぜか、サンプルサイズ30のところで変な挙動になっているが、原因はわからない。
Rで行うカイ二乗検定
では、具体的データをもとにカイ二乗検定を行ってみよう。今、サイコロを振って、1か、それ以外かを集計したとする。92回1が出現し、208回それ以外が出現した、以下のようなデータが得られたとしよう。
------------------------------------------------------
tab <- c(208, 92)
------------------------------------------------------
さて、このサイコロは、1の出現に関してイカサマダイスといえるだろうか？ それを確認できるのがカイ二乗検定で、chisq.test関数で以下のように確認できる。このとき、サイコロを仮定しているから、1の期待確率は1/6、それ以外は5/6のはずだ。
------------------------------------------------------
chisq.test(tab, p = c(5/6, 1/6))
##
## Chi-squared test for given probabilities
##
## data: tab
## X-squared = 42.336, df = 1, p-value = 7.686e-11
------------------------------------------------------
検定の結果、p-valueを確認すると、p < 0.05だから、有意水準5%で期待確率とは有意な違いがあると言える。つまり、このサイコロは1の出現率に関してイカサマダイスだ。
それぞれの値を手計算しよう。まず、X-squaredはχ^2のことで、以下の通り。
------------------------------------------------------
chisq <- ((tab[1] - sum(tab)*5/6)^2)/(sum(tab)*5/6*1/6)
chisq
## [1] 42.336
------------------------------------------------------
dfは自由度で、1。そして、自由度1のカイ二乗分布に従うと考えているので、p-valueは以下の通り。
------------------------------------------------------
pchisq(chisq, df = 1, lower.tail = F)
## [1] 7.686459e-11
------------------------------------------------------
ちなみに、カイ二乗検定ではなく、母比率の検定という名前としてのprop.test関数もあって、以下のように同じ結果を得られる。
------------------------------------------------------
prop.test(as.table(tab), p = 5/6, correct = F)
##
## 1-sample proportions test without continuity correction
##
## data: as.table(tab), null probability 5/6
## X-squared = 42.336, df = 1, p-value = 7.686e-11
## alternative hypothesis: true p is not equal to 0.8333333
## 95 percent confidence interval:
## 0.6389838 0.7427942
## sample estimates:
## p
## 0.6933333
------------------------------------------------------
同様の値がいくつか見られるだろう。alternative hypothesisでは帰無仮説で母比率が5/6であることを仮定したことがかかれている。そして、その下に95%信頼区間、そしてsample estimateは208/300を表している。
そして、χ^2の平方根は標準正規分布に従うから、標準正規分布を使ってp値を計算できて以下の通り。
------------------------------------------------------
pnorm(sqrt(chisq), lower.tail = F)*2
## [1] 7.686459e-11
------------------------------------------------------
以下は小標本の例を出す。カイ二乗検定だけでなく、二項検定も行った。
------------------------------------------------------
tab <- c(6, 4)
chisq.test(tab, p = c(5/6, 1/6))
## Warning in chisq.test(tab, p = c(5/6, 1/6)): Chi-squared approximation may be
## incorrect
##
## Chi-squared test for given probabilities
##
## data: tab
## X-squared = 3.92, df = 1, p-value = 0.04771
chisq <- ((tab[1] - sum(tab)*5/6)^2)/(sum(tab)*5/6*1/6)
chisq
## [1] 3.92
pchisq(chisq, df = 1, lower.tail = F)
## [1] 0.04771488
prop.test(as.table(tab), p = 5/6, correct = F)
## Warning in prop.test(as.table(tab), p = 5/6, correct = F): Chi-squared
## approximation may be incorrect
##
## 1-sample proportions test without continuity correction
##
## data: as.table(tab), null probability 5/6
## X-squared = 3.92, df = 1, p-value = 0.04771
## alternative hypothesis: true p is not equal to 0.8333333
## 95 percent confidence interval:
## 0.3126738 0.8318197
## sample estimates:
## p
## 0.6
pnorm(sqrt(chisq), lower.tail = F)*2
## [1] 0.04771488
binom.test(tab, p = 5/6)
##
## Exact binomial test
##
## data: tab
## number of successes = 6, number of trials = 10, p-value = 0.06973
## alternative hypothesis: true probability of success is not equal to 0.8333333
## 95 percent confidence interval:
## 0.2623781 0.8784477
## sample estimates:
## probability of success
## 0.6
------------------------------------------------------
カイ二乗検定では有意になったが、二項検定では有意にならなかった。それぞれの検定の特徴が出ていると言えるだろう。
次に、1標本で3カテゴリ以上の場合だ。例えば、サイコロを300回ふって、すべての出目を数えたとしよう。1～6まで、それぞれ35、46、48、45、48、78回出現した。この時、カイ二乗検定の結果も示す。期待確率はサイコロだから、すべて1/6だ。
------------------------------------------------------
tab <- c(35, 46, 48, 45, 48, 78)
chisq.test(tab, p = rep(1/6, 6))
##
## Chi-squared test for given probabilities
##
## data: tab
## X-squared = 21.16, df = 5, p-value = 0.0007556
chisq <- sum(((tab - sum(tab*rep(1/6, 6)))^2)/(sum(tab)*rep(1/6, 6)))
chisq
## [1] 21.16
pchisq(chisq, df = 5, lower.tail = F)
## [1] 0.0007555674
------------------------------------------------------
すると、5%水準で有意差があり、期待確率とは異なっていることがわかる。ゆえにこのサイコロはイカサマダイスだ。ただし、どこが期待される確率と異なるかまではわからない。それを確認する場合は、多重比較が必要になる。
小標本の場合は、以下の通り。二項検定の拡張版、多項検定をrstatixライブラリの、multinom_test関数を使って行った。
------------------------------------------------------
tab <- c(2, 1, 3, 2, 3, 9)
chisq.test(tab, p = rep(1/6, 6))
## Warning in chisq.test(tab, p = rep(1/6, 6)): Chi-squared approximation may be
## incorrect
##
## Chi-squared test for given probabilities
##
## data: tab
## X-squared = 12.4, df = 5, p-value = 0.0297
chisq <- sum(((tab - sum(tab*rep(1/6, 6)))^2)/(sum(tab)*rep(1/6, 6)))
chisq
## [1] 12.4
pchisq(chisq, df = 5, lower.tail = F)
## [1] 0.02969946
multinom_test(tab, p = rep(1/6, 6), detailed = T)
## # A tibble: 1 x 3
## p p.signif method
## * <dbl> <chr> <chr>
## 1 0.0754 ns Exact multinomial test
------------------------------------------------------
カイ二乗検定では有意になったが、多項検定では有意にならなかった。それぞれの検定の特徴が出ていると言えるだろう。
次は、2×2のクロス集計表を考える。例えば、架空のアンケート調査で男性300人、女性600人に「甘いもの」が好きか尋ねたとして、以下のような結果が得られたとしよう。男性では嫌いに250人、好きに50人、女性では嫌いに399人、好きに201人投票したような状況だ。このとき、性別は「甘いもの」が好きかどうかに影響を与えているだろうか。より正確には、男性での甘いもの好き率は、女性と同じだろうか。
------------------------------------------------------
tab <- matrix(c(250,50,399,201), nrow = 2, byrow = T)
chisq.test(tab, correct = F)
##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 28.18, df = 1, p-value = 1.106e-07
n <- sum(tab)
xt <- sum(tab[1,])
yt <- sum(tab[2,])
Ft <- sum(tab[,1])
Tt <- sum(tab[,2])
etab <- matrix(c(Ft*xt/n, Tt*xt/n, Ft*yt/n, Tt*yt/n), nrow = 2, byrow = T)
etab #期待値
## [,1] [,2]
## [1,] 216.3333 83.66667
## [2,] 432.6667 167.33333
chisq <- sum(((tab - etab)^2)/etab)
chisq
## [1] 28.17973
pchisq(chisq, df = 1, lower.tail = F)
## [1] 1.105563e-07
pnorm(sqrt(chisq), lower.tail = F)*2
## [1] 1.105563e-07
------------------------------------------------------
結果は、5%水準で有意差ありで、男性の甘いもの好き率は、女性の甘いもの好き率と等しくないと言えそうだ。
最後に小標本を場合を考える。正確検定のFisherの正確確率検定をfisher.test関数で行った。
------------------------------------------------------
tab <- matrix(c(3,7,8,2), nrow = 2, byrow = T)
chisq.test(tab, correct = F)
## Warning in chisq.test(tab, correct = F): Chi-squared approximation may be
## incorrect
##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 5.0505, df = 1, p-value = 0.02462
n <- sum(tab)
xt <- sum(tab[1,])
yt <- sum(tab[2,])
Ft <- sum(tab[,1])
Tt <- sum(tab[,2])
etab <- matrix(c(Ft*xt/n, Tt*xt/n, Ft*yt/n, Tt*yt/n), nrow = 2, byrow = T)
etab
## [,1] [,2]
## [1,] 5.5 4.5
## [2,] 5.5 4.5
chisq <- sum(((tab - etab)^2)/etab)
chisq
## [1] 5.050505
pchisq(chisq, df = 1, lower.tail = F)
## [1] 0.02461876
pnorm(sqrt(chisq), lower.tail = F)*2
## [1] 0.02461876
fisher.test(tab)
##
## Fisher's Exact Test for Count Data
##
## data: tab
## p-value = 0.06978
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.007870555 1.133635839
## sample estimates:
## odds ratio
## 0.1226533
------------------------------------------------------
カイ二乗検定では有意になったが、Fisherの正確確率検定では有意にならなかった。それぞれの検定の特徴が出ていると言えるだろう。
補足：Fisherの正確確率検定
任意のクロス集計表において、カイ二乗のような漸近検定ではなく、正確検定を行うことができる。それが、Fisherの正確確率検定Fisher's exact testである。二項検定と同じく、保守的な方法であるが、観測値の少ないクロス集計表では、やはり重要な解析手法となる。
解析の名前の通り、この方法は、手元の標本から考えられるすべてのクロス集計表が得られる確率を具体的に求める。そして、手元の標本が得られる確率、およびそれよりも珍しい標本が得られる確率の総和をp値とするのである。まずに、任意のクロス集計表が得られる確率は以下で与えられる。!は階乗を表す。
m = 2、n = 2の場合は以下のようになる。
解析の名前の通り、この方法は、手元の標本から考えられるすべてのクロス集計表が得られる確率を具体的に求める。そして、手元の標本が得られる確率、およびそれよりも珍しい標本が得られる確率の総和をp値とするのである。まずに、任意のクロス集計表が得られる確率は以下で与えられる。!は階乗を表す。
そして、以下の手順でp値を求める。
では、具体的に解析してみよう。まずは、上記の例を再掲する。
------------------------------------------------------
tab <- matrix(c(3,7,8,2), nrow = 2, byrow = T)
tab
## [,1] [,2]
## [1,] 3 7
## [2,] 8 2
------------------------------------------------------
fisher.test関数を使って、Fisherの正確確率検定を行うことができる。
------------------------------------------------------
fisher.test(tab)
##
## Fisher's Exact Test for Count Data
##
## data: tab
## p-value = 0.06978
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.007870555 1.133635839
## sample estimates:
## odds ratio
## 0.1226533
------------------------------------------------------
そのときのp値は以下のように求める。その際、各テーブルが得られる確率も図示しよう。
------------------------------------------------------
n <- sum(tab)
r1 <- sum(tab[1,])
r2 <- sum(tab[2,])
c1 <- sum(tab[,1])
c2 <- sum(tab[,2])
P <- NULL
for(i in -tab[1,2]:tab[1,1]){
P <- c(P, choose(n = c1, k = tab[1,1] - i) * choose(n = c2, k = tab[1,2] + i)/choose(n = n, k = r1))
}
P0 <- choose(n = c1, k = tab[1,1]) * choose(n = c2, k = tab[1,2])/choose(n = n, k = r1)
col <- rep("pink", length = length(-tab[1,2]:tab[1,1]))
col[P > P0] <- "grey"
barplotn(P, col.fill = col)#図9の描画
overdraw(abline(h = P0))
------------------------------------------------------
図9 確率の分布。水平線は手元の標本が得られる確率。ピンク色部分は手元の標本が得られる確率よりも低い部分。
もう一例、解析例を提示しておく。
------------------------------------------------------
tab <- matrix(c(1,4,4,6), nrow = 2, byrow = T)
tab
## [,1] [,2]
## [1,] 1 4
## [2,] 4 6
fisher.test(tab)
##
## Fisher's Exact Test for Count Data
##
## data: tab
## p-value = 0.6004
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.006118131 6.509053681
## sample estimates:
## odds ratio
## 0.3989622
n <- sum(tab)
r1 <- sum(tab[1,])
r2 <- sum(tab[2,])
c1 <- sum(tab[,1])
c2 <- sum(tab[,2])
P <- NULL
for(i in -tab[1,2]:tab[1,1]){
P <- c(P, choose(n = c1, k = tab[1,1] - i) * choose(n = c2, k = tab[1,2] + i)/choose(n = n, k = r1))
}
P0 <- choose(n = c1, k = tab[1,1]) * choose(n = c2, k = tab[1,2])/choose(n = n, k = r1)
col <- rep("pink", length = length(-tab[1,2]:tab[1,1]))
col[P > P0] <- "grey"
barplotn(P, col.fill = col)#図10の描画
overdraw(abline(h = P0))
------------------------------------------------------
図10 確率の分布。水平線は手元の標本が得られる確率。ピンク色部分は手元の標本が得られる確率よりも低い部分。
余談：クロス集計表に関する統計の個人的に不思議に思っている性質
私がこの解析を知って驚いたのは、クロス集計表は基本的に、行方向と列方向では異なる性質であるのに、最終的な確率計算では、行方向と列方向のデータは平等に扱われる点だ。
上記の例でいえば、男性と女性の調査人数は制御可能なデータだ。一方で、好きか嫌いかという回答は、調査によって制御はできない。この時、とてもじゃないが、男性、女性人数と、好き、嫌いと回答した人数を平等に扱えるわけないと思うのだが、カイ二乗検定では区別していない点が、少し興味深い。