In [1]:
pkgs <- c("pipeR", "dplyr", "tidyr", "ggplot2", "readr", "readxl")
lapply(pkgs, require, character.only = TRUE)
In [2]:
d <- read_csv("teaching_methods.csv")
head(d)
In [3]:
str(d)
In [4]:
names(d) <- c("id", "name", "sex","math", "stat", "psy_test", "stat_test1", "stat_test2", "teaching_method")
In [5]:
library(rlist)
In [6]:
library(purrr)
- 指導法ごとの統計テスト2の点数
In [7]:
d %>>%
# dplyr::select(stat_test2, teaching_method) %>>%
split(.$teaching_method) %>>%
purrr::map(~ .[["stat_test2"]]) %>>%
rlist::list.cbind()
7.1 一元配置分散分析（対応なし）¶
- 帰無仮説: 4群の母平均は等しい
- 検定統計量$F$: 帰無仮説が正しいとき，自由度$df_1$，$df_2$のF分布に従う
In [8]:
?df
In [9]:
options(repr.plot.width = 4, repr.plot.height = 4)
In [10]:
ggplot(data.frame(x = c(0:5)), aes(x)) +
stat_function(fun = df, args = list(df1 = 3, df2 = 16))
- 有意水準: 5%，片側検定
- 検定統計量の実現値
In [11]:
dd <- d %>>%
select(stat_test2, teaching_method) %>>%
mutate(teaching_method = factor(teaching_method))
str(dd)
oneway.test()
In [12]:
oneway.test(dd$stat_test2 ~ dd$teaching_method, var.equal = TRUE)
aov()
In [13]:
aov(dd$stat_test2 ~ dd$teaching_method)
In [14]:
aov(dd$stat_test2 ~ dd$teaching_method) %>>% summary
anova()
In [15]:
lm(dd$stat_test2 ~ dd$teaching_method) %>>% anova
7.1.2 平方和の分解¶
- 群内平方和
In [16]:
intra_group_sq_sum <- dd %>>%
group_by(teaching_method) %>>%
mutate_each(funs(my_fun = stat_test2 - mean(stat_test2))) %>>%
ungroup %>>%
mutate(stat_test2 = stat_test2^2) %>>%
summarise(res = sum(stat_test2)) %>>%
(.[["res"]])
intra_group_sq_sum
- 群間平方和
In [17]:
all_mean <- dd$stat_test2 %>>% mean
inter_group_sq_sum <- dd %>>%
group_by(teaching_method) %>>%
summarise_each(funs(mean, n())) %>>%
mutate(dev_sq_sum = (mean - all_mean)^2 * n) %>>%
summarise(res = sum(dev_sq_sum)) %>>%
(.[["res"]])
inter_group_sq_sum
- 全体平方和
In [18]:
all_sq_sum <- dd %>>%
mutate(dev_sq_sum = (stat_test2 - all_mean)^2) %>>%
summarise(res = sum(dev_sq_sum)) %>>%
(.[["res"]])
all_sq_sum
分散分析表との対応を見る
In [19]:
aov(dd$stat_test2 ~ dd$teaching_method) %>>% summary
- 全体平方和 = 群間平方和 + 群内平方和
In [20]:
inter_group_sq_sum + intra_group_sq_sum
自由度
- 群間: 群の数 - 1
In [21]:
df_inter <- (unique(dd$teaching_method) %>>% length) - 1
df_inter
- 郡内: sum(各群のデータ数 - 1)
In [22]:
df_intra <- dd %>>%
group_by(teaching_method) %>>%
summarise(n = n()) %>>%
ungroup %>>%
summarise(res = sum(n - 1)) %>>%
(.[["res"]])
df_intra
- 全データ: 全データ数 - 1
In [23]:
df_all <- length(dd$stat_test2) - 1
df_all
全データ = 群間 + 群内
In [24]:
df_inter + df_intra
群間平均平方
In [25]:
inter_mean_sq <- inter_group_sq_sum / df_inter
inter_mean_sq
群内平均平方
In [26]:
intra_mean_sq <- intra_group_sq_sum / df_intra
intra_mean_sq
全体平方和を全体の自由度で割った平均平方は、全データの不偏分散となる
In [27]:
all_sq_sum / df_all
In [28]:
var(d$stat_test2)
- $F$ 値
In [29]:
f_value <- inter_mean_sq / intra_mean_sq
f_value
多重比較（Tukey）¶
- どの群とどの群に差があるかどうか
- データ数 n が等しく，各群の母分散も等しいと仮定して検定統計量$q$を求める
In [30]:
combn(c("A", "B", "C", "D"), 2)
In [31]:
dd %>>% group_by(teaching_method) %>>% summarise(n = n())
In [32]:
n <- dd %>>% group_by(teaching_method) %>>% summarise(n = n()) %>>% (.[["n"]][1])
n
In [33]:
dd %>>%
group_by(teaching_method) %>>%
summarise(gr_mean = mean(stat_test2)) %>>%
(function(x){
ret <- x$gr_mean
names(ret) <- x$teaching_method
ret
})() %>>%
(function(x){
names(x) %>>%
combn(2) %>>%
apply(2, function(y){
a <- y[1]
b <- y[2]
q <- abs(x[a] - x[b]) / sqrt(intra_mean_sq / n)
cat(a, " and ", b, ": ", q, "\n")
})
invisible()
})()
In [34]:
ls.str()
In [35]:
?qtukey
- 検定統計量$q$の実現値
qtukey()
In [36]:
ngroup <- dd$teaching_method %>>% unique %>>% length
ngroup
In [37]:
qtukey(0.95, ngroup, df_intra)
- A-C, A-D は棄却域に入った
TukeyHSD()を用いると，
In [38]:
aov(dd$stat_test2 ~ dd$teaching_method) %>>% TukeyHSD
7.2 一元配置分散分析（対応あり）¶
例題¶
In [39]:
d2 <- read_csv("chap07_2.csv", locale=locale(encoding="CP932")) # this csv file was generated by Excel
d2
In [40]:
names(d2) <- c("student", "algebra", "calculus", "probability")
str(d2)
7.2.1 対応なしとした場合¶
In [41]:
d2 %>>%
gather(subject, point , -student) %>>% head()
In [42]:
d2 %>>%
gather(subject, point , -student) %>>%
{aov(point ~ subject, data = .)} %>>%
summary
7.2.2 対応あり¶
$$ F = \frac{\mbox{条件平方和} / \mbox{条件の自由度}}{\mbox{残差平方和} / \mbox{残差の自由度}}$$In [43]:
d2 %>>% gather(subject, point, -student) %>>%
{aov(point ~ subject + student, data = .)} %>>%
summary
棄却域は，
In [44]:
qf(0.05, 2, 8, lower.tail = FALSE)
- Tukey多重比較
In [45]:
d2 %>>% gather(subject, point, -student) %>>%
{aov(point ~ subject + student, data = .)} %>>%
TukeyHSD()
In [46]:
d2 %>>% gather(subject, point, -student) %>>%
ggplot(aes(x = subject, y = point)) +
geom_boxplot()
7.2.4 平方和の分解と自由度の計算¶
- 対応なし: 全体平方和 $=$ 群間平方和 $+$ 群内平方和
- 対応あり: 全体平方和 $=$ 条件平方和 $+$ 個人差平方和 $+$ 残差平方和
全体平方和
In [47]:
d2.long <- d2 %>>% gather(subject, point, -student)
In [48]:
all_mean <- mean(d2.long$point)
all_mean
In [49]:
all_sq_sum <- d2.long %>>%
mutate(dev_sq_sum = (point - all_mean)^2) %>>%
summarise(res = sum(dev_sq_sum)) %>>%
(.[["res"]])
all_sq_sum
条件平方和
In [50]:
subject_sq_sum <- d2.long %>>%
group_by(subject) %>>%
summarise_each(funs(mean, n()), point) %>>%
mutate(dev_sq_sum = (mean - all_mean)^2 * n) %>>%
summarise(res = sum(dev_sq_sum)) %>>%
(.[["res"]])
subject_sq_sum
個人差平方和
In [51]:
student_sq_sum <- d2.long %>>%
group_by(student) %>>%
mutate_each(funs(my_fun = mean(point) - all_mean), point) %>>%
ungroup %>>%
mutate(point = point^2) %>>%
summarise(res = sum(point)) %>>%
(.[["res"]])
student_sq_sum
残差平方和
In [52]:
residuals_sq_sum <- d2.long %>>%
group_by(student) %>>%
mutate(mean.student = mean(point) - all_mean) %>>%
group_by(subject) %>>%
mutate(mean.subject = mean(point) - all_mean) %>>%
ungroup %>>%
mutate(point = point - all_mean - mean.student - mean.subject) %>>%
mutate(point = point^2) %>>%
summarise(res = sum(point)) %>>%
(.[["res"]])
residuals_sq_sum
分解の確認
In [53]:
all_sq_sum
subject_sq_sum + student_sq_sum + residuals_sq_sum
all_sq_sum == (subject_sq_sum + student_sq_sum + residuals_sq_sum)
identical(all_sq_sum, (subject_sq_sum + student_sq_sum + residuals_sq_sum))
all.equal(all_sq_sum, (subject_sq_sum + student_sq_sum + residuals_sq_sum))
残差の自由度 $=$ 条件自由度 $\times$ 個人差の自由度
$F$値は，
In [54]:
(subject_sq_sum / (3 - 1)) / (residuals_sq_sum / ((3 - 1) * (5 - 1)))
7.3 二元配置分散分析（対応なし）¶
In [55]:
d3 <- read_csv("chap07_3.csv", locale = locale(encoding="CP932"))
d3
In [56]:
d3.long <- d3 %>>%
gather(key, taste) %>>%
separate(key, into = c("water", "condition"), sep = "\\.")
head(d3.long)
7.3.2¶
- 帰無仮説と対立仮設
要因A，要因B，AとBの交互作用効果のそれぞれに関して設定する
In [57]:
aov(taste ~ condition * water, data = d3.long) %>>% summary
In [58]:
aov(taste ~ condition + water + water:condition, data = d3.long) %>>% summary
平方和の分解の確認
In [59]:
d3.res <- aov(taste ~ condition * water, data = d3.long) %>>% summary
str(d3.res)
In [60]:
all_sq_sum <- d3.long %>>%
mutate(dev_sq_sum = (taste - mean(taste))^2) %>>%
summarise(res = sum(dev_sq_sum)) %>>%
(.[["res"]])
all_sq_sum
In [61]:
sum(d3.res[[1]]$`Sum Sq`)
In [62]:
all.equal(all_sq_sum, sum(d3.res[[1]]$`Sum Sq`))
In [63]:
library(Cairo)
In [64]:
CairoFonts(regular = "IPAexGothic")
In [65]:
options(repr.plot.width = 8, repr.plot.height = 4)
In [66]:
Cairo(type="raster")
par(family = "IPAexGothic")
par(cex.axis = 0.8)
layout(matrix(c(1, 2), 1, 2, byrow = TRUE))
with(d3.long, {
interaction.plot(condition, water, taste)
interaction.plot(water, condition, taste)
})
dev.off()
In [67]:
library(grid)
library(gridExtra)
In [68]:
Cairo(type = "raster")
gp1 <- d3.long %>>%
ggplot(aes(x = condition, y = taste, group = water, linetype = water)) +
stat_summary(fun.y=mean, geom="line") +
theme(legend.position = c(.85, .85),
legend.text = element_text(family="IPAexGothic"))
gp2 <- d3.long %>>%
ggplot(aes(x = water, y = taste, group = condition, linetype = condition)) +
stat_summary(fun.y=mean, geom="line") +
theme(legend.position = c(.9, .9),
legend.text = element_text(family="IPAexGothic"))
grid.arrange(gp1, gp2, ncol = 2)
dev.off()
7.3.3 温度の違いのみに注目する（銘柄の違いを無視する）¶
In [69]:
aov(taste ~ condition, data = d3.long) %>>% summary()
- 残差平方和が大きくなる
- 残差平方和は，二元配置の場合の銘柄，銘柄:温度の交互作用，残差 の和と一致
In [70]:
d3.res[[1]]$`Sum Sq`[1]
In [71]:
sum(d3.res[[1]]$`Sum Sq`[2:4])
7.4 二元配置分散分析(2要因とも対応あり)¶
In [72]:
d4 <- d3 %>>% mutate(person = c("村松","川崎","井口","松中","城島"))
d4
In [73]:
d4.long <- d4 %>>%
gather(key, taste, -person) %>>%
separate(key, into = c("water", "condition"), sep = "\\.")
str(d4.long)
In [74]:
aov(taste ~ condition*water +
Error(person + person:condition + person:water + person:condition:water), data = d4.long) %>>%
summary
7.5 二元配置分散分析(1要因のみ対応あり)¶
In [75]:
d5 <- c("A", "B") %>>%
lapply(function(x){
switch(x,
"A" = {person <- c("村松","川崎","井口","松中","城島")},
"B" = {person <- c("斉藤","和田","寺原","杉内","新垣")})
d3 %>>%
select(ends_with(x)) %>>%
mutate(person = person) %>>%
gather(key, taste, -person) %>>%
separate(key, into = c("water", "condition"))
}) %>>%
bind_rows
glimpse(d5)
In [76]:
aov(taste ~ condition * water +
Error(person:condition + person:condition:water), data = d5) %>>% summary
In [77]:
aov(taste ~ condition * water +
Error(person + person:condition + person:water + person:condition:water), data = d5) %>>% summary
In [78]:
aov(taste ~ condition * water +
Error(person), data = d5) %>>% summary
In [79]:
d.ex1 <- read_csv("chap07_ex1.csv", locale = locale(encoding = "cp932"), col_names = FALSE)
d.ex1
In [80]:
d.ex1 %>>%
(function(x){
faculty <- x[[1]]
ret <- x[,-1] %>>% t %>>% data.frame
names(ret) <- faculty
ret
})() %>>%
gather(fuculty, point) %>>%
(~ d.ex1.long) %>>%
{aov(point ~ fuculty, data = .)} %>>%
summary
In [81]:
aov(point ~ fuculty, data = d.ex1.long) %>>% TukeyHSD
（2）¶
In [82]:
d.ex2 <- read_csv("chap07_ex2.csv", locale = locale(encoding = "cp932"), col_names = FALSE, col_types = "cccccccc")
d.ex2
In [83]:
d.ex2 %>>%
(function(x){
header <- enc2native(x[[1]]) # avoid UTF-8 Encoding on Windows
student <- x[1,] %>>% unlist
ret <- x[-1,-1] %>>% t %>>% apply(2, as.numeric) %>>% data.frame()
ret[[student[1]]] <- student[2:length(student)]
ret <- ret[c(4, 1:3)]
names(ret) <- header
ret
})() %>>%
gather(type, point, -学生) %>>%
rename("student" = 学生) %>>%
(~ d.ex2.long) %>>%
{aov(point ~ type + student, data = .)} %>>%
summary
In [84]:
aov(point ~ type + student, data = d.ex2.long) %>>% TukeyHSD
In [85]:
devtools::session_info()