myscratchpad

1要因反復測定の分散分析をaov, nlmeパッケージのlme関数、lme4パッケージのlmer関数でやってみた。
lmeは線型混合モデル、lmerは混合効果モデルでの分析

dat <- read.table("http://personality-project.org/r/datasets/R.appendix3.data", header=T)
colnames(dat) <- c("no", "sbj", "x", "y") # 列名を短くする
aggregate(dat$y, list(dat$x), mean, na.rm=TRUE) # 平均
aggregate(dat$y, list(dat$x), sd, na.rm=TRUE) # 標準偏差
length(levels(dat$x)) # 水準数を調べる
levels(dat$x) # 水準の中身を調べる

# aovでやってみた
aovres <- aov(y ~ x + Error(sbj/x), dat)
## aov(y ~ x + sbj + Error(x:sbj), dat) # これでも同じ。
summary(aovres)

# lmeでやってみた
require(nlme)
lmeres <- lme(y ~ x, random = ~1|sbj/x, data=dat)
summary(lmeres)
anova(lmeres)

# lmerでやってみた
require(lme4)
lmerres <- lmer(y ~ x + (1|sbj), data=dat)
summary(lmerres)
anova(lmerres)
## lmerres2 <- lmer(y ~ x + (x|sbj), data=dat) # ランダム切片/傾きモデル

当然だが、全て同じ結果を返す。lme, lmerについてはanova(...) でF値がわかる。


# 多重比較
#holm法
pairwise.t.test(dat$y, dat$x, p.adjust.method="holm", paired=TRUE)
# Tukey法。対応ありのデータでTukey法を使うのはだめ、と誰かが言ってた気がするが
require(multcomp)
summary(glht(aovres,linfct=mcp(x="Tukey"))) # aov関数の結果は受け取れず、エラーを返す
summary(glht(lmeres,linfct=mcp(x="Tukey")))
summary(glht(lmerres,linfct=mcp(x="Tukey")))


参考: http://gribblelab.org/2009/03/09/repeated-measures-anova-using-r/
愚痴: 分散分析というのはマイナーな手法なのか、解説がない。aovを使ったもの (タイプ1平方和) は多いが、タイプ3では計算できない。carパッケージのAnova関数を使えばタイプ3で実行できるが、今度は誤差項を指定できない。混合デザインに使えなくもないが。事後検定についても書いてあるところはないし、aov関数のように汎用性が高くないし、分析結果に変量効果 (random effect) が入らない (たぶん計算はしてる) 。
ほとんどのページがしれっとaovでのやり方を書いてるけどそれありなの？