---
title: "Two-Way ANOVA"
subtitle: |
  | Additional Topic 
  | Math 219
author: "Yurk"
format: 
  revealjs:
    theme: beige
editor: source
---

```{r}
#| include: false
#| echo: false

library(tidyverse)
library(infer)
library(openintro)
library(broom)
library(HH)
library(GLMsData)
```

## Experiments with two factors

- We have considered experiments with one treatment factor and one blocking factor 
- Now we will consider experiments with two treatment factors
- The approach generalizes to more than two factors as well

## Experimental Designs

- Randomized complete block design (RCBD): 1 treatment factor with $t$ levels, 1 blocking factor with $b$ levels, exactly $t$ experimental units in each block (one for each treatment level)
- **Factorial design**: 2 or more factors, 1 or more observations per cell (**replication**)
- A factorial design is called **balanced** if there is the same number of observations in each cell
- Our focus will be on balanced factorial designs

## Iron Content and Pot Type

```{r}
#| include: false
#| echo: false

data(toothbrush)

toothbrush <- toothbrush |>
  mutate(plaque_rem = Before - After)

lm(plaque_rem ~ Sex + Toothbrush, data = toothbrush) |>
  anova() |>
  tidy()

lm(plaque_rem ~ Toothbrush + Sex, data = toothbrush) |>
  anova() |>
  tidy()

data(flowers)

lm(Flowers ~ Light*Timing, data = flowers) |>
  anova() |>
  tidy()

data(ironpot)

lm(Iron ~ Pot*Food, data = ironpot) |>
  anova() |>
  tidy()
```

- Anemia, caused by iron deficiency, is a common form of malnutrition in developing countries
- Does iron content in cooked food depend on the type of pot?
- Do the results depend on the type of food?
- Data from Adish et al. (1999) [^1]
- `ironpot` dataset from `HH` package

[^1]: Adish, A.A., Esrey, S.A., Gyorkos, T.W., Jean-Baptiste, J., Rojhani, A., 1999. Effect of consumption of food cooked in iron pots on iron status and growth of young children: a randomized trial, Lancet 353, Pp. 712-716.

---

- Response: `Iron` (mg per 100 g of food)
- Two Treatment factors:

  - `Pot` with 3 levels ($a=3$): "Aluminum", "Clay", "Iron"
  - `Food` with 3 levels ($b=3$): "meat", "legumes", "vegetables"
  
- Meat, legume, and vegetable dishes cooked according to recipes from Ethiopian Nutritional Institute
- Each dish was cooked 4 times in each type of pot (4 experimental units per cell, 9 cells)

## Additive Statistical Model

**Additive model** for the $k$th observation at the $i$th level of factor $A$ and the $j$the level of factor $B$

$$y_{ijk}=\mu+\alpha_i+\beta_j+\varepsilon_{ijk}$$

- $\mu$ is the overall/grand mean
- $\alpha_i$ is the differential effect of the $i$th level of treatment $A$
- $\beta_j$ is the differential effect of the $j$th level of treatment $B$
- $\varepsilon_{ijk}\sim N(0,\sigma^2)$ represents the error

## Hypothesis Tests

- There are two null-hypotheses of interest
- For type of pot:

  - $H_0: \alpha_1 = \alpha_2 = \alpha_3 = 0$
  - $H_A:$ At least one $\alpha_i$ is different
  
- For type of food:

  - $H_0: \beta_1 = \beta_2 = \beta_3 = 0$
  - $H_A:$ At least one $\beta_i$ is different

## EDA

::: panel-tabset
### Pot Type

```{r}
#| include: TRUE
#| echo: FALSE
#| fig-cap: "Iron content in food cooked in different types of pots."

ironpot |>
  ggplot(aes(Pot, Iron)) +
  geom_boxplot() +
  theme_minimal()
```

### Food Type

```{r}
#| include: TRUE
#| echo: FALSE
#| fig-cap: "Iron content of different types of cooked food."

ironpot |>
  ggplot(aes(Food, Iron)) +
  geom_boxplot() +
  theme_minimal()
```

### Both

```{r}
#| include: TRUE
#| echo: FALSE
#| fig-cap: "Iron content of different types of food cooked in different types of pots."

ironpot |>
  ggplot(aes(Food, Iron, color = Pot)) +
  geom_boxplot() +
  theme_minimal()
```
:::

## ANOVA Table (Additive Model)

```{r}
#| include: true
#| echo: true

lm(Iron ~ Pot + Food, data = ironpot) |>
  anova() |>
  tidy()
```


## Does order matter?

```{r}
#| include: true
#| echo: true

lm(Iron ~ Pot + Food, data = ironpot) |>
  anova() |>
  tidy()
```

```{r}
#| include: true
#| echo: true

lm(Iron ~ Food + Pot, data = ironpot) |>
  anova() |>
  tidy()
```

---

- Due to the experimental design, there is no association between the treatment factors
- There are the same number of replicates of each dish cooked with each type of pot
- As a result, the order that the treatment variables enter the model does not matter

## Conclusions based on additive model

- Using the additive model, we would conclude that there is convincing evidence of an effect for both food type and pot type
- However, the additive model assumes that the effect of pot type is the same for each food type
- This does not appear to be the case

---

```{r}
#| include: TRUE
#| echo: FALSE
#| fig-cap: "Iron content of different types of food cooked in different types of pots."

ironpot |>
  ggplot(aes(Food, Iron, color = Pot)) +
  geom_boxplot() +
  theme_minimal()
```

## Statistical Model with Interaction

Statistical model with **interaction**

$$y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\varepsilon_{ijk}$$

- $\mu$ is the overall/grand mean
- $\alpha_i$ is the differential effect of the $i$th level of treatment $A$
- $\beta_j$ is the differential effect of the $j$th level of treatment $B$
- $(\alpha\beta)_{ij}$ models the possible interaction between treatment $A$ and $B$
- $\varepsilon_{ijk}\sim N(0,\sigma^2)$ represents the error

---

- For example, if cooking a dish in an iron pot ($i=3$) raises the iron level more for meat ($j=2$) than it does for legumes or vegetables, then we would expect $(\alpha\beta)_{3,2}>0$
- Treatments $A$ and $B$ **interact** if the difference in response between two levels of $A$ depends on the level of $B$
- If an interaction is present, it is inappropriate to use an additive model

## Properties of Interaction Model

- The interaction coefficients satisfy $$\sum_{i=1}^a(\alpha\beta)_{ij}=\sum_{j=1}^b(\alpha\beta)_{ij}=0$$
- The model can be written more simply as
$$y_{ijk}=\mu_{ij}+\varepsilon_{ijk}$$
where $\mu_{ij}$ is the mean for level $i$ of treatment $A$ and level $j$ of treatment $B$

## Testing for an Interaction

- Test for an interaction using the hypotheses


  - $H_0: (\alpha\beta)_{ij} = 0,$ for all $i=1\ldots a$, $j=1\ldots b$
  - $H_A:$ At least one $(\alpha\beta)_{ij}$ is different
  
- If interaction is present, we do not test for main effects
- If there is not a significant interaction, we can drop the interaction and use the additive model instead
- We can only test for an interaction if there is replication, otherwise there are not enough degrees of freedom

## ANOVA table with interaction


```{r}
#| include: true
#| echo: true

lm(Iron ~ Pot * Food, data = ironpot) |>
  anova() |>
  tidy()
```

ANOVA table key ($n$ = number of observations in each cell)

::: {style="font-size: 18px"}

| term | df | sumsq | meansq | statistic |
|--|:--:|:--:|:--:|:--:|
| *Treatment A* | $df_A=a-1$ | $SSA$ | $MSA=SSA/df_A$ | $F=MSA/MSE$ |
| *Treatment B* | $df_B=b-1$ | $SSB$ | $MSB=SSB/df_B$ | $F=MSB/MSE$ |
| *Interaction AB* | $df_{AB}=(a-1)(b-1)$ | $SSAB$ | $MSAB=SSAB/df_{AB}$ | $F=MSAB/MSE$ |
| Residuals (error) | $df_E=ab(n-1)$ | $SSE$ | $MSE=SSE/df_E$ |  |

:::

---

- There is a significant interaction between food type and plot type
- There is convincing evidence that the response of iron content to pot type depends on the type of food

## Interaction sum of squares

- The model predicts the response (iron content) for an observation with level $i$ for treatment $A$ and level $j$ for treatment $B$ using the corresponding cell sample mean
$$\widehat{y_{ijk}}=\bar{y}_{ij}$$
- The sum of squares for the interaction term is
$$SSAB = n\sum_{i=1}^{a}\sum_{j=1}^{b}\left(\bar{y}_{ij}-\bar{\bar{y}}\right)^2-SSA - SSB$$

## Cell Means

- When an interaction is present we report the cell means
- We can make a table of the cell means in R

```{r}
#| include: true
#| echo: true

ironpot |>
  group_by(Pot, Food) |>
  summarize(mean = mean(Iron)) |>
  pivot_wider(names_from = Pot, names_prefix = "Pot_",
              values_from = mean)
```

## Interaction plot

- An **interaction plot** shows how the cell means vary with respect to the levels of one of the treatments
- Separate line plots for each level of the other treatment

---

```{r}
#| include: TRUE
#| echo: FALSE
#| fig-cap: "Interaction plot for food type and pot type."

ironpot |>
  group_by(Pot, Food) |>
  summarize(mean = mean(Iron)) |>
  ggplot(aes(Pot, mean, color = Food, group = Food)) +
  geom_point() +
  geom_line() +
  labs(y = "Mean Iron Content") +
  theme_minimal()
```

## Interaction Effects

- **Interaction effects** compare cell means across levels of both factors
- For example, $$(\bar{y}_{12}-\bar{y}_{32})-(\bar{y}_{13}-\bar{y}_{33})=(2.06-4.68)-(1.23-2.79) = -1.06$$ is the interaction effect that compares the difference in cell means for aluminum ($i=1$) and iron ($i=3$) pots across the meat ($j=2$) and vegetable ($j=3$) food types

## Toothbrushes

- Is there a difference in the effectiveness of different types of toothbrushes in removing plaque?
- Do the results depend on sex of the toothbrusher?
- `tootbrush` dataset from `GLMsData` package

---

- Response: `plaque_rem` (difference in a plaque index after brushing - before brushing)
- Independent variables:

  - `Sex` with 2 levels ($a=2$): "F", "M"
  - `Toothbrush` with 2 levels ($b=2$): "Conventional", "Hugger"
  
---

- There is not a significant interaction
- We should consider an additive model instead

```{r}
#| include: true
#| echo: true

lm(plaque_rem ~ Sex * Toothbrush, data = toothbrush) |>
  anova() |>
  tidy()
```

---

- There is convincing evidence of an association between sex and the amount of plaque removed
- However, we are unable to reject the null hypothesis for the effect of toothbrush type

```{r}
#| include: true
#| echo: true

lm(plaque_rem ~ Sex + Toothbrush, data = toothbrush) |>
  anova() |>
  tidy()
```

## Main effects

- The **main effects** compare the marginal means for one factor
- For example, $\bar{y}_{1.}-\bar{y}_{2.}$ compares means for females ($i=1$) and males ($i=2$)
- And, $\bar{y}_{.1}-\bar{y}_{.2}$ compares means for conventional toothbrushes ($j=1$) and hugger toothbrushes ($j=2$)
- It is inappropriate to consider main effects when there is an interaction

---

- We can use the `TukeyHSD` function calculate the main effects, including confidence intervals
- This is also appropriate for follow-up tests comparing pairs of marginal means (when there are more than 2 levels)

---

```{r}
#| include: true
#| echo: true

aov(plaque_rem ~ Sex + Toothbrush, data = toothbrush) |>
  TukeyHSD()
```

