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

# Complete Block Design

```{r}
tree <- c(rep("tree1",3), rep("tree2",3), 
          rep("tree3", 3),rep("tree4", 3), rep("tree5", 3))
treatment <- rep(c("water", "spores", "oil"), 5)
bugs_change <- c(-9, -4, 4, 18, 29, 29, 10, 4, 14, 9, -2, 14, -6, 11, 7)

mealybugs <- data.frame(tree, treatment, bugs_change)
```

## Formal ANOVA

```{r}
mod <- lm(bugs_change ~ treatment + tree, data = mealybugs)

anova(mod)
```

## Informal Analysis

```{r}
#testing for S assumption
mealybugs %>%
  group_by(treatment) %>%
  summarise(sd = sd(bugs_change))
```

```{r}
mealybugs_wide <- mealybugs %>%
  spread(treatment, bugs_change)

ggplot(mealybugs_wide, aes(x = spores, y = oil)) +
  geom_point() +
  geom_abline(intercept = 8, slope = 1) +
  geom_smooth(method = "lm", se = 0)
```

```{r}
ggplot(mealybugs_wide, aes(x = spores, y = water)) +
  geom_point() +
  geom_abline(intercept = -1, slope = 1) +
  geom_smooth(method = "lm", se = 0)
```

```{r}
ggplot(mealybugs_wide, aes(x = oil, y = water)) +
  geom_point() +
  geom_abline(intercept = -7, slope = 1) +
  geom_smooth(method = "lm", se = 0)
```

# Split Plot/Repeated Measures

To illustrate the analysis of the split plot/repeated measures design, we will use our familiar hamster data. Instead of entering in the data, we can load it from the URL below.

```{r}
hamster <- read_csv("https://randilgarcia.github.io/sds290interterm21/KellysHamsters.csv")
```

We can use the `glimpse()` function from the `dplyr` package to get a first look at the data and the variables, as well as the variable types.

```{r}
glimpse(hamster)
```

All of our factor variables are saved as `int` type variables, so first we need to identify factor variables as categorical. We will use the `factor()` function in combination with the `mutate()` function to accomplish this.

```{r}
hamster <- hamster %>%
  mutate(day_length = factor(day_length, levels=c(0,1),labels=c("short","long")),
         organ = factor(organ, levels=c(0,1), labels=c("brain","heart")),
         interact = factor(interact, levels=c(1,2,3,4), labels=c("LH","SH","LB","SB")),
         id = factor(id)) 

glimpse(hamster)
```

## Informal analyses

### Parallel dot graph

We can then move on to the informal analyses stage.

```{r}
ggplot(hamster, aes(x = day_length, y = conc, color = organ)) +
  geom_jitter(height = 0, width = 0.03, alpha = .7)
```

Check group means and standard deviations. First for organ. Note that here we are using the `group_by()` and `summarize()` function from `dplyr` instead of the `favstats()` function from `mosaic`. There no real reason to use one over the other, it's just your preference.

```{r}
hamster %>%
  group_by(organ) %>%
  summarise(mean = mean(conc),
            sd = sd(conc))
```

Then for day length.

```{r}
hamster %>%
  group_by(day_length) %>%
  summarise(mean = mean(conc),
            sd = sd(conc))
```

We can also do this for interaction cell means.

```{r}
hamster %>%
  group_by(interact) %>%
  summarise(mean = mean(conc),
            sd = sd(conc))
  
hamster %>%
  group_by(organ, day_length) %>%
  summarise(mean = mean(conc),
            sd = sd(conc))
```

### Transformation

The variances are not equal, so we can try a transformation.

```{r}
hamster <- hamster %>%
  mutate(lconc = log10(conc)*100) 

hamster %>%
  group_by(interact) %>%
  summarise(lmean = mean(lconc),
            lsd = sd(lconc))
```

### Parallel dot graphs

Draw plots using transformed data. It looks better

```{r}
ggplot(hamster, aes(x = day_length, y = lconc, color = organ)) +
  geom_jitter(height = 0, width = 0.03, alpha = .7)
```

Try a side-by-side boxplot too if you'd like!

### Interaction Plots

Next, because the SP/RM design also has crossed factors, we should make an interaction graph.

```{r}
ggplot(hamster, aes(x = organ, y = lconc, 
                    group = day_length, 
                    linetype = day_length, 
                    shape = day_length,
                    color = day_length)) +
  geom_jitter(height = 0, width = 0.03, alpha = .7) +
  geom_smooth(method = "lm", se = 0)
```

Try flipping the aesthetic mappings and having day length on the x-axis and organ as color.

### Scatter plots

Whenever we have a within block factor we can check the additivity assuption (A) by making scatterplot(s) and looking for the `x=y` regression line. That is, the `45 degree`, or `slope = 1` line. We will need to first restructure data for scatter plots using `spread()` from `tidyr`.

```{r}
hamster %>%
  select(-conc, -interact) %>%
  spread(organ, lconc) %>%
  ggplot(aes(x = heart, y = brain, shape = day_length)) +
  geom_point() 
```

## Formal Modeling with the ANOVA

Split-plot ANOVA. Things get a little weird with the SP/RM ANOVA because we need to incorporate the error term defined by the blocks nested within the between block factor. The `aov()` function wanted to know about your blocking factor by adding `+ Error(id)` in the equation. Here, `id` is the hamster ID variable in your data set.

```{r}
mod <- aov(lconc ~ day_length*organ + Error(id), data = hamster)
summary(mod)
qplot(x = residuals(mod), bins = 6)
qplot(sample = residuals(mod)) + stat_qq_line()
```

Or we can use a linear model, but we'll need the `lmer()` function from the `lme4` package. The `lmer()` function wants the blocking factor added as `+ (1|id)`.

```{r}
library(lme4)
library(lmerTest) #using this package in addition to lme4 gives us more information in the ANOVA table

mod <- lmer(lconc ~ day_length*organ + (1|id), data = hamster)

anova(mod) #note that the F-ratios are based on df calculated with a different method (Satterthwaite) than we've been learning in class. Please use the aov() function for the ANOVA table
summary(mod) #use the lmerTest package and lmer() function if you desire test of differences between specific conditions. Note that the df are calculated with Satterthwaite
confint(mod)
```

Note that `lmer()` doesn't give p-values.
