```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, fig.align="center")
```

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

```{r}
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)
```

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

```{r,message=FALSE}
# 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.

```{r}
g.out <- glht(aov.out, linfct=mcp(trt="Dunnett"))
summary(g.out)
```

By default, R uses the first level of the specified factor as the comparison group, which in our case was 'C', the control group.

```{r}
str(sugar)
levels(sugar$trt)
```

If you (for example) want to compare all groups to the mixed sugar group, we first have to re-order the factor levels.

```{r}
sugar$trt <- relevel(sugar$trt, "G+F")
str(sugar)
levels(sugar$trt)
```

Re-run the analysis.

```{r}
aov.out <- aov(rsp ~ trt, data=sugar)
summary(aov.out)
g.out <- glht(aov.out, linfct=mcp(trt="Dunnett"))
summary(g.out)
```