---
title: "ANCOVA"
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)
```



## Hot Dogs

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

data(hotdog)
```

-   Are hot dogs made with some types of meat healthier than others?
-   We will consider both calorie content and sodium levels for beef, poultry, and meat (mostly pork and beef) hot dogs
- Based on an example from text by Heiberger and Holland (2015)
- `hotdog` dataset from `HH` package

---

- Hot dogs made from poultry are often lower in calories
- Do manufacturers add sodium to enhance flavor to make up for the lower fat content?
- We might start to approach this question by comparing sodium content between hot dog types using a standard one-way ANOVA

## Sodium Levels

::: panel-tabset
### Ridge Plot

```{r}
#| include: true
#| echo: false
#| fig-cap: "Ridge plot showing distributions of sodium content in different types of hot dogs."

library(ggridges)
hotdog |>
  ggplot(aes(x = Sodium, y = Type, fill = Type)) +
  geom_density_ridges() +
  theme_minimal()
```

### Summary Statistics

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

hotdog |>
  group_by(Type) |>
  summarize(n = n(), mean = mean(Sodium), sd = sd(Sodium))
```

| Type   |  n  | mean |  sd  |
|---------|:---:|:----:|:----:|
| Beef   | 20  | 401 | 102 |
| Meat  | 17 | 419 | 93.9 |
| Poultry   | 17  | 459 | 84.7 |

:::

## One-Way ANOVA

Statistical model
$$y_{ij}=\mu + \alpha_i + \varepsilon_{ij}$$
ANOVA table

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

hds_lm <- lm(Sodium ~ Type, data = hotdog)

hds_lm |>
  anova() |>
  tidy()
```

---

- Based on this analysis we do not have convincing evidence that at least one of the mean sodium levels is different
- However, we have not accounted for the dependence of sodium on calorie content
- Next we will conduct an analysis that accounts for this covariate

## Analysis of Covariance (ANCOVA)

- Statistical model for the $j$th observation in the $i$th group:
$$y_{ij}=\mu + \alpha_i + \beta(X_{ij} - \bar{\bar{X}})+\varepsilon_{ij}$$
- Like before, $\mu$ is the overall population mean, $\alpha_i$ is the differential effect of group $i$, $\beta$ is the slope
- $\bar{\bar{X}}$ is the overall mean of the $X_{ij}$s
- Each group has a different intercept ($\mu+\alpha_i$)
- All groups have a common slope $\beta$

## Linear model

- We can fit a linear model as we have in the past
- One categorical predictor, one numeric predictor

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

hdsc_lm <- lm(Sodium ~ Calories + Type, data = hotdog)

hdsc_lm |>
  tidy()
```

---

- The linear model that R fit to the data is
$$\begin{array}{rcl}\widehat{Sodium} &=& -113 + 3.28\times Calories\\ && +11\times TypeMeat + 183\times TypePoultry\end{array}$$
- We can recast the standard ANCOVA model in a similar form
$$y_{ij} = (\mu-\beta\bar{\bar{X}}+\alpha_1) + (\alpha_i - \alpha_1) + \beta X_{ij} + \varepsilon_{ij}$$
- By identifying terms in this model with the regression output, we can estimate the coefficients in the standard model

---

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

hotsum <- hotdog |> 
  summarize(SodMean = mean(Sodium),
          CalMean = mean(Calories))

b2 <- hdsc_lm$coefficients["TypeMeat"]
b3 <- hdsc_lm$coefficients["TypePoultry"]
int <- hdsc_lm$coefficients["(Intercept)"]
b <- hdsc_lm$coefficients["Calories"]
xbar <- hotsum$CalMean
bxbar <- b*xbar

m <- (3*int + 3*bxbar + b2 + b3)/3
a1 <- int + bxbar - m
a2 <-  b2 + a1
a3 <-  b3 + a1
```

Prediction model in standard form
$$\widehat{Sodium} = 428+3.28\times (Calories-145) + \left\{\begin{array}{ll}-65 & \text{if } Beef\\
-53 & \text{if } Meat \\ 118 & \text{if } Poultry\end{array}\right.$$

---

```{r}
#| include: true
#| echo: false
#| fig-cap: "Sodium vs Calories, faceted by Type with linear model"

regline <- tibble(slope = rep(3.28, 3), intercept = 428-3.28*145 + c(-65, -53, 118), Type = c("Beef", "Meat", "Poultry"))

hotdog |>
  ggplot(aes(Calories, Sodium, color = Type)) +
  geom_point() +
  geom_abline(data = regline, aes(slope = slope, intercept = intercept)) +
  facet_wrap(~Type, nrow = 1) +
  theme_bw()
```

## Hypotheses

- We will test the hypotheses
  - $H_0: \alpha_1=\alpha_2=\alpha_3=0$
  - $H_A:$ at least one alpha is different
- However, this time our analysis (ANCOVA) will take into account the relationship between `Sodium` and `Calories`

## ANOVA table

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

hdsc_lm |>
  anova() |>
  tidy()
```

## A different conclusion

- When we take into the covariate (`Calories`) into account, we come to a different conclusion
- We reject the null hypothesis. There is an association between sodium and hotdog type
- The ANCOVA compared the intercepts of the three lines
- We found that the vertical distance between the lines is significantly different from 0

## Adjusting for Calories

We can adjust `Sodium` for `Calories` by subtracting $b(X_{ij}-\bar{\bar{X}})$ from each $y_{ij}$, where $b$ is the estimate of the slope

```{r}
#| include: true
#| echo: false
#| fig-cap: "Sodium Adjusted for Calories vs Calories, faceted by Type with adjusted linear model"

regline <- tibble(intercept = 428 + c(-65, -53, 118), Type = c("Beef", "Meat", "Poultry"))

hotdog |>
  mutate(Sodium.Adjusted = Sodium - b*Calories + 3.28*145) |>
  ggplot(aes(Calories, Sodium.Adjusted, color = Type)) +
  geom_point() +
  geom_hline(data = regline, aes(yintercept = intercept)) +
  facet_wrap(~Type, nrow = 1) +
  theme_bw()
```

## Sequential sums of squares

 - R computes sums of square sequentially by default
 - First, the sums of squares for calories is calculated (as a regression sum of squares) $$SS_{Calories}=\sum_{i=1}^n(\hat{y}_i-\bar{y})^2$$
 - $\hat{y}$ is based on a model that does not account for hot dog type
 
---

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

lm(Sodium ~ Calories + Type, data = hotdog) |>
  anova() |>
  tidy()
```

Compare to $SS_{Calories}$ without `Type`

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

lm(Sodium ~ Calories, data = hotdog) |>
  anova() |>
  tidy()
```

---

- Next, the sums of squares for hot dog type is calculated, accounting for calories
$$SS_{Type}=\left(\sum_{i=1}^a\sum_{j=1}^{n_i}(\hat{y}_{ij}-\bar{\bar{y}})^2\right)-SS_{Calories}$$
- Here, the prediction $\hat{y}_{ij}$ uses the full model: a different intercept for each type of hot dog (but same slope)
- This is the sum of squares that is accounted for by the full model that is not accounted for by calories alone

---

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

Sodium_mean <- hotdog |>
  summarize(mean = mean(Sodium)) |>
  pull()

hdsc_lm |>
  augment() |>
  mutate(ssij = (.fitted - Sodium_mean)^2) |>
  summarize(SS = sum(ssij)) |>
  pull()

# subtract SS_{Calories} = 106269.7
# to get SS_{Type}
```
Compare $SS_{Type}$ accounting for `Calories`

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

lm(Sodium ~ Calories + Type, data = hotdog) |>
  anova() |>
  tidy()
```

To $SS_{Type}$ without accounting for `Calories`

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

lm(Sodium ~ Type, data = hotdog) |>
  anova() |>
  tidy()
```

## Sequential Sums of Squares: Order

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

lm(Sodium ~ Calories + Type, data = hotdog) |>
  anova() |>
  tidy()
```

Compare to `Type` first

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

lm(Sodium ~ Type + Calories, data = hotdog) |>
  anova() |>
  tidy()
```

---

If there is one factor of interest (`Type`), but we want to account for another variable (`Calories`), the factor of interest should enter the model *last*