Statistics for Laboratory Scientists ( 140.615 )
Analysis of Variance - Multiple Comparisons - Advanced
Example from class: pea sections and sugar
Imagine you want to compare all treatment levels to the control group only. There is a procedure for that called Dunnett’s test, that does that efficiently.
First, read in the data, make a plot, and carry out the ANOVA.
x <- c(75,67,70,75,65,71,67,67,76,68,
57,58,60,59,62,60,60,57,59,61,
58,61,56,58,57,56,61,60,57,58,
58,59,58,61,57,56,58,57,57,59,
62,66,65,63,64,62,65,65,62,67)
treat <- factor(rep(c("C","G","F","G+F","S"), rep(10,5)))
sugar <- data.frame(rsp=x, trt=treat)
set.seed(1)
par(las=1)
stripchart(sugar$rsp ~ sugar$trt, method="jitter", pch=1, xlab="response", ylab="group")
aov.out <- aov(rsp ~ trt, data=sugar)
summary(aov.out)## Df Sum Sq Mean Sq F value Pr(>F)
## trt 4 1077.3 269.33 49.37 6.74e-16 ***
## Residuals 45 245.5 5.46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The procedure for Dunnett’s test is not in the base distribution, but in an add-on package called ‘multcomp’.
https://cran.r-project.org/web/packages/multcomp/index.html
When you use this package for the first time, uncomment the first line to install the package. The ‘library’ function then adds the package to your workspace.
# install.packages("multcomp")
library("multcomp")We then use the ‘glht’ function (general linear hypotheses testing) and specify we want Dunnett’s test for the multiple comparisons procedure.
g.out <- glht(aov.out, linfct=mcp(trt="Dunnett"))
summary(g.out)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: aov(formula = rsp ~ trt, data = sugar)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## F - C == 0 -11.900 1.045 -11.392 < 1e-06 ***
## G - C == 0 -10.800 1.045 -10.339 < 1e-06 ***
## G+F - C == 0 -12.100 1.045 -11.584 < 1e-06 ***
## S - C == 0 -6.000 1.045 -5.744 1.7e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)By default, R uses the first level of the specified factor as the comparison group, which in our case was ‘C’, the control group.
str(sugar)## 'data.frame': 50 obs. of 2 variables:
## $ rsp: num 75 67 70 75 65 71 67 67 76 68 ...
## $ trt: Factor w/ 5 levels "C","F","G","G+F",..: 1 1 1 1 1 1 1 1 1 1 ...levels(sugar$trt)## [1] "C" "F" "G" "G+F" "S"If you (for example) want to compare all groups to the mixed sugar group, we first have to re-order the factor levels.
sugar$trt <- relevel(sugar$trt, "G+F")
str(sugar)## 'data.frame': 50 obs. of 2 variables:
## $ rsp: num 75 67 70 75 65 71 67 67 76 68 ...
## $ trt: Factor w/ 5 levels "G+F","C","F",..: 2 2 2 2 2 2 2 2 2 2 ...levels(sugar$trt)## [1] "G+F" "C" "F" "G" "S"Re-run the analysis.
aov.out <- aov(rsp ~ trt, data=sugar)
summary(aov.out)## Df Sum Sq Mean Sq F value Pr(>F)
## trt 4 1077.3 269.33 49.37 6.74e-16 ***
## Residuals 45 245.5 5.46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1g.out <- glht(aov.out, linfct=mcp(trt="Dunnett"))
summary(g.out)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: aov(formula = rsp ~ trt, data = sugar)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## C - G+F == 0 12.100 1.045 11.584 <1e-04 ***
## F - G+F == 0 0.200 1.045 0.191 0.999
## G - G+F == 0 1.300 1.045 1.245 0.543
## S - G+F == 0 6.100 1.045 5.840 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)