myscratchpad

Holm, Scheffe, TukeyHSD, Dunnet あたりがメジャー。あと、ペリの方法が検定力が高いらしい。
基本的に、群間の独立性を仮定している。すなわち、対応なしの検定である。
ただし、Holm等、有意水準調整型の場合は対応ありでも一応使える

dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), result = c(9,7,8,8,12,11,8,13, 6,5,6,3,6,7,10,9, 10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14))

# 平均と合計
aggregate(dat[2], list(dat[,1]), mean)
aggregate(dat[2], list(dat[,1]), sum)

# lmで分析
lmres <- lm(result~a, dat)
library(car)
Anova(lmres)

## 色々な多重比較

# Holm法、プールされていない標準偏差
## こっちと同じ
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275792703")
val <- dat$result
idv <- dat$a
cbn <- combn(levels(idv), 2)
tresdat <- tobj(t.test(rnorm(10), rnorm(10), paired=F)); tresdat[1,] <- NA # 結果格納用データフレーム
for (i in 1:ncol(cbn)) {
tres <- t.test(val[which(idv==cbn[1,i])], val[which(idv==cbn[2,i])], paired=F)
tob <- tobj(tres)
tresdat[i,] <- tob
}
rownames(tresdat) <- apply(cbn, 2, function(x) paste(x, collapse=","))
tresdat
## Holm法でp値の調整をする
pholms <- p.adjust(tresdat[,"p.vl"], "holm")
(tresdat <- data.frame(tresdat, pholms))
colnames(tresdat) <- c(colnames(tresdat[-ncol(tresdat)]), "Pr(>|t|).Holm")
print.anova(tresdat) ## 有意の星をつける。文字型の変数は01にされる

# Holm法、プールされた標準偏差
## pairwise.t.test
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275849015")
pairwise.t.test(dat$result, dat$a, p.adj="holm")
ptest(dat$result, dat$a)
## 青木先生の関数。感謝
source("http://aoki2.si.gunma-u.ac.jp/R/src/Bonferroni.R", encoding="euc-jp")
Bonferroni(dat$result, dat$a, method="Holm")
print.anova(Bonferroni(dat$result, dat$a, method="Holm")$result2)
x <- Bonferroni(dat$result, dat$a, method="Holm")$result2[,2]
round(p.adjust(x, "holm"), 5)

# scheffe
## 青木先生のコードより。感謝
source("http://aoki2.si.gunma-u.ac.jp/R/src/scheffe.R", encoding="euc-jp")
ns <- tapply(dat$result, dat$a, length)
ms <- tapply(dat$result, dat$a, mean)
us <- tapply(dat$result, dat$a, var)
scheffe(ns, ms, us, 1,2)
scheffe(ns, ms, us, 1,3)
scheffe(ns, ms, us, 1,4)
scheffe(ns, ms, us, 2,3)
scheffe(ns, ms, us, 2,4)
scheffe(ns, ms, us, 3,4)
## こういう関数も見つけた。
http://www.biw.kuleuven.be/vakken/statisticsbyr/ANOVAbyRr/multiplecompJIMRC.htm
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276112608")
aovres <- aov(result~a, dat)
summary(aovres)
nis <- tapply(dat$result, dat$a, length); nis
ms <- tapply(dat$result, dat$a, mean); ms
mse.df <- aovres$df.residual; mse.df
mse <- sum(aovres$residuals^2)/mse.df; mse
scheffeCI(ms, nis, mse.df, mse, conf=.95)
## ちょっと改造してF値とp値を出すようにした
scres <- scheffeCI2(ms, nis, mse.df, mse, conf=.95)
print.anova(scres)



# Tukey
library(multcomp)
summary(glht(lmres,linfct=mcp(a="Tukey"))) # これはなんかおかしい
summary(glht(aov(lmres),linfct=mcp(a="Tukey"))) # 教科書と合うのはこっち
TukeyHSD(aov(lmres))
## 小塩先生のspssでの分析例をやってみた。感謝
cond <- factor(rep(1:3, each=7))
result <- c(4,1,3,2,2,4,3,6,8,5,9,8,7,7,4,3,4,6,5,5,5)
aovres <- aov(result~cond)
TukeyHSD(aovres)
library(multcomp)
summary(glht(aovres, linfct=mcp(cond="Tukey")))
## TukeyHSDとglhtで微妙に結果が違うのはなんでだろーなー

## Dunnet 対照群は一番上の水準 ("a1")
summary(glht(aov(lmres),linfct=mcp(a= "Dunnett")))
## Williams なんだかよくわからない。水準間に順序が仮定できるときに使うらしい。
summary(glht(aov(lmres),linfct=mcp(a= "Williams")))

## そのうちやろう
# 2要因以上の分散分析で交互作用が出たときはプールされた分散をつかい統計量を算出する
# glhtの使い方。特にTukey
# 各群のサンプルサイズが違う場合の多重比較。

参考
http://www.ibaraki-kodomo.com/toukei/posthoc.html
http://www.gen-info.osaka-u.ac.jp/testdocs/tomocom/tazyu.html
http://home.hiroshima-u.ac.jp/keiroh/maeda/statsarekore/posthoc.html