#-------------------------------------------------------------------# # 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 # Listing 9.1 - 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 with 95% CI") detach(cholesterol) # Listing 9.2 - 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="Q-Q Plot", labels=FALSE) # Assessing homogeneity of variances bartlett.test(response ~ trt, data=cholesterol) # Assessing outliers library(car) outlierTest(fit) # Listing 9.3 - 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) # Listing 9.4 - 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))) # Listing 9.5 - Testing for homegeneity of regression slopes library(multcomp) fit2 <- aov(weight ~ gesttime*dose, data=litter) summary(fit2) # Visualizing a one-way ANCOVA library(HH) ancova(weight ~ gesttime + dose, data=litter) # Listing 9.6 - Two way ANOVA attach(ToothGrowth) table(supp,dose) aggregate(len, by=list(supp,dose), FUN=mean) aggregate(len, by=list(supp,dose), FUN=sd) dose <- factor(dose) 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) # Listing 9.7 - Repeated measures ANOVA with one between and within groups factor CO2$conc <- factor(CO2$conc) 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) # Listing 9.8 - One-way MANOVA library(MASS) attach(UScereal) shelf <- factor(shelf) y <- cbind(calories, fat, sugars) aggregate(y, by=list(shelf), FUN=mean) cov(y) fit <- manova(y ~ shelf) summary(fit) summary.aov(fit) # Listing 9.9 - 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 # Listing 9.10 - Robust one-way MANOVA library(rrcov) Wilks.test(y,shelf, method="mcd") # this can take a while # Listing 9.11 - A regression approach to the Anova problem fit.lm <- lm(response ~ trt, data=cholesterol) summary(fit.lm) contrasts(cholesterol$trt)