# template do-file for two-way ANOVA
# assume a dataset mydata with outcome y and factors factA and factB

# two-way ANOVA model with interaction (two equivalent versions)
mydata.full <- lm(y ~ as.factor(factA)*as.factor(factB), data=mydata)
mydata.full <- lm(y ~ as.factor(factA)+as.factor(factB)+as.factor(factA):as.factor(factB), data=mydata)
# model information and estimates
summary(mydata.full)
# ANOVA table
anova(mydata.full)

# residuals analysis
mydata.diag <- mydata # all diagnostics to be stored in adv.diag (for convenience, not at all necessary)
mydata.diag$fit <- fitted(mydata.full) 
mydata.diag$stdres <- rstandard(mydata.full)
plot(stdres~fit, data=mydata.diag)
hist(mydata.diag$stdres, nclass=10) # no overlaid normal distribution, note manual choice of no. bins (10)
qqnorm(mydata.diag$stdres)
qqline(mydata.diag$stdres)
# one example of normality test (recall that P-values are only approximate)
shapiro.test(mydata.diag$stdres)

# interaction plot (factor A on x-axis)
library(lattice)
xyplot(fit ~ factA, groups=factB, type="b", data=mydata.diag)

# post-ANOVA inference is cumbersome without using a dedicated packages:
# many exist, e.g. the companion package faraway to books by Julian Faraway

# means for factors, using doBy package to facilitate processing by groups
library(doBy)
summaryBy( y ~ factA, data=mydata, FUN=c(length,mean,sd))
summaryBy( y ~ factB, data=mydata, FUN=c(length,mean,sd))
summaryBy( y ~ factA+factB, data=mydata, FUN=c(length,mean,sd))

# confidence intervals for model parameters (i.e. those listed by summary())
confint(mydata.full)

# pairwise comparisons: require general vector/matrix manipulations, 
# see Exercises 12.54 and 13.31 for examples

# additional model specifications:
# additive model
mydata.red <- lm(y ~ as.factor(factA)+as.factor(factB)+as.factor(factA), data=mydata)

# full model parametrised as one-way ANOVA for combined factors
mydata.full <- lm(y ~ as.factor(as.factor(factA):as.factor(factB), data=mydata)
mydata.full <- lm(y ~ as.factor(as.factor(factA):as.factor(factB)-1, data=mydata) # without intercept
# alternatively one may create a combined factor, such as
mydata$factAB <- as.factor(mydata$factA)+0.1*as.factor(mydata$factB)
mydata.full <- lm(y ~ as.factor(factAB), data=mydata)