Lab-11-Repeated-Measures-ANOVA-With-Code.Rmd

---
title: "Lab-11-Repeated-Measures-ANOVA"
author: "Matt Defenderfer"
date: "10/5/2019"
output: 
  html_document:
    toc: TRUE
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Overall Progression for Lab 11
1. Load Data and Packages
2. Assumptions for Repeated Measures ANOVA
3. Running a Repeated Measures ANOVA
4. Reporting a Repeated Measures ANOVA

## Load Data and Packages

We will be using the `car`, `pastecs`, `ez`, and `tidyverse` packages for the same functions we needed for the independent one-way ANOVA.

```{r Load Packages, message = FALSE}
library(car)
library(pastecs)
library(ez)
library(tidyverse)

theme_set(theme_bw())
```

The dataset we will be using is in "Repeated-Measures-ANOVA.Rdata". This file contains a single dataframe, `noise`, describing performance on a perceptual task with varying levels of background noise. `noise` has 5 columns:

- `ID` (factor, 20 levels): ID number for each subject
- `none` (numeric): score on the task with no background noise
- `low` (numeric): score on the task with low level of background noise
- `med` (numeric): score on the task with medium level of background noise
- `high` (numeric): score on the task with high level of background noise

```{r Load the data}
load('/Users/Matt/Google Drive/Grad School/GR 770 Statistics/R Labs/Data/Repeated-Measures-ANOVA.RData')
```

`none`, `low`, `med`, and `high` all are measures of the same task in the same participants under different conditions and constitute a repeated-measures design. We know that some of our functions, like all of our plots, behave best when there is only one observation per row, so we should gather the different conditions and different scores into single columns.

```{r gather} 
noise.g <- gather(noise, key = 'noise.level', value = 'score', none:high)

noise.g$noise.level = factor(noise.g$noise.level, levels = c('none','low','med','high'), ordered = TRUE)
```

Now that we have both wide and long formats for our data, we can move on to the assumptions

## Assumptions for Repeated Measures ANOVA

1. Normality within Groups (test before running the ANOVA)
2. Sphericity (test while running the ANOVA)

### Normality
As always, we will be testing normality with histograms, boxplots, q-q plots, and statistical tests

#### Histograms
```{r histogram}
ggplot(noise.g,aes(x = score, fill = noise.level)) +
  geom_histogram(bins = 5, color = "black") +
  facet_wrap(~noise.level, scales = "free")
```

Overall, there isn't much to worry about. None and Low look fine since the largest bar is in the center and dropoff is consistent on both sides of the plot. Med is slightly worrying since the bars increase in size until the last one, but nothing overtly skewed. High is slightly concerning since the middle bar is flanked by two higher bars, however the end bars are both quite low. Nothing about the histogram immediately indicates skew for any single group

#### Q-Q Plots
```{r Q-Q}
ggplot(noise.g,aes(sample = score)) +
  geom_qq() +
  geom_qq_line() +
  facet_wrap(~noise.level, scales = "free")
```

None, Low, and Med are all pretty much perfect in the Q-Q plots. All the dots are close to the line, and although there are slight curves in the dot pattern, they are nothing to worry too much about. High is definitely not uniformly shaped, but none of the dots are too far away from the line. Special attention should be paid to this group in the statistical tests considering its histogram result. 

#### Boxplots
```{r Boxplot}
ggplot(noise.g,aes(x = noise.level, y = score)) +
  geom_boxplot()
```

None and Med are both perfectly fine. Median line is close to the center of the IQR in both cases and the whiskers are fairly even. Low has a single outlier, but otherwise also looks great. The outlier isn't extremely far away either, so it should be alright. The whiskers for High are even, but the median line is fairly low. Overall, though it looks more normal in this graph than in either the histogram or the Q-Q plot.

#### Statistical Tests

```{r stats desc}
list.out <- tapply(noise.g$score, noise.g$noise.level, stat.desc, basic = FALSE, norm = TRUE)
norm.stats <- as.data.frame(do.call(rbind, list.out))
t(norm.stats) # the t() function transposes the dataframe to make the output look nicer
```

Skewness is less than 0.5 for all of our groups and skew.2SE is very far away from the 1.58 cutoff. None of the SW tests came out significant either indicating normality. Taking the smaller sample size into account, more weight should be given to the statistical tests than especially the histogram. The Q-Q plot is somewhat concerning for the High group, but it looks normal elsewhere. Overall, I would say all of the groups are normal enough to continue without any transformations

### Sphericity
Sphericity is testing automatically when running both `Anova` and `ezANOVA` so there is no reason to run it separately. There is also not any specific requirement that the dataset be spherical. If the dataset isn't spherical, then all that changes is the degrees of freedom at the end. We can skip this for now.

## Running a Repeated Measures ANOVA

### car Anova Method
Unfortunately, repeated measure ANOVAs using `Anova` are quite a bit more annoying to set up than a typical between-subjects Anova. To begin with, you want your data in wide format, where all of the repeated measures have their own separate column. For this reason, we are going to use the original `noise` data frame as opposed to the gathered version. Essentially, the different levels of the repeated measure are going to be put in as different dependent variables in the model and bound together using a function called `cbind`. Any between-subjects independent factors will be input as usual. If there are no between-subjects factors, put a on the right hand side of the formula instead.

For example:

```{r car model}
# Use noise dataframe instead
noise.lm <- lm(cbind(none,low,med,high) ~ 1, noise)
```

Next, we need to set up a design matrix giving each unique combination of dependent variable levels. So if we had two within-subject variables, we would create a dataframe where each row was a unique combination of each level in those two factors. When there is only a single within-subjects factor, we just need to make a dataframe with each unique level of that factor. We will do that here, store the dataframe in `idata` and name the created variable `noise.level`. 

```{r idata}
# set up design matrix
idata <- data.frame(noise.level = factor(c('none','low','med','high')))
```

We will be passing this in to the `Anova` function. The other new input to the `Anova` function is the `idesign` variable. This will set up the within-subjects effects we are testing. We are saying we want a design as a function of the within-subjects variable(s) and are also saying whether we want to test interactions or just main effects of each within-subjects variable. So here, since we only have a single within-subjects variable, we will just set `idesign` to be a function of that variable. We will also be setting our sum of squares to be type 3 and storing the output in `noise.Anova`

```{r run the Anova}
# run the Anova
noise.Anova <- Anova(noise.lm, idata=idata, idesign = ~noise.level, type = 3)
```

The last step is actually viewing the results. This time, we are going to be using the `summary` function passing in the output from `Anova`. However, we also need to set another variable, `multivariate`, to be false. This method is actually also how you do a multivariate ANOVA, so by default `summary` will output both univariate and multivariate results. However, we know we only want univariate results, so we will set `multivariate` to be FALSE. 

```{r View Anova results}
summary(noise.Anova, multivariate = FALSE)
```

Here, you see the output table for a repeated measures ANOVA along with tests for sphericity and the corresponding sphericity correction factors. 

The first part of the table will give you sums of squares, numerator and denominator degrees of freedom, the F-statistic, and the p-value for each independent variable in the ANOVA. As you can see, mean score was significantly different across all levels of `noise.level`. 

The second part of the table gives output for Mauchly's Test for Sphericity. It gives you the test statistic and the p-value. If the p-value is significant, that means sphericity has been violated. This tells us that the variances of differences in the levels are significantly different from one another. Fortunately, you do not need to go back and change anything about your data if sphericity is violated. If the outcome of sphericity is significant, you will move to the third part of the table, the correction factors for departure from sphericity.

By default `Anova` outputs both Greenhouse-Geisser and Huynh-Feldt corrections for each main and interaction effect with each within-subjects variable. If sphericity is violated, do the following:

1. Look the at the Greenhouse-Geisser $\epsilon$ value (the `GG eps` output in the `Anova` table)
2. If the GG $\epsilon$ value is < 0.75, multiply both the numerator and denominator degrees of freedom by the Greenhouse-Geiser $\epsilon$ value when you report the results.
3. If the GG $\epsilon$ value is > 0.75, instead multiply the numerator and denominator degrees of freedom by the Huynh-Feldt $\epsilon$ value when you report the results.

IF THE P-VALUE FOR MAUCHLY'S TEST IS NON-SIGNIFICANT, DO NOT ALTER THE REPORTED DEGREES OF FREEDOM!!!!!!

Now, you can take your results, throw them into a table and write your report sentence. We will talk about the form this will take here in just a minute. 

### ezANOVA Method
Luckily, the ezANOVA method is much simpler and still gives us the same outputs, albeit in an uglier format.

Like we did before, we will take the gathered dataframe where each dependent and independent variable has its own column and throw it into the `ezANVOA` function. This time, I am wrapping the name of the dataframe in the `data.frame` command. This is to get around a weird and pesky error that will tell me that `score` is not a numeric variable when it clearly is. Doing it this way does not change anything meaningful for your output, but the error crops up a decent amount and this is the way around it.

So we will set `wid` to ID, `dv` to score, and then use the `within` input and set it to `noise.level`. Everything else is pretty much the same as before in the between-subjects ANOVA lab. This time, I'm setting `detailed` to TRUE to output a little more information and not returning an `aov` object. I'm not returning an `aov` object here because I'm not doing anything with the residuals or performing a post-hoc test right now.

```{r ez}
noise.ez <- ezANOVA(data = data.frame(noise.g),
                    wid = ID,
                    dv = score,
                    within = noise.level,
                    type = 3,
                    detailed = TRUE,
                    return_aov = TRUE)

print(noise.ez)
```

So if we look at the table, we will see essentially the same things we saw from the `car` `Anova` output. We have the overall results from the ANOVA at the top with sums-of-squares, degrees of freedom, F, and p-values. You also have effect size outputs which we will talk about in a later lab. The results table is exactly the same as from `Anova` which is reassuring.

The test for sphericity is also performed here and comes out non-significant just like before.

We also have the sphericity corrections below that. This time, both GG and HF corrections are in the same row. The Greenhouse-Geisser $\epsilon$ value is `GGe` in the table, same for `HFe`. All the results are the same, which is great. How you choose to perform the ANOVA is up to you, but for this case, ezANOVA is easier in my opinion.

## Reporting Results from One-Way Repeated-Measures ANOVA
Using the example data from this lab, we would write the results like so:

"Mauchly's test indicated that the assumption of sphericity had been preserved for the main effect of noise level, $\chi^2$(2) = 0.7133, p = 0.308. There was a significant main effect of noise level on score, F(3,57) = 482.81. Post-hoc tests revealed that....."


If sphericity had been violated, you would have multiplied the numerator and denominator degrees of freedom by either the GG or HF $\epsilon$ value in your report, replaced the p-value with the corrected p-value next to the test you ran, and added a sentence saying such: 

"Therefore degrees of freedom were corrected using the <insert GG or HF here, depending on which you used> estimates of sphericity ($\epsilon$ = \<insert GGe value\> for the main effect of \<insert name of your main effect\>)."