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).