In [25]:
binom.test(c(13, 7), p=1/2)
In [26]:
binom.test(c(65, 35), p=1/2)
12.2 $\chi^2$ 適合度検定¶
$$ X^2 = \sum\limits^n_{i=1}\frac{(O_i - E_i)^2}{E_i} $$- メンデルの法則（分離の法則）
In [28]:
d2 <- matrix(c(152, 39, 53, 6, 140.625, 46.875, 46.875, 15.625), nrow=2, ncol=4, byrow=TRUE)
d2
$X^2$ 値
In [30]:
sum((d2[1,] - d2[2,])**2 / d2[2,])
棄却域
In [31]:
qchisq(0.95, df=3)
In [32]:
chisq.test(d2[1,], p= d2[2,], rescale=TRUE)
有意確率
In [33]:
pchisq(8.9724, df=3, lower.tail = FALSE)
小標本の場合は $X^2$ と $\chi^2$分布 が合わないので，補正をする
- イェーツの連続性補正（コクランの条件）
- 正確$\chi^2検定$
独立性検定ならば
- フィッシャーの正確確率検定
例題: アズキゾウムシ
In [37]:
d3 <- matrix(c(0:7,
0, 0, 60, 101, 84, 48, 7, 0,
9.35, 32.39, 56.20, 65.00, 56.39, 39.14, 22.63, 19.06), nrow=3, ncol=8, byrow=TRUE)
print(d3)
In [40]:
library(dplyr)
In [46]:
d3[c(1:2),] %>% apply(2, function(x){rep(x[1], x[2])}) %>% unlist() %>%
hist(br=seq(-0.5, 9.5, 1), plot=FALSE)
In [48]:
options(repr.plot.width=4, repr.plot.height=4)
In [49]:
d3[c(1:2),] %>% apply(2, function(x){rep(x[1], x[2])}) %>% unlist() %>%
hist(br=seq(-0.5, 9.5, 1), plot=TRUE)
In [62]:
egg <- d3[c(1:2),] %>% apply(2, function(x){rep(x[1], x[2])}) %>% unlist()
plt.d3 <- egg %>% hist(br=seq(-0.5, 9.5, 1), plot=TRUE)
egg.expected <- length(egg) * dpois(c(0:9), mean(egg))
lines(plt.d3$mids, egg.expected)
In [82]:
res.d3 <- rbind(plt.d3$counts, egg.expected)
res.d3
In [83]:
rownames(res.d3) <- c("O", "E")
colnames(res.d3) <- plt.d3$mids
res.d3
In [84]:
res.d3[,8]
In [85]:
res.d3[,c(8:10)]
In [86]:
res.d3[, 8] <- apply(res.d3[,c(8:10)], 1, sum)
res.d3
In [87]:
res.d3 <- res.d3[, c(-9, -10)]
res.d3
In [88]:
chisq.test(res.d3["O",], p=res.d3["E",], rescale=TRUE)
chisq.test() は外部仮説
今回は 内部仮説
自由度が変わるので，有意確率は別に求める
本文では e-21 になっているが
In [90]:
pchisq(105.8636, df=6, lower.tail = FALSE)
12.3 正確$\chi^2$検定: モンテカルロシミュレーション¶
"引数のp=expは期待度数（expected）を意味し..." と書いてあるが，exp は前のページで作った変数（206ページで作ったポアソン乱数）
In [93]:
chisq.test(res.d3["O",], p=res.d3["E",], rescale=TRUE, simulate.p.value=TRUE, B=10000)
12.4 $\chi^2$独立性検定¶
In [94]:
d4 <- matrix(c(25, 15, 10, 30), nrow=2, ncol=2, byrow=TRUE)
print(d4)
In [96]:
rowSums(d4)
In [97]:
colSums(d4)
In [98]:
sum(d4)
In [104]:
d4.expected <- matrix(rep(rowSums(d4) *colSums(d4) / sum(d4), 2), nrow=2, ncol=2, byrow=TRUE)
print(d4.expected)
$X^2$ 値
In [105]:
sum((d4 - d4.expected)**2 / d4.expected)
棄却値
In [109]:
qchisq(0.05, df=1, lower.tail = FALSE)
In [110]:
chisq.test(d4)
補正なし
In [111]:
chisq.test(d4, correct=FALSE)
12.5 $\chi^2$独立性検定 vs. 分割表の対数尤度比検定（G検定）¶
In [112]:
hair <- matrix(c(32, 43, 16, 9, 55, 65, 64, 16), nrow=2, ncol=4, byrow=TRUE)
print(hair)
In [114]:
chisq.test(hair)
$$ G = 2\left[ \sum\limits^m_{i=1}\sum\limits^n_{j=1}O_n\ln\left(\frac{O_{ij}}{E_{ij}}\right) \right] $$
In [115]:
devtools::dev_mode(TRUE, "Chap12")
In [119]:
library("Deducer")
In [120]:
likelihood.test(hair)
In [122]:
devtools::dev_mode(FALSE)
12.6 フィッシャーの正確確率法¶
- 超幾何分布
例題: 淡水性巻貝2種の水流への耐性
In [123]:
shellfish.a <- matrix(c(12, 7, 2, 9), 2, 2, TRUE)
shellfish.a
In [124]:
shellfish.b <- matrix(c(13, 6, 1, 10), 2, 2, TRUE)
In [125]:
shellfish.c <- matrix(c(14, 5, 0, 11), 2, 2, TRUE)
In [126]:
fisher.test(shellfish.a, alternative = "greater")
In [127]:
fisher.test(shellfish.b, alternative = "greater")
In [128]:
fisher.test(shellfish.c, alternative = "greater")
例題: 髪の色
In [129]:
fisher.test(hair)
12.7 オッズ比¶
オッズ比（$\hat{\theta}$）
$$ \hat{\theta} = \frac{p/(1-p)}{q/(1-q)} $$2 $\times$ 2 分割表では，
$$ \hat{\theta} = \frac{a_{11}a_{22}}{a_{12}a_{21}} $$対数オッズ比: 正規分布に収束
$$ \ln\hat{\theta} = \ln a_{11} + \ln a_{22} - \ln a_{12} - \ln a_{21}$$標準誤差は，
$$ SE(\ln \hat{\theta}) = \sqrt{\frac{1}{a_{11}} + \frac{1}{a_{22}} + \frac{1}{a_{12}} + \frac{1}{a_{21}}} $$In [130]:
ls()
In [131]:
d5 <- matrix(c(54, 10325, 25, 51790), 2, 2, TRUE)
d5
オッズ比
In [133]:
d5[1] * d5[4] / (d5[2] * d5[3])
対数オッズ比
In [137]:
log(d5[1]) + log(d5[4]) - log(d5[2]) - log(d5[3])
対数オッズ比の標準誤差
In [138]:
sqrt(1/d5[1] + 1/d5[2] + 1/d5[3] + 1/d5[4])
95%信頼区間
In [158]:
qnorm(c(0.025, 0.975))
In [157]:
2.382737 + qnorm(c(0.025, 0.975)) * 0.242145
In [156]:
exp(2.382737 + qnorm(c(0.025, 0.975)) * 0.242145)
手計算は Agresti (2013)らしい
In [146]:
fisher.test(d5, or=TRUE, conf.level = 0.95)
In [159]:
devtools::session_info()