Anova-Twoway
# twoway.r
# code for two way ANOVA
food <- read.csv('food.csv', as.is=T)
food$type.f <- factor(food$type)
food$sex.f <- factor(food$sex)
# can fit using lm, but more helper functions
# available with aov()
diet.aov <- aov(y ~ type.f + sex.f + type.f:sex.f,data=food)
# note : specifies the interaction
# also have all the usual lm() helper functions
# a shortcut * specifies all main effects and interaction
diet.aov <- aov(y ~ type.f*sex.f,data=food)
# equivalent to first model
anova(diet.aov)
# gives sequential (type I) SS
# but same as type III for balanced data
model.tables(diet.aov)
# tables of means
rat <- read.csv('ratweight.csv',as.is=T)
rat$amount.f <- factor(rat$amount)
rat$type.f <- factor(rat$type)
replications(rat)
# gives number of replicates for each factor
table(rat$amount,rat$type)
# 2 x 3 table of counts for each treatment
rat.aov <- aov(gain ~ amount.f*type.f,data=rat)
# BEWARE: type I (sequential SS)
# to get type III SS, need to declare othogonal contrasts
# can do that factor by factor, but the following does it for all
options(contrasts=c('contr.helmert','contr.poly'))
# first string is the contrast for unordered factors
# the second for ordered factors
rat.aov2 <- aov(gain ~ amount.f*type.f,data=rat)
drop1(rat.aov2,~.)
# drop each term from full model => type III SS
# second argument specifies all terms
drop1(rat.aov, ~.)
# rat.aov() was fit using default contr.treatment
# very different and very wrong numbers if
# forget to use an orthogonal parameterization
# getting marginal means is gruesome!
# model.tables() gives you the wrong numbers
# They are not the lsmeans and not the raw means
# I haven't taken the time to figure out what they are
# easiest way I know is to fit a cell means model
# and construct your own contrast matrices
rat.aov3 <- aov(gain ~ -1 + amount.f:type.f, data=rat)
# a cell means model (no intercept, one X column
# for each combination of amount and type
coef(rat.aov3)
# There is at least one R package that tries to
# calculate lsmeans automatically, but I know
# one case where the computation is wrong.
# (but appears correct).