---
title: "Analysis of variance"
editor: 
  markdown: 
    wrap: 72
---

## Packages

```{r}
library(tidyverse)
library(smmr)
library(PMCMRplus)
```

## Jumping rats

-   Link between exercise and healthy bones (many studies).
-   Exercise stresses bones and causes them to get stronger.
-   Study (Purdue): effect of jumping on bone density of growing rats.
-   30 rats, randomly assigned to 1 of 3 treatments:
    -   No jumping (control)
    -   Low-jump treatment (30 cm)
    -   High-jump treatment (60 cm)
-   8 weeks, 10 jumps/day, 5 days/week.
-   Bone density of rats (mg/cm$^3$) measured at end.

## Jumping rats 2/2

-   See whether larger amount of exercise (jumping) went with higher
    bone density.
-   Random assignment: rats in each group similar in all important ways.
-   So entitled to draw conclusions about cause and effect.

## Reading the data

Values separated by spaces:

```{r inference-5-R-9}
my_url <- "http://ritsokiguess.site/datafiles/jumping.txt"
rats <- read_delim(my_url," ")
```

## The data

```{r inference-5-R-10}
# rats %>% slice_sample(n=10)
rats
```

## Boxplots

```{r inference-5-R-11, fig.height=5}
ggplot(rats, aes(y=density, x=group)) + geom_boxplot()
```

## Or, arranging groups in data (logical) order

```{r inference-5-R-12, fig.height=5}
ggplot(rats, aes(y=density, x=fct_inorder(group))) +
  geom_boxplot()
```

```{r}
rats %>% mutate(group = fct_inorder(group)) -> rats
ggplot(rats, aes(y=density, x=group)) +
  geom_boxplot()

```

## Analysis of Variance

-   Comparing > 2 groups of independent observations (each rat only
    does one amount of jumping).
-   Standard procedure: analysis of variance (ANOVA).
-   Null hypothesis: all groups have same mean.
-   Alternative: "not all means the same", at least one is different
    from others.

## Testing: ANOVA in R

```{r inference-5-R-13}
rats.aov <- aov(density ~ group, data=rats)
summary(rats.aov)
```

-   Usual ANOVA table, small P-value: significant result.
-   Conclude that the mean bone densities are not all equal.
-   Reject null, but not very useful finding.

## Which groups are different from which?

-   ANOVA really only answers half our questions: it says "there are
    differences", but doesn't tell us which groups different.
-   One possibility (not the best): compare all possible pairs of
    groups, via two-sample t.
-   First pick out each group:

```{r inference-5-R-14}
rats %>% filter(group=="Control") -> controls
rats %>% filter(group=="Lowjump") -> lows
rats %>% filter(group=="Highjump") -> highs
```

## Control vs. low

```{r inference-5-R-15}
t.test(controls$density, lows$density)
```

No sig. difference here.

## Control vs. high

```{r inference-5-R-16}
t.test(controls$density, highs$density)
```

These are different.

## Low vs. high

```{r inference-5-R-17}
t.test(lows$density, highs$density)
```

These are different too.

## But...

-   We just did 3 tests instead of 1.
-   So we have given ourselves 3 chances to reject $H_0:$ all means
    equal, instead of 1.
-   Thus $\alpha$ for this combined test is not 0.05.

## John W. Tukey

::: columns
::: {.column width="40%"}
![](John_Tukey.jpg){width="200"}
:::

::: {.column width="60%"}
-   American statistician, 1915--2000
-   Big fan of exploratory data analysis
-   Popularized boxplot
-   Invented "honestly significant differences"
-   Invented jackknife estimation
-   Coined computing term "bit"
-   Co-inventor of Fast Fourier Transform
:::
:::

## Honestly Significant Differences

-   Compare several groups with one test, telling you which groups
    differ from which.
-   Idea: if all population means equal, find distribution of highest
    sample mean minus lowest sample mean.
-   Any means unusually different compared to that declared
    significantly different.

## Tukey on rat data

```{r inference-5-R-18, echo=F}
width <- getOption("width")
options(width = 60)
```

```{r inference-5-R-19}
rats.aov <- aov(density~group, data = rats)
TukeyHSD(rats.aov)
```

```{r inference-5-R-20, echo=F}
options(width=width)
```

-   Again conclude that bone density for highjump group significantly
    higher than for other two groups.

## Why Tukey's procedure better than all t-tests

Look at P-values for the two tests:

```         
Comparison        Tukey    t-tests
----------------------------------
Highjump-Control 0.0016     0.0021
Lowjump-Control  0.4744     0.2977
Lowjump-Highjump 0.0298     0.0045
```

-   Tukey P-values (mostly) higher.
-   Proper adjustment for doing three t-tests at once, not just one in
    isolation.

## Checking assumptions

```{r inference-5-R-21}
#| fig.height = 4
ggplot(rats,aes(y = density, x = fct_inorder(group)))+
  geom_boxplot()
```

Assumptions:

-   Normally distributed data within each group
-   with equal group SDs.

## Normal quantile plots by group

```{r inference-5-R-22, fig.height=3.5}
ggplot(rats, aes(sample = density)) + stat_qq() + 
  stat_qq_line() + facet_wrap( ~ group)
```

## The assumptions

-   Normally-distributed data within each group
-   Equal group SDs.
-   These are shaky here because:
    -   control group has outliers
    -   highjump group appears to have less spread than others.
-   Possible remedies (in general):
    -   Transformation of response (usually works best when SD increases
        with mean)
    -   If normality OK but equal spreads not, can use Welch ANOVA.
        (Regular ANOVA like pooled t-test; Welch ANOVA like
        Welch-Satterthwaite t-test.)
    -   Can also use Mood's Median Test (see over). This works for any
        number of groups.

## Mood's median for multiple groups

-   Find median of all bone densities, regardless of group

-   Count up how many observations in each group above or below overall
    median
    
- Test association between group and above/below

-   Mood's `median_test` (over).

## Mood's median test here

```{r inference-5-R-27}
median_test(rats, density, group)
```

## Comments

-   No doubt that medians differ between groups (not all same).
-   This test is equivalent of $F$-test, not of Tukey.
-   To determine which groups differ from which, can compare all
    possible pairs of groups via (2-sample) Mood's median tests, then
    adjust P-values by multiplying by number of 2-sample Mood tests done
    (Bonferroni):

```{r inference-5-R-28}
pairwise_median_test(rats, density, group)
```

-   Now, lowjump-highjump difference no longer significant.

## Welch ANOVA

-   For these data, Mood's median test probably best because we doubt
    both normality and equal spreads.
-   When normality OK but spreads differ, Welch ANOVA way to go.
-   Welch ANOVA done by `oneway.test` as shown (for illustration):

```{r inference-5-R-29}
oneway.test(density~group, data=rats)
```

-   P-value very similar, as expected.
-   Appropriate Tukey-equivalent here called Games-Howell.

## Games-Howell

-   Lives in package `PMCMRplus`. Install first.

```{r games-howell, warning=F}
gamesHowellTest(density ~ group, data = rats)
# gamesHowellTest(density ~ factor(group), data = rats)
```

The line below appears to be no longer true:

Careful: explanatory must be `factor` (so commented-out line does not work).

## Deciding which test to do

For two or more samples:

![](testflow.png)
