In [1]:
library(dplyr)
library(tidyr)
library(ggplot2)
library(readr)
library(gridExtra)
library(Cairo)
5.1 2標本の平均値間の“有意な差”とは？¶
In [2]:
data1 <- data_frame(
m = c(143.1, 140.9, 147.2, 139.8, 141.3, 150.7, 149.4, 145.6, 146.5, 148.5, 141.2, 136.5, 145.8, 148.1, 144.3),
f = c(138.7, 142.8, 150.3, 148.4, 141.7, 149.5, 156.5, 144.6, 144.4, 145.7, 148.3, 140.8, 146.2, 149.9, 144.1)
)
head(data1)
In [3]:
data2 <- data_frame(
m = c(142.3, 142.5, 145.7, 143.5, 144.2, 145.1, 145.9, 145.2, 146.8, 145.7, 145.4, 144.6, 144.2, 145.9, 142.1),
f = c(143.5, 144.6, 143.4, 146.6, 145.3, 147.7, 147.2, 147.8, 145.3, 145.7, 147.5, 147.2, 148.8, 147.9, 143.3)
)
head(data2)
In [4]:
data1 %>% summarise_all(.funs = mean)
In [6]:
data2 %>% summarise_all(.funs = mean)
- 「小学6年生だと女子の方が男子よりも身長が高い」のか？
In [19]:
options(repr.plot.height = 4, repr.plot.width = 6)
In [24]:
gp1 <- data1 %>% gather(gender, height) %>%
ggplot(aes(x = gender, y = height)) +
geom_point(alpha = 0.5) + ylim(c(130, 160))
gp2 <- data2 %>% gather(gender, height) %>%
ggplot(aes(x = gender, y = height)) +
geom_point(alpha = 0.5) + ylim(c(130, 160))
grid.arrange(gp1, gp2, ncol=2)
5.2 帰無仮説と対立仮説: “第1種の過誤”と“第2種の過誤”¶
5.3 帰無仮説検定: $t$検定の概念と計算法¶
5.4 $t$ 検定の考え方とRを使った計算¶
In [25]:
surgery <- data_frame(new = c(2, 3, 6, 7, 4, 5, 6, 3),
old = c(5, 7, 5, 8, 9, 7, 7, 6))
surgery
In [28]:
surgery %>% {
t.test(.$new, .$old, var=TRUE)
}
In [32]:
my.t.test <- function(x, y){
n.x <- length(x)
n.y <- length(x)
mean.x <- mean(x)
mean.y <- mean(y)
ss.x <- sum((x - mean.x)^2)
ss.y <- sum((y - mean.y)^2)
sp2 <- (ss.x + ss.y) / (n.x + n.y - 2)
se_mean <- sqrt(sp2 * (n.x + n.y) / (n.x * n.y))
t_val <- (mean.x - mean.y) / se_mean
return(t_val)
}
surgery %>% {
my.t.test(.$new, .$old)
}
In [33]:
qt(0.025, df = 14)
In [34]:
qt(0.975, df = 14)
In [40]:
ggplot(data_frame(x = c(-3, 3)), aes(x)) + stat_function(geom="line", fun = dt, args = list(df = 14)) +
geom_vline(xintercept = -2.144, linetype = 2) +
geom_vline(xintercept = 2.144, linetype = 2)
$t$ 検定のいくつかの事例¶
In [42]:
data1 %>% {
t.test(.$m, .$f, var = TRUE)
}
In [43]:
data2 %>% {
t.test(.$m, .$f, var=TRUE)
}
In [44]:
growth1 <- data_frame(
t = c(5.5, 4.2, 3.7, 5.1, 4.4, 4.3),
c = c(3.9, 4.1, 3.8, 3.2, 4.5, 3.8)
)
growth1 %>% {
t.test(.$t, .$c, var=TRUE)
}
In [45]:
growth1 %>% {
t.test(c(.$t, .$t), c(.$c,.$c), var=TRUE)
}
5.6 検出力を高めるもう一つの方法: 対応のある$t$ 検定¶
In [46]:
foods <- data_frame(
tomato = c(6, 16, 10, 14, 24, 8),
banana = c(1, 10, 5, 15, 20, 3)
)
foods %>% {
t.test(.$tomato, .$banana, var=TRUE)
}
In [49]:
foods %>% mutate(id = c(1:6), d = tomato - banana) %>% select(id, tomato, banana, d)
In [55]:
?summarise_at
In [57]:
foods %>% mutate(id = c(1:6), d = tomato - banana) %>% select(id, tomato, banana, d) %>%
summarise_at(vars(-id), mean)
In [61]:
foods %>% mutate(id = as.character(c(1L:6L))) %>% gather(m, dec, -id) %>%
ggplot(aes(x = id, y = dec, group = m, colour = m)) + geom_point()
In [64]:
?t.test
In [62]:
foods %>% mutate(d = tomato - banana) %>% {
t.test(.$d, mu = 0)
}
- 対応のある$t$検定をするには，
paired = TRUEとすればよい
In [65]:
foods %>% {
t.test(.$tomato, .$banana, paired = TRUE)
}
- 正規性の検定（シャピーロ・ウィルク検定）
- 帰無仮説: 正規分布に従う母集団から標本が抽出された
In [66]:
?shapiro.test
In [69]:
data1$m %>% shapiro.test()
In [70]:
data1$f %>% shapiro.test()
- 等分散性の検定（バートレット検定）
- 帰無仮説: 分散は等しい
var.test()？
In [71]:
?bartlett.test
In [77]:
?as.factor
In [78]:
data1 %>% gather(sex, score) %>% mutate(sex = factor(sex, levels = c("m", "f"))) %>% str()
In [79]:
data1 %>% gather(sex, score) %>% mutate(sex = factor(sex, levels = c("m", "f"))) %>% {
bartlett.test(score ~ sex, data = .)
}
二つの標本が等分散でないときの有意差検定: ウェルチの検定¶
In [82]:
ls()
In [83]:
data3 <- data_frame(
x1 = c(23, 20, 20, 24, 17, 19, 26, 22, 19, 21),
x2 = c(17, 23, 25, 34, 25, 28, 20, 31, 26, 36)
)
In [84]:
options(repr.plot.width = 3, repr.plot.height = 3)
In [86]:
data3 %>% gather(vars, vals) %>% ggplot(aes(x = vars, y = vals)) + geom_boxplot()
In [87]:
data3 %>% summarise_all(.funs = funs(var))
In [89]:
data3 %>% gather(group, score) %>% mutate(group = factor(group, levels = c("x1", "x2"))) %>% {
bartlett.test(score ~ group, data = .)
}
In [90]:
data3 %>% {
t.test(.$x1, .$x2, var = FALSE)
}
In [91]:
growth2 <- data_frame(
t = c(25, 15, 19, 17, 22, 20),
c = c(16, 18, 17, 11, 19, 14)
)
growth2
In [124]:
growth2 %>% gather(type, growth) %>% ggplot(aes(x = type, y = growth)) + geom_boxplot()
In [92]:
growth2$t %>% shapiro.test()
In [93]:
growth2$c %>% shapiro.test()
In [97]:
growth2 %>% gather(type, growth) %>% mutate(type = factor(type)) %>% {
bartlett.test(growth ~ type, data = .)
}
In [98]:
growth2 %>% {
t.test(.$t, .$c, var=TRUE)
}
ex1-2¶
In [126]:
bind_rows(growth2, growth2) %>% gather(type, growth) %>% ggplot(aes(x = type, y = growth)) + geom_boxplot()
In [99]:
growth2 %>% {
t.test(rep(.$t, 2), rep(.$c, 2), var=TRUE)
}
ex.2¶
In [107]:
mental <- list(
pos = c(12, 16, 11, 9, 18, 17, 14, 16, 10, 11),
neg = c(9, 10, 11, 8, 7, 9, 6, 10, 11, 13, 12)
)
mental
In [123]:
names(mental) %>% sapply(function(x){
data_frame(type = x, score = mental[[x]])
}, simplify=FALSE) %>% bind_rows() %>%
ggplot(aes(x = type, y = score)) + geom_boxplot()
In [103]:
mental$pos %>% shapiro.test()
In [104]:
mental$neg %>% shapiro.test()
In [113]:
names(mental) %>% sapply(function(x){
data_frame(type = x, score = mental[[x]])
}, simplify=FALSE) %>% bind_rows() %>% {
bartlett.test(score ~ type, data = .)
}
In [114]:
mental %>% {
t.test(.$pos, .$neg, var=TRUE)
}
ex.3¶
- 問題文が途中で切れているようだ
In [116]:
diet <- data_frame(
id = as.character(c(1L:7L)),
abdomen = c(1.2, 2.0, 1.3, 1.7, 2.0, 1.9, 1.3),
jogging = c(2.3, 2.6, 1.9, 2.8, 1.5, 1.9, 2.2)
)
diet
In [122]:
diet %>% gather(method, loss, -id) %>% ggplot(aes(x = method, y = loss)) + geom_boxplot()
In [118]:
diet$abdomen %>% shapiro.test()
In [119]:
diet$jogging %>% shapiro.test()
In [120]:
diet %>% gather(method, loss, -id) %>% {
bartlett.test(loss ~ method, data = .)
}
In [121]:
diet %>% {
t.test(.$abdomen, .$jogging, var = TRUE, paired = TRUE)
}
問題文にはペアであると書かれているが，解答はペアでない結果のようだ
In [128]:
diet %>% {
t.test(.$abdomen, .$jogging, var = TRUE, paired = FALSE)
}
In [127]:
devtools::session_info()