---
title: "Analysis of variance"
author: "Jonathan Harrington, Johanna Cronenberg"
output: 
  bookdown::html_document2:
    number_sections: TRUE
    toc: true
    theme: flatly
    highlight: pygments
---

<style>
div.col {background-color: #ffd5c0; border-radius: 5px; padding: 20px;}
body {font-size: 16pt;}
h1 {font-size: 24pt;}
h2 {font-size: 22pt;}
p.caption {font-size: 12pt; text-align: justify;}
code.sourceCode {font-size: 16pt;}
</style>

If not already done, carry out parts 1-5 of the setup,
[as described here](https://www.phonetik.uni-muenchen.de/studium_lehre/lehrmaterialien/R_speech_processing/setup.html).

Load these libraries and data-frames:

```{r, message=FALSE}
library(tidyverse)
library(broom)
library(afex)
library(emmeans)
urla = "https://www.phonetik.uni-muenchen.de/"
urlb = "studium_lehre/lehrmaterialien/R_speech_processing/Rdf/"
url = paste0(urla, urlb)
e.df = read.table(file.path(url, "e.txt"), 
                  stringsAsFactors = T)
form.df = read.table(file.path(url, "form.df.txt"),
                     stringsAsFactors = T)
dg.df = read.table(file.path(url, "dg.txt"), 
                   stringsAsFactors = T)
ssb.df = read.table(file.path(url, "ssb.txt"), 
                   stringsAsFactors = T)
```

# Analysis of variance: general considerations.

An analysis of variance (ANOVA) is an extension of a *t*-test to cases involving several levels e.g.:

- 'Do /p, t, k/ differ in voice onset time?'. 
Factor: place of articulation, **three levels** /p, t, k/.

and/or to multiple factors e.g.

- 'Do /p, b/ in stressed and unstressed syllables differ in voice onset time?' which is the same as: 'Is voice onset time affected by stop voicing and stress'. 
**Factor 1: stop voicing**, levels /p, b/; **factor 2: stress**, levels 'stressed', 'unstressed'.

<div class="col">
An ANOVA is used to calculate the relationship between two type of variation: 

- the systematic variation **between** levels. 

- the random variation within levels.

This is schematically outlined in the example figure below of first formant (F1) frequencies of three vowels (i.e. three levels) /ɪ,ɛ,a/

![](img/anova.png)


- F1 varies systematically **between** the three vowels. This systematic variation comes about for acoustic-phonetic reasons by which F1 increases with phonetic vowel height through  /ɪ,ɛ,a/ respectively.

- F1 varies **within** each level (in this case vowel): that is,  there is random F1-variation in /ɪ/, random F1-variation in /ɛ/, and random F1-variation in /a/. 

In an ANOVA, the *F*-ratio, which is the ratio of these two quantities, is tested for significance.

$$ Fratio = \frac{systematic-variation-BETWEEN-the-levels}{random-variation-WITHIN-the-levels} $$

- If there is a high degree of systematic separation between the vowels, then the numerator in the above equation has a high value (causing an *increase* in the *F*-ratio).

- If the random variation within /ɪ/ and/or within /ɛ/, and/or within /a/ is high, then the the denominator is also high (causing a *decrease* in the *F*-ratio).

Thus, the greater the systematic differences between the vowels, and the less the random variation within vowels, the greater the *F*-ratio.

ANOVA then tests whether the *F-ratio* is significantly greater than chance (significantly greater than 1) using an *F*-test.
</div>

# Relationship to t-test.

As when running a *t*-test, it is essential to understand how
[within and between subjects factors are different](https://www.phonetik.uni-muenchen.de/studium_lehre/lehrmaterialien/R_speech_processing/stat_phonetics_intro/topics/02_statistics.html#between-and-within-subjects-factors)

## ANOVA and a paired t-test.
Recall from
[here](https://www.phonetik.uni-muenchen.de/studium_lehre/lehrmaterialien/R_speech_processing/stat_phonetics_intro/topics/02_statistics.html#a-paired-t-test)
that a paired *t*-test was used when for a within-subject factor. When testing whether `F2` was influenced by the within-subject factor`Stress`,  the command was:

```{r}
form.df %>%
  pivot_wider(names_from = Stress, 
              values_from = F2) %>%
  mutate(F2diff = stressed - unstressed) %>%
  pull(F2diff) %>%
  t.test()
```
which was reported as:

- There was a significant influence of stress on F2 (t[11] = 4.4, p < 0.01).

The corresponding command with an ANOVA is:

```{r}
form.aov = 
  aov_4(F2 ~ Stress  + (Stress|Subject), data = form.df)
form.aov
# or with dot notation
form.df %>%
  aov_4(F2 ~ Stress  + (Stress|Subject), .)
  
```
<div class="col">
`(Stress|Subject)` is the syntax for indicating that `Stress` is a within-subjects factors in the `aov_4()` function.
</div>

The most important parts of the output are:

- two types of degrees of freedom.
- an *F*-value. 
- a significance level and probability value.

This should be reported as:

- There was a significant influence of stress on F2 (F[1, 11] = 18.96, p < 0.01).

<div class="col">
The choice of p < 0.01 is because of the ** which means that the probability under `p.value` is between 0.001 and 0.01. The `p.value` in this output is rounded to three decimal places. To get the actual value, pass the output of `aov_4()` to `anova()` i.e. in this case:

```{r}
anova(form.aov)
```

The *F*-value obtained from the ANOVA is the *t*-statistic squared from the paired *t*-test. Thus:

```{r}
# square of the t-statistic
(-4.3543)^2
# (same as the F-value. other than rounding differences).
```
Note also that (rounding differences apart) the probability values from the *t*-statistic (`p-value` = 0.001147) and from `Pr(>F)` in the ANOVA are the same.

The output of  the above ANOVA also gives a column `ges` which is *generalised-eta-squared*. It can be used to give an indication of the size of the effect found. More details on this statistic
[can be found here](https://link.springer.com/content/pdf/10.3758%2FBF03192707.pdf).
The paper suggests that a generalised-eta-squared value of over 0.26 is large, and less than 0.02 is small. A small effect might be one in which e.g. a study showed a significant change in F2 of only 5 Hz. This effect size is very small in relation to the typically value of F2 (between 220 and 3000 Hz) and is probably meaningless even though significant (because F2 changes less than 10 Hz are probably not perceptible). In this hypothetical case, the `ges` values is likely to be small.
</div>


## ANOVA and an unpaired t-test.
Recall also from
[here](https://www.phonetik.uni-muenchen.de/studium_lehre/lehrmaterialien/R_speech_processing/stat_phonetics_intro/topics/02_statistics.html#the-unpaired-t-test)
that a unpaired *t*-test was used for a between-subject factor. When testing whether `F2` was influenced by the between-subject factor `Sprache`, the command was:

```{r}
e.df %>%
t.test(F2 ~ Sprache, .) %>%
  tidy()
```

and reported as:

- F2 was significantly influenced by language (t[21.1] = 2.3, p < 0.05).

The corresponding command with an ANOVA is:

```{r}
e.df %>%
aov_4(F2 ~ Sprache + (1|Vpn), data = .)
```

and can be reported as:

- F2 was significantly influenced by language (F[1,25] = 5.36, p < 0.05).

<div class="col">
Information is reported *Contrasts set to contr.sum for the following variables: Sprache*. This is a type of coding of the levels that is commonly used in ANOVAs often known as *deviation coding*: for this example, it compares the mean F2 for each level of `Sprache` ito the overall mean of the F2 in entire data-frame.  More details
[are given here](https://stats.oarc.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/#DEVIATION).

Different results are produced for the unpaired *t*-test and ANOVA. This is because in an unpaired *t*-test,  the 
[variances of the levels are *not* assumed to be equal](https://www.phonetik.uni-muenchen.de/studium_lehre/lehrmaterialien/R_speech_processing/stat_phonetics_intro/topics/02_statistics.html#the-unpaired-t-test)
(resulting in fractional degrees of freedom) whereas in an ANOVA they are (hence integer degrees of freedom). The same result is obtained if the *t*-test is run under the assumption of equal variance which
[as explained in the final note here](https://www.phonetik.uni-muenchen.de/studium_lehre/lehrmaterialien/R_speech_processing/stat_phonetics_intro/topics/02_statistics.html#the-unpaired-t-test)
is given by:

```{r}
e.df %>%
t.test(F2 ~ Sprache, ., var.equal=T)
```
The *t*-statistic squared, degrees of freedom, and probability value from this unpaired *t*-test under the assumption of equal variance are now the same as those obtained with the above ANOVA.

</div>

# Running an ANOVA with more than one factor.
The question to be resolved for the data-frame `dg` is whether `F2` is influenced by region (`Region`) and speaker sex (`Gen`). Here then there are *two* (between-subject) factors: `Region` has three levels and `Gen` has two levels.

<div class="col">
There are therefore two reasons why this problem has to be accomplished with an ANOVA (and why a *t*-test is insufficient):

- there is more than one factor.
- one of the factors has more than two levels.
</div>

```{r}
summary(dg.df)
```
The steps in carrying out the ANOVA are:

1. identify as before whether these are between- or within-factors.

2. make a plot of the data.

3.  assess whether there are *interactions* (which is dealt with separately below).

4. run `aov_4()`

5. if there is an interaction, run so-called *post-hoc tests*.

1-4 are considered below; 5 is further explored at the end of this document.

## Between and within factors.
Both are likely to be between factors, because they deal with attributes of the speaker. For the two-factor case, this can be most easily seen as follows:

```{r}
dg.df %>%
mutate(Factors = interaction(Region, Gen)) %>%
  select(Vpn, Factors) %>%
  table()
```


Because for all speakers there is a single 1 and remaining zeros, these are between-subject factors (i.e. a speaker is e.g. from region A and is male and therefore not from the other two regions nor female).

## A plot of the data.

```{r}
dg.df %>%
  ggplot+
  aes(y = F2, x = Region, col = Gen) +
  geom_boxplot()
```

In an analysis with two factors, three questions are considered:

- is F2 influenced by `Gen`?
- is F2 influenced by `Region`?
- is there an interaction between `Gen` and `Region`?

For the first two questions, the plot suggests:

- F2 is higher for females than for males
- F2 is about the same in `A` and `B` and both of these have a higher F2 than `C`.

## Checking for an interaction.
The third question - which is to do with the interaction - means either:

- is the difference between regions the same for males and females?

and:

- is the difference between males and females the same in all three regions?

The plot suggests that the answer to this third question is 'no'. This is because the male-female difference is much less in region `C` than in either of the other two regions. Therefore, the statistics are likely to show  a significant interaction. No interaction between the factors would occur if the male-female differences in regions A, B, and C were about the same (which is not the case).

## Running `aov_4()` with two factors.
The syntax is much as before except that two factors `A` and `B` are now specified as `A * B` (or for an equivalent result as `A + B`). The specification for the subject factor follows from the earlier examples:

- If (as in this example) both are between factors:  `(1|Vpn)`.
- If one factor is within: `(A|Vpn)`.
- If both factors are within: `(A + B|Vpn)` or `(A * B|Vpn)`

For the present example, therefore:

```{r}
dg.df %>% 
  aov_4(F2 ~ Region * Gen + (1|Vpn), .)
```

The results show:

- F2 was significantly influenced by `Region` (F[2,54] = 119.64, p < 0.001).
- F2 was significantly influenced by `Gen` (F[1,54] = 106.15, p < 0.001).
- There was a significant (F[2,54] = 12.08, p< 0.001) interaction between these factors.

All of these results were expected from the plot.

<div class="col">
AN ANOVA can only be run if for each speaker there is:

- ONE value for EACH level of each within factor.
- ONE value for ONE level of each between factor.

This makes it possible to calculate in advance exactly how many observations are input to an ANOVA, given the number of speakers, and a knowledge of the factors (whether between or within) and the number of levels. For example, suppose there are 12 speakers, 7 old and 5 young who produce three different words. There is, then, one between factor (Age: two levels, young, old) and one within factor (Word: three levels, word1, word2, word3). Therefore there must be 12 x 1 (Between) x 3 (within) = 36 observations in total (not more, not less) for the ANOVA to be applicable. Phonetics experiments usually deviate from this principle because they typically collect *multiple repetitions* of the same factor or factor combinations. An example of this type of data is in the data-frame `ssb.df`. Here younger and older speakers (between factor `Alter`: two levels) produced `swoop`, `used` `who'd` (within factor `Wort`: three levels) *up to 10 times* and the second formant frequency (`F2`) was measured.

```{r}
ssb.df %>%
  mutate(both = interaction(Wort, Alter)) %>%
  select(Vpn, both) %>%
  table()
```

Data such as these have to be aggregated (averaged) so that there is only one (and not 10) value for each level of the within factor (and one value for one of the levels of the between factor):

```{r}
ssb.m.df =
  ssb.df %>%
  group_by(Alter, Wort, Vpn) %>%
  summarise(F2 = mean(F2)) %>%
  ungroup()
```

Now the data are ready for the ANOVA because there is one value for each level of the within factor per speaker:

```{r}
ssb.m.df %>%
  mutate(both = interaction(Wort, Alter)) %>%
  select(Vpn, both) %>%
  table()
```

In fact, the `aov_4()` function will do the averaging, as shown here:
```{r}
ssb.df %>%
aov_4(F2 ~  Wort * Alter + (Wort|Vpn), .)
```
and give the warning "More than one observation per design cell, aggregating data using `fun_aggregate = mean`." Because `aov_4()` has taken care of the averaging, the above is exactly equivalent to:

```{r}
ssb.m.df %>%
aov_4(F2 ~  Wort * Alter + (Wort|Vpn), .)
```

An ANOVA can also not be run if there is missing data. Suppose for example, we remove observation 5 from the averaged data: then speaker `mapr` will only have produced `used` and `who'd` but not `swoop`. This is evident from the zero in `mapr` for `swoop`:

```{r}
ssb.m.df %>%
  # remove observation 5
  slice(-5) %>%
  mutate(both = interaction(Wort, Alter)) %>%
  select(Vpn, both) %>%
  table()
```

If `aov_4()` is run on this data with the missing value, then it removes *all* of that speaker's data so that there are once again no missing values:

```{r}
 ssb.m.df %>%
  # remove observation 5
  slice(-5) %>%
  aov_4(F2 ~  Wort * Alter + (Wort|Vpn),  .)
```

On a different point, the output also reports (in black) "Sphericity correction method: GG". The abbreviation "GG" stands for *Greenhouse-Geisser*. This correction only ever applies if a within-factor has more than two levels (as in the case of `Wort` here). The correction is to compensate for the different variances of the within-factor levels. The compensation is done by changing the degrees of freedom often to a fractional value (which is why `Wort` and `Alter:Wort` both have degrees of freedom [1.37, 13.72]) and adjusting the probabilities upwards (so that it is slightly more difficult to achieve significance). Note that this correction never applies to any between-factor.

 
## An introduction to the `emmeans()` function and post-hoc tests.
Following a significant interaction, further tests should be carried out to test which **combination(s) of levels** is/are significant. The package of choice for doing so is `emmeans`. As its name suggests, the post-hoc test is by default based on **e**stimated **m**arginal **m**eans that are also often called least squares or adjusted means.

<div class="col">
Estimated marginal means are based on so-called *marginal means* (MM) which are the means of the *levels* of factors. As explained
[in further details here](https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html)
they are used to correct for differences in the number of observations between levels. In a simple data-frame such as the following modelled on
[an example from here](https://webpages.uidaho.edu/CALS-StatProg/sas/workshops/glm/handout2.1_glm.pdf):

```{r}
x = c(4, 6, 2, 7, 3, 5, 4, 2, 3)
dialect = c(rep("A", 3), rep("B", 3), rep("C", 3))
x.df = data.frame(x, dialect)
x.df
```
the data-frame is 'balanced' because there are three observations for each level of dialect. In this special case, the arithmetic and marginal mean (MM) are the same. The arithmetic mean is just the mean of the numerical variable across all observations:
```{r}
x.df %>% 
  pull(x) %>%
  mean()
```
and the MM is the mean of each of the three means, one per level:
```{r}
x.df %>%
  # means of "A", "B", "C"
  group_by(dialect) %>%
  summarise(m = mean(x)) %>%
  # mean of these three means
  pull(m) %>%
  mean()
```
That is, both are 4. Now consider the same calculation but after removing observation 4. In this case, the data-frame is 'imbalanced' because whereas levels `A` and `C` each have three observations, `B` now has only two. Because this data-frame is imbalanced, the arithmetic mean and MM will be different:
```{r}
# same calculations as above but remove obs. 4
x.df %>% 
  slice(-4) %>%
  pull(x) %>%
  # Arithmetic mean
  mean()

x.df %>%
  slice(-4) %>%
  group_by(dialect) %>%
  summarise(m = mean(x)) %>%
  pull(m) %>%
  # MM
  mean()
```
In post-hoc tests with the EMM package, the marginal means are *estimated* (hence EMM) by compensating for such imbalances. According to
[the vignette](https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html),
using EMMs "is a much fairer way to compute marginal means, in that they are not biased by imbalances in the data. We are, in a sense, estimating what the marginal means *would be*, had the experiment been balanced'. For this reason, EMMs are typically based on a *model* of the marginal means (and not on the actual means themselves).
</div>

For the present example, there are nine tests to be carried out. 

1. Gender: Do males and females differ within each of the three regions? (3 tests).

- male vs. female in region A
- male vs. female in region B
- male vs. female in region C

2. Region: Do regions differ within each sex?
(6 tests).

- region A vs. region B in males
- region A vs. region C in males
- region B vs. region C in males
- region A vs. region B in females
- region A vs. region C in females
- region B vs. region C in females

These nine tests can be carried out as follows
```{r}
dg.df %>% 
  # Anova
  aov_4(F2 ~ Region * Gen + (1|Vpn), .) %>%
  # emmeans to test for interaction
  emmeans(., ~ Region * Gen) %>%
  # the results of the post-hoc tests
  pairs(., simple = "each", combine=T)
```
The first command with `emmeans()` calculates the estimated marginal means and their 95% confidence intervals. The second command with `pairs()` runs the required nine tests. The first line in the output shows the difference between regions A and B for male speakers. The last line shows the difference between male and female speakers for region C. 

An appropriate summary of the output might be as follows.

Post-hoc Bonferroni corrected tests showed a significant difference between regions A and C (p < 0.001) and between regions B and C (p < 0.001) for both male and female speakers. There were no differences between regions A and B for either male nor female speakers. The post-hoc tests also shows significant differences between male and female speakers within regions A (p < 0.001) and B (p < 0.001) but not within C. 

<div class="col">
Some further information on the output and `emmeans()` is as follows.

1. **Bonferroni correction for multiple tests**.
The information in the last line of the output *P value adjustment: bonferroni method for 9 tests* means that the probabilities have been adjusted to take into account that multiple tests have been run. The adjustment is made because running multiple tests increases the likelihood that one or more of them will be significant by chance i.e., of wrongly rejecting the null hypothesis (sometimes known as a Type I error). The correction is to multiply the resulting probability by the number of tests that are made - or equivalently to reset the level at which the null hypothesis is rejected from 1-*α* (e.g. 0.05) to (1-*α*)/9 (0.0055) for 9 tests as in the above example. Everything else remains the same. Notice that when the same test is run *without Bonferroni correction* which can be done with `adjust=none`, then the resulting probabilities are 1/9th of those as when Bonferroni corrected:

```{r}
dg.df %>% 
aov_4(F2 ~ Region * Gen + (1|Vpn), .) %>%
  # emmeans to test for interaction
  emmeans(., ~ Region * Gen) %>%
  # no Bonferroni correction
  pairs(., simple = "each", combine=T, adjust="none")
```

As discussed 
[here](https://cran.r-project.org/web/packages/emmeans/vignettes/confidence-intervals.html), 
the Bonferroni test is quite conservative meaning that it is often more difficult to achieve significance than with other post-hoc tests. It can, therefore, be helpful to try other methods for post-hoc testing, such as `adjust=mvt` or `adjust=holm`. There is no unambiguous answer as to which post-hoc test is the correct one to use.

2. **The arguments to `pairs()`**.

Leaving out the last two arguments causes all the possible pairwise tests to be listed. In this case, there are 15 such tests (15 = 6!/4!2! where ! is factorial; 6! because there are 3 levels in `Dialect` x 2 levels in `Gen`):
```{r}
dg.df %>% 
aov_4(F2 ~ Region * Gen + (1|Vpn), .) %>%
  # emmeans to test for interaction
  emmeans(., ~ Region * Gen) %>%
  # no Bonferroni correction
  pairs(.)
```
However, combinations in which levels of both factors are different such as comparing male-B with female-C are not usually needed. Such combinations are removed with `simple="each"`:
```{r}
dg.df %>% 
aov_4(F2 ~ Region * Gen + (1|Vpn), .) %>%
  # emmeans to test for interaction
  emmeans(., ~ Region * Gen) %>%
  # no Bonferroni correction
  pairs(., simple="each")
```
to give 9 tests. The last argument `combine=T` gives the appropriate Bonferroni correction (thus changing the probabilities) but is in all other respects the same.
</div>

