# R in Action (2nd ed): Chapter 9
# Analysis of variance
# requires packages multcomp, gplots, car, HH, effects, 
#                   rrcov, mvoutlier to be installed
# install.packages(c("multcomp", "gplots", "car", "HH", "effects", 
                     "rrcov", "mvoutlier"))
#------------------------------------------------------------------

par(ask=TRUE)
opar <- par(no.readonly=TRUE) # save original parameters

# One-way ANOVA
library(multcomp)
attach(cholesterol)
table(trt)     
aggregate(response, by=list(trt), FUN=mean) 
aggregate(response, by=list(trt), FUN=sd) 
fit <- aov(response ~ trt)                                  
summary(fit)
library(gplots)
plotmeans(response ~ trt, xlab="Treatment", ylab="Response", 
          main="Mean Plot\nwith 95% CI")
detach(cholesterol)

# Tukey HSD pairwise group comparisons
TukeyHSD(fit)
par(las=2)
par(mar=c(5,8,4,2)) 
plot(TukeyHSD(fit))
par(opar)

# Multiple comparisons the multcomp package
library(multcomp)
par(mar=c(5,4,6,2))
tuk <- glht(fit, linfct=mcp(trt="Tukey"))
plot(cld(tuk, level=.05),col="lightgrey")
par(opar)

# Assessing normality
library(car)
qqPlot(lm(response ~ trt, data=cholesterol), 
          simulate=TRUE, main="QQ Plot", labels=FALSE)

# Assessing homogeneity of variances
bartlett.test(response ~ trt, data=cholesterol)

# Assessing outliers
library(car)
outlierTest(fit)

# One-way ANCOVA
data(litter, package="multcomp")
attach(litter)
table(dose) 
aggregate(weight, by=list(dose), FUN=mean)
fit <- aov(weight ~ gesttime + dose)                             
summary(fit)

# Obtaining adjusted means
library(effects)
effect("dose", fit)

#  Multiple comparisons using user supplied contrasts
library(multcomp)
contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))
summary(glht(fit, linfct=mcp(dose=contrast)))

# Testing for Homegeneity of Regression Slopes
fit2 <- aov(weight ~ gesttime*dose)
summary(fit2)

# Visualizing a one-way ANCOVA
library(HH)
ancova(weight ~ gesttime + dose, data=litter)

# Two way ANOVA
attach(ToothGrowth)
table(supp,dose)
aggregate(len, by=list(supp,dose), FUN=mean)
aggregate(len, by=list(supp,dose), FUN=sd)
fit <- aov(len ~ supp*dose)
summary(fit)

# plotting interactions
interaction.plot(dose, supp, len, type="b", 
                 col=c("red","blue"), pch=c(16, 18),
                 main = "Interaction between Dose and Supplement Type")
library(gplots)
plotmeans(len ~ interaction(supp, dose, sep=" "),
          connect=list(c(1, 3, 5),c(2, 4, 6)), 
          col=c("red","darkgreen"),
          main = "Interaction Plot with 95% CIs", 
          xlab="Treatment and Dose Combination")
library(HH)
interaction2wt(len~supp*dose)

#  Repeated measures ANOVA with one between and within groups factor	
w1b1 <- subset(CO2, Treatment=='chilled')
fit <- aov(uptake ~ (conc*Type) + Error(Plant/(conc)), w1b1)
summary(fit)
par(las=2)
par(mar=c(10,4,4,2))
with(w1b1, 
     interaction.plot(conc,Type,uptake, 
                      type="b", col=c("red","blue"), pch=c(16,18),
                      main="Interaction Plot for Plant Type and Concentration"))
boxplot(uptake ~ Type*conc, data=w1b1, col=(c("gold","green")),
        main="Chilled Quebec and Mississippi Plants", 
        ylab="Carbon dioxide uptake rate (umol/m^2 sec)")
par(opar)

# One-way MANOVA
library(MASS)
attach(UScereal)
y <- cbind(calories, fat, sugars)
aggregate(y, by=list(shelf), FUN=mean)
cov(y)
fit <- manova(y ~ shelf)
summary(fit)
summary.aov(fit)

#  Assessing multivariate normality
center <- colMeans(y)
n <- nrow(y)
p <- ncol(y)
cov <- cov(y)
d <- mahalanobis(y,center,cov)
coord <- qqplot(qchisq(ppoints(n),df=p),
                d, main="QQ Plot Assessing Multivariate Normality",
                ylab="Mahalanobis D2")
abline(a=0,b=1)
identify(coord$x, coord$y, labels=row.names(UScereal))

# multivariate outliers
library(mvoutlier)
outliers <- aq.plot(y)
outliers

# Robust one-way MANOVA
library(rrcov)
Wilks.test(y,shelf, method="mcd")