myscratchpad

by関数とAnovaのerror引数を使ってプールした誤差項による検定を行う
こちらのページを参考にした。感謝。卒論でやったという。最近の卒論はレベルが高いのう
データはテクニカルブックp95
2要因とも対応がない場合

# データ入力
dat <- data.frame(
s=factor(paste("s", 1:40, sep="")),
a=factor(c(rep("a1", 20), rep("a2", 20))),
b=factor(rep(c(rep("b1",5), rep("b2",5), rep("b3",5), rep("b4",5)),2)),
result = c(3,3,1,3,5, 4,3,4,5,7, 6,6,6,4,8, 5,7,8,7,9, 3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,5, 2,3,3,4,6))

# 合計を出す
aggregate(dat[4], list(dat[,3], dat[,2]), sum)

# lmとcarのAnovaで分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(result~a*b, data=dat)
library(car)
Anovares <- Anova(lmres, type=3)

# a要因についてプールした誤差項で検定
by(dat, dat$a, function(x) Anova(lm(result~b, data=x), type=2, error=lmres))
# b要因についてプールした誤差項で検定
by(dat, dat$b, function(x) Anova(lm(result~a, data=x), type=2, error=lmres))
## セルサイズが等しいため、テクニカルブックp101どおりの結果になる
## by関数がこんなに便利とは、Anovaで誤差項が指定できるとは知らなかった

## 問題
上記の分析で、type=3とするとプールした誤差項を使わず、そのまんま水準別で検定したのと同じ結果になる というか、linear.hypothesisでerror.SSやerror.dfの情報を使っていないような…
car:::linear.hypothesis.lm
car:::linear.hypothesis.default

…分散分析ばっかりふえるのう…