(引き続きハンバーガー統計学にようこそ!の第3章をRでやってみたものです。)
3.1 チキンの売上は少ないのか
このようなタイプのデータを「とりあえずプロットする」方法として、モザイクプロットというものがある。
sales <- matrix(c(435, 265, 165, 135), nrow = 2) rownames(sales) <- c("wakuwaku", "mogumogu") colnames(sales) <- c("potato", "chicken") mosaicplot(sales, col=TRUE)
モグモグのチキン売上割合が高そうなのがなんとなく見て取れるだろう。
もし差が無かったとした場合の度数(期待度数)は、ちょっとズルいやり方だけど、chisq.test()
の中の$expected
を覗けば分かる。
> (expected <- chisq.test(sales)$expected) potato chicken wakuwaku 420 180 mogumogu 280 120
3.2 カイ二乗値とカイ二乗分布
カイ二乗値を計算してみよう。
> dif <- sales - expected # 期待度数との差を求める > mean(dif) # 平均はゼロ😌 [1] 0 > sum(dif^2) # 自乗してから足せば…🤔 [1] 900 > sum(dif^2 / expected) # 期待度数で割ってから足せば…🤔🤔 [1] 4.464286
ピンポン玉実験の部分もやってみよう。
> x <- 5:9 > (x - 5)^2 / 5 + (10 - x - 5)^2 / 5 [1] 0.0 0.4 1.6 3.6 6.4
自由度1のカイ二乗分布をプロットしてみよう。
chisq <- seq(0, 8, by = 0.01) plot(chisq, dchisq(chisq, df = 1), type = "l", ylim = c(0, 1), xlab = expression(chi^2), ylab = "p")
なおこの例、カイ二乗分布から計算されるp値は近似値なので注意。「オレンジ50:白50の箱から10個ピンポン玉を取り出す」という事象は、箱の外と中のピンポン玉の数を調べることで2x2の分割表で表現でき、従って対応するp値は厳密に計算できる。しかもセルの中には期待度数が5というものが含まれるので、結構ずれるんじゃないかと思われる。
先程描画した自由度1のカイ二乗分布の上に、先程の定義で計算されたカイ二乗値とフィッシャーの正確確率検定で計算されるp値の対応をポイントで示してみよう。
fisher_p <- function(N, ...){ sapply(N, function(n){ fisher.test(matrix(c(50-n, n, 50-(10-n), 10-n), nrow=2), ...)$p.value } ) } chisq_point <- (x - 5)^2 / 5 + (10 - x - 5)^2 / 5 points(chisq_point, fisher_p(5:9))
では次に自由度別のカイ二乗分布を示そう(グラフの描き方は正規分布など好きな関数の曲線を描きたい - Qiitaが参考になった)。
df = c(1, 3, 5) ggplot(data.frame(x=c(0, 8)), aes(x=x)) + mapply(function(d, co){stat_function(fun=dchisq, args=list(df=d), aes_q(color=co))}, df, factor(df)) + scale_color_discrete(name = "df") + labs(x = expression(chi^2), y = "p") + theme_classic() + theme(legend.position = c(0.8, 0.9), legend.direction = "horizontal")
自由度が大きいというのは、分割表で言えばセルの数が多いということに相当する。セルの数が増えれば、それだけ期待値からの相対的なずれの合計も大きくなりやすいし、ずれのばらつきも大きくなるということが見て取れる。
3.3 カイ二乗検定
最初に計算したカイ二乗値と、自由度1のカイ二乗分布の1%点と5%点を比較してみよう。
> sum(dif^2/expected) # そもそもカイ二乗値はいくつだっけ… [1] 4.464286 > qchisq(p=0.95, df = 1) # 左から積分していって0.95になる点 [1] 3.841459 > qchisq(p=0.99, df = 1) [1] 6.634897
5%水準なら帰無仮説を棄却できるし、1%ならできないということが分かった。
が、もちろん最初の分割表をそのままchisq.test()
に放り込んでいい。
> chisq.test(sales) Pearson's Chi-squared test with Yates' continuity correction data: sales X-squared = 4.1716, df = 1, p-value = 0.04111
今回は期待度数が大きいのでカイ二乗分布で近似してもさほど支障はないが、fisher.test()
によってフィッシャーの正確確率検定を実施しても良い。
> fisher.test(sales) Fisher's Exact Test for Count Data data: sales p-value = 0.041 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 1.010939 1.782292 sample estimates: odds ratio 1.342623
ハンバーガーを入れてみる
唐突なタイトル回収だ。
> sales2 <- cbind(sales, hamburger = c(650, 350)) # 唐突なタイトル回収 > sales2 potato chicken hamburger wakuwaku 435 165 650 mogumogu 265 135 350 > mosaicplot(sales2, col = TRUE)
> chisq.test(sales2) Pearson's Chi-squared test data: sales2 X-squared = 9.9048, df = 2, p-value = 0.007067 > fisher.test(sales2) Fisher's Exact Test for Count Data data: sales2 p-value = 0.007509 alternative hypothesis: two.sided
3.9 通過テスト
> d <- matrix(c(24, 8, 18, 18), nrow = 2) > rownames(d) <- c("kokugo", "sansu") > colnames(d) <- c("sakura", "momo") > d sakura momo kokugo 24 18 sansu 8 18 > mosaicplot(d, col=TRUE)
この問題は有意水準1%だと帰無仮説は棄却できず、有意水準5%だとカイ二乗検定かフィッシャーの正確確率検定かで判断が異なる例となっている。そんな些細な差で判断を分けるよりおとなしく追試をしたほうが良いだろうが。
> chisq.test(d) Pearson's Chi-squared test with Yates' continuity correction data: d X-squared = 3.4874, df = 1, p-value = 0.06184 > fisher.test(d) Fisher's Exact Test for Count Data data: d p-value = 0.04643 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.9600903 9.7553314 sample estimates: odds ratio 2.950439
モザイクプロットでは結構差があるように見えるのにp値がそこそこ大きくて、カイ二乗検定とフィッシャーの正確確率検定で結果が異なるのは、要するにサンプルサイズが小さいのだ。
それよりどうもこの問題、ここまでのテキストを参照するに「帰無仮説を採択する」と答えさせようとしているように思えるのだが、帰無仮説は採択できないのでその点には注意してもらいたい(cf. 「帰無仮説を採択」? | Okumura’s Blog)。