Rによるハンバーガー統計学 #3

(引き続きハンバーガー統計学にようこそ!の第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) 

f:id:Rion778:20170802220451p:plain

モグモグのチキン売上割合が高そうなのがなんとなく見て取れるだろう。

もし差が無かったとした場合の度数(期待度数)は、ちょっとズルいやり方だけど、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")

f:id:Rion778:20170802221255p:plain

なおこの例、カイ二乗分布から計算される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))

f:id:Rion778:20170802231028p:plain

では次に自由度別のカイ二乗分布を示そう(グラフの描き方は正規分布など好きな関数の曲線を描きたい - 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")

f:id:Rion778:20170802231300p:plain

自由度が大きいというのは、分割表で言えばセルの数が多いということに相当する。セルの数が増えれば、それだけ期待値からの相対的なずれの合計も大きくなりやすいし、ずれのばらつきも大きくなるということが見て取れる。

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)

f:id:Rion778:20170802232123p:plain

> 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)

f:id:Rion778:20170802232937p:plain

この問題は有意水準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)。