```{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()
```

```{r}
ggplot(mealybugs_wide, aes(x = oil, y = water)) +
  geom_point()
```

# Latin Square

```{r}
subject <- c(rep("1",3), rep("2",3), rep("3", 3))
treatment <- c("solution", "capsule", "tablet",
               "capsule", "tablet", "solution",
               "tablet", "solution", "capsule")
period <- rep(c("1", "2", "3"),3)
concentration <- c("A 1799", "C 1846", "B 2147", "C 2075", "B 1156", "A 1777",
                   "B 1396", "A 868", "C 2291")

bioequivalence <- data.frame(subject, treatment, period, concentration) %>%
  separate(concentration, c("group", "c_curve"))
```

## Informal ANOVA

```{r}
ggplot(bioequivalence, aes(x = treatment, y = c_curve)) +
  geom_point()
```

```{r}
sd <- favstats(c_curve ~ treatment, data = bioequivalence)[,8]

max(sd)/min(sd)
```

## Formal ANOVA

```{r}
ls_mod <- lm(c_curve ~ treatment + period + subject, data = bioequivalence)

anova(ls_mod)
```

```{r}
bioequivalence <- bioequivalence %>%
  mutate(fitted = fitted(ls_mod), 
         residuals = residuals(ls_mod))

ggplot(bioequivalence, aes(x = fitted, residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red")
```

# 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)
```

Or we can use a linear model, but we'll need the `lmer()` function from the `lme4` package. The `lmer()`
functon 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.

