---
title: 'Hybridization and invasiveness in Eurasian watermilfoil (*Myriophyllum spicatum*):
  Is prioritizing hybrids in management justified?'
author: "Hannah K. Hoff and Ryan A. Thum"
date: "07/16/2021"
output:
  pdf_document: default
  word_document: default
---

# 1. Load the data into R.
```{r}
library(readr)
ExptData <- read.table("2,4-D_Assay_Data_Full.txt", header = TRUE, sep = "\t")
print.data.frame(ExptData)
```

```{r, include = FALSE}
head(ExptData)
tail(ExptData)
Genotype <- as.factor(ExptData$Genotype)
Tank <- as.factor(ExptData$Tank)
Treatment <- as.factor(ExptData$Treatment)
Taxon <- as.factor(ExptData$Taxon)
Dry.Weight <- ExptData$Dry.Weight
Trial <- as.factor(ExptData$Trial)
Rep <- as.factor(ExptData$Rep)
Gen_Code <- as.factor(ExptData$Gen_Code)
Trt_Code <- as.factor(ExptData$Trt_Code)
library(car)
library(lmerTest)
##interaction plot E-H##
interaction.plot(x.factor = Treatment, trace.factor = Taxon, response = Dry.Weight,ylab = "Average Dry Weight (g)", xlab = "Treatment", col = 3:10, lty = 1, pch = c(16, 24, 22), type = "l")
interaction.plot(x.factor = Treatment, trace.factor = Genotype, response = Dry.Weight,ylab = "Average Dry Weight (g)", xlab = "Treatment", col = 3:10, lty = 1, pch = c(16, 24, 22), type = "l")
```

# 2. Determine best fit linear mixed effects model.
```{r}
Expt.lmer1 <- lmer(Dry.Weight ~ Genotype %in% Taxon * Taxon * Treatment  + (1|Trial),
                   data = ExptData)

Expt.lmer2 <- lmer(Dry.Weight ~ Genotype %in% Taxon * Taxon * Treatment + (1 | Trial:Tank) +
                     (1 | Tank), data = ExptData)

Expt.lmer3 <- lmer(Dry.Weight ~ Genotype %in% Taxon * Taxon * Treatment + (1 | Trial:Tank:Rep) +
                     (1|Tank:Rep) + (1|Rep), data = ExptData)

anova(Expt.lmer1, Expt.lmer2, Expt.lmer3)
```

```{r, include = FALSE}
qqnorm(resid(Expt.lmer1))
qqline(resid(Expt.lmer1)) 
Plot.Model.Variances <- plot(Expt.lmer1) 
Plot.Model.Variances 

logDW <- log10(Dry.Weight)
Expt.lmer.log <- lme4::lmer(logDW ~ Genotype %in% Taxon * Taxon * Treatment  + (1|Trial),
                     data = ExptData)

Plot.Model.Variances <- plot(Expt.lmer.log)
Plot.Model.Variances

qqnorm(resid(Expt.lmer.log))
qqline(resid(Expt.lmer.log))
shapiro.test(logDW)
Plot.Model.Variances.log <- plot(Expt.lmer.log) 
Plot.Model.Variances.log

library(afex)
options(contrasts = c("contr.sum","contr.poly"))
Expt.Anov.log <- Anova(Expt.lmer.log, type = 2, test.statistic = "F")
Expt.Anov.log <- Anova(Expt.lmer.log, type = 2, REML = TRUE, test.statistic = "F")
Expt.Anov.log

library(lmerTest)
Expt.lmer.log.1 <- lmer(logDW ~ Genotype %in% Taxon * Taxon * Treatment  + (1|Trial),
     data = ExptData)
```

# 3. ANOVA of the best fit linear mixed effects model (overall data were log-transformed to meet ANOVA assumptions; see Materials and Methods). 
```{r}
library(lmerTest)
log.anova <- anova(Expt.lmer.log.1, type = 2, ddf = "Satterthwaite")
log.anova
```

# 4. Calculating pairwise least-square means and contrasts between Eurasian and hybrid watermilfoil in control and treatment. 
```{r}
library(emmeans)
Expt.contrast <- emmeans(Expt.lmer.log.1, ~ Taxon | Treatment)
pairs(Expt.contrast)
```

# 5. Calculating pairwise comparisons between genotypes in control and treatment. 
```{r}
library(devtools)
library(emmeans)
Expt.tukey = emmeans(Expt.lmer.log, pairwise~ Genotype * Treatment, adjust="tukey")
Expt.tukey
```
