---
title: "Piglets Data Two-Way Factorial"
author: "Randi Garcia"
---

```{r global-options, include=FALSE}
knitr::opts_chunk$set(eval = FALSE)
library(tidyverse)
library(dae)
library(Stat2Data)
```

# Explore the Data

```{r}
data("PigFeed")
```

First we get the cell counts.

```{r}
PigFeed %>%
  group_by(Antibiotic, B12) %>%
  count()
```

Then we look at it visually with boxplots.

```{r}
ggplot(PigFeed, aes(x = interaction(B12, Antibiotic), y = WgtGain)) + 
  geom_boxplot()
```

Or with `fill`.

```{r}
ggplot(PigFeed, aes(x = B12, y = WgtGain, fill = Antibiotic)) +
  geom_boxplot()
```

Create a parallel dot plot for `PigFeed` data

```{r}
ggplot(PigFeed, aes(x = B12, y = WgtGain, color = Antibiotic)) +
  geom_point()
```

# Check for Equal Variances

```{r}
ds <- PigFeed %>%
  group_by(Antibiotic, B12) %>%
  summarise(n = n(),
            mean = mean(WgtGain),
            sd = sd(WgtGain))
ds

ds %>% 
  ungroup() %>%
  summarise(max(sd)/min(sd))
```

Note variances unequal. Try log transformation.

```{r}
PigFeed <- PigFeed %>%
  mutate(logWgtGain = log(WgtGain + 1)) 

ggplot(PigFeed, aes(x = interaction(B12, Antibiotic), y = logWgtGain)) + 
  geom_boxplot()
```

Doesn't help much!

```{r}
PigFeed %>%
  group_by(Antibiotic, B12) %>%
  summarise(n = n(),
            mean = mean(logWgtGain),
            sd = sd(logWgtGain)) %>% 
  ungroup() %>%
  summarise(max(sd)/min(sd))
```

# Formal Analysis

Run two-way ANOVA without the interaction term using the original data.

```{r}
mod1 <- lm(WgtGain ~ Antibiotic + B12 + Antibiotic*B12, data = PigFeed)

anova(mod1)
```

## Interactions

Check for interactions: Create an interaction plot.

```{r}
ggplot(PigFeed, aes(x = B12, y = WgtGain, 
                    group = Antibiotic, 
                    color = Antibiotic)) +
  geom_point() +
  geom_smooth(method = "lm", se = 0)
```

Conduct a two-way ANOVA

```{r}
mod2 <- lm(WgtGain ~ Antibiotic * B12, data = PigFeed)

anova(mod2)
```

# Residual Analysis

Check conditions on the residuals

```{r}
qplot(x = mod2$residuals, bins = 6)
plot(mod2, which = 2)
plot(mod2, which = 1)
```

# Confidence intervals and effect sizes

```{r}
SD = sqrt(36.25) 
t = qt(.975, 8)

#Conditional averages
mean_B12_antiNO <- 22 - 19

mean_B12_antiYES <- 54 - 3

LL_B12_antiNO = mean_B12_antiNO - SD*t*sqrt(1/3+1/3)
UL_B12_antiNO = mean_B12_antiNO + SD*t*sqrt(1/3+1/3)

LL_B12_antiYES = mean_B12_antiYES - SD*t*sqrt(1/3+1/3)
UL_B12_antiYES = mean_B12_antiYES + SD*t*sqrt(1/3+1/3)

LL_B12_antiNO
UL_B12_antiNO

LL_B12_antiYES
UL_B12_antiYES
```

```{r}
mean_B12_antiNO/SD
mean_B12_antiYES/SD
```