```{r echo = FALSE, results = 'hide', message = FALSE}
library(tidyverse)
library(mosaic)
library(ggthemes)
theme_set(theme_solarized())                 # change theme for ggplot2
trellis.par.set(theme=theme.mosaic()) # change default color scheme for lattice

knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    tidy    = FALSE)
```

## ANOVA Models and Pairwise Comparisons of Means

Let's do an experiment to examine what happens when we do lots of comparisons between pairs of means in the setting of a predictor variable with lots of categories.

First, let's suppose that we have a categorical predictor variable with 12 levels (that is, our data comes from 12 groups), and a quantitative response variable.

In reality, **the predictor is unrelated to the response**. 

In other words, **in the population / in the long run**, the mean of the response variable is **the same regardless of which of the twelve groups the case is in**.

However, when conducting a study, we don't know this reality going in, and want to use the data to learn about possible differences among the groups. That is, we want to test against a null hypothesis that says that the 12 population/long run means are identical.

We collect data from this population, with 20 cases in each of the 12 groups.


### Generating Simulated Data from a Known Population Model

The following R code will create a simulated dataset from the reality described above.

```{r}
## First we set a random seed for reproducibility
set.seed(47)
## Then we generate 12 sets of 20 observations, all from a common N(100,15) distribution
FakeData <- tibble(
  Group    = rep(c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"), each = 20),
  Response = rnorm(20 * 12, mean = 100, sd = 15))
```

Let's take a look at the dataset to verify that it looks as expected.

```{r}
FakeData
```


### Descriptive Visualization and Statistics

We can plot the data by group using side-by-side box plots, adding the means to the plot, and get quantitative summary statistics from the dataset. 

```{r, warning = FALSE}
gf_boxplot(Response ~ Group, data = FakeData) +
  stat_summary(fun = mean, color = "darkorange")
```

We can also examine a table of summary statistics, and arrange the groups in order of the sample means:

```{r}
favstats(Response ~ Group, data = FakeData) %>% 
  format(digits = 4) %>%
  arrange(mean)
```

It looks like group D had the largest sample mean, at 108.44, and group B had the smallest, at 95.13.

### Fitting an ANOVA Model

Since we have a dataset with a categorical predictor (`Group`) and a quantitative response (`Response`), it makes sense to fit an ANOVA model to the data. Recall that this model has the form:

$$ Response_i = \mu_{Overall} + \alpha_{Group_i} + \epsilon_i $$

```{r}
aovModel <- aov(Response ~ Group, data = FakeData)
```

Here are the estimates of the $\alpha$ parameters

```{r}
aovModel %>%
  model.tables(type = "effects")
```

We can construct the ANOVA table and perform an overall F-test to see how strong the evidence is that _any_ of the group means differ:

```{r}
summary(aovModel)
```


### Pairwise Comparisons of Means

Suppose instead of fitting a model to all the data, we had first looked at the data and then decided to do a hypothesis test to see how strong the evidence was for a difference in population means between groups B and D, since these appeared to be the groups with the most promising evidence for a difference. What would happen?

```{r}
FakeData %>%
  filter(Group %in% c("B","D")) %>%
  lm(Response ~ Group, data = .) %>%
  summary()
```

We appear to have "strong evidence" that the means of the populations that gave us the data in groups B and D are different.

We can actually get $P$-values for all possible pairwise comparisons of pairs of group means as follows. 

Since there are 12 groups, and each group can be compared to each other group, there are 12 * 11 = 132 ordered pairs of means; but since this counts each pair twice, there are actually 66 unique pairs.

```{r}
PairwiseTests <- with(
  FakeData, 
  pairwise.t.test(
    x = Response, g = Group, 
    p.adjust.method = 'none'))
```

```{r}
PairwiseTests
```


Let's make this a bit easier to read by rounding the P-values:

```{r}
PairwiseTests %>% 
  pluck("p.value") %>%
  round(digits = 3)
```

Looking over the P-values for these pairwise tests, it looks like there are several "significant differences" at the 5% significance level between groups:

A and D
B and D
B and F
B and H
B and K
D and L

with a several more pairs that have P-values just a bit above 0.05 (A and H, D and E, D and J, J and H, H and L).

So we have two results: an overall F test that says there isn't good evidence that the population means are any different (which is a desirable result in this case, since in fact the process producing the data did not differ _at all_ across groups), and a collection of t-tests, six of which yielded significant evidence of a difference between means.

**Pairs:** 

1. How should we make sense of this discrepancy?

```{r}
linearModel <- lm(Response ~ Group, data = FakeData)
```

```{r}
summary(linearModel)
```


2. Given that we know exactly the "population" model that generated the data, what can we say about the number of Type I and Type II errors (false discoveries and missed discoveries) in the results?

3. There are 12 * 11 / 2 = 66 sets of pairs of groups and therefore 66 different t-tests being performed here. If you had a dataset like this handed to you and did all pairwise tests, how many times do you **expect** you would reject $H_0$ mistakenly?


### Adding a genuinely different mean

Suppose we add a 13th group whose population mean really _is_ different from the others: 110 instead of 100.

```{r}
NewData <- tibble(
  Group    = rep("M", 20),
  Response = rnorm(20, mean = 110, sd = 15))
NewData
```

Let's join this dataset with the original one.

```{r}
CombinedData <- FakeData %>%
  bind_rows(NewData)
```

Repeating the boxplot code:

```{r}
gf_boxplot(Response ~ Group, data = CombinedData) +
  stat_summary(fun = mean, color = "darkorange")
```

Let's refit the ANOVA model on this new dataset.

```{r}
aovModel2 <- aov(Response ~ Group, data = CombinedData)
```

```{r}
summary(aovModel2)
```

Now we have significant evidence of a difference between some population mean and some other population mean. And now, it's true that some population mean is different!

So suppose we had this data did this test. What would we do next? How do we know _which_ means differ from which other means?

If we repeat our collection of pairwise tests from before:

```{r}
PairwisePValues <- with(
  CombinedData, 
  pairwise.t.test(
    x = Response, g = Group, 
    p.adjust.method = 'none')) %>%
  pluck("p.value") %>%
  round(digits = 3) 
PairwisePValues
```

```{r}
ifelse(PairwisePValues < 0.05, "*", "")
```

The same six comparisons that were significant before still are. In addition we have significant evidence of a difference between the population mean for group M and the population means for groups A, B, C, E, G, I, J, and L (but not D, F, H, and K)

**Pairs:**

4. How many Type I and Type II Errors are there here?

## Using adjusted P-values and Confidence Levels

If we employ a Bonferroni correction to the P-values and confidence levels...

```{r}
library(DescTools)
PostHocTest(aovModel2, method = "bonferroni")
```

none of the spurious differences come out significant, but only one of the 12 real differences does.

Using Fisher's LSD, we get the same results we got before by just doing the t-tests, since our only "protection" was the overall F-test, which was significant.

```{r}
PostHocTest(aovModel2, method = "lsd")
```

Finally, using Tukey's HSD:

```{r}
PostHocTest(aovModel2, method = "hsd")
```

we get the same results in this case as Bonferroni (although the P-vales are smaller in general).