install.packages('ggplot2')
install.packages('lme4') 
install.packages('lmerTest')
install.packages('emmeans')
install.packages('MCMCglmm)
install.packages('plotMCMC')
install.packages('tidybayes')

library(dplyr)
library(ggplot2)
library(lme4)
library(lmerTest)
library(emmeans)
library(corrplot)
library(ggfortify)
library(reshape2)
library(car)
library(MCMCglmm)
library(tidyverse)
library(reshape2)
library(plotMCMC)
library(tidybayes)

#issues are that:
# 1) coolers weren't selected randomly at each timepoint; the same coolers were used for TP=1-4 and then 
#    5-7; would have been better to randomly select but I don't think this would change the model we can fit
# 2) ideally would have sampled plants (in pairs) rather than tanks

#read in data 
rm(list = ls(all=TRUE))
TS1_dat <- read.csv("TimeSeries_Expt1.csv")

#make factors categorical 
TS1_dat$timepoint <- as.factor(TS1_dat$timepoint)
TS1_dat$coolerno <- as.factor(TS1_dat$coolerno)
TS1_dat$tankno <- as.factor(TS1_dat$tankno)

#subset total data for leaf 2 -- includes PT data & Severity sampled from leaf 2 & totalSeverity summed for all leaves
leaf2_dat <- TS1_dat[TS1_dat$leafno==2,]

###SEVERITY

#make new variable for severity, to make non-zero 
leaf2_dat$lntotalSeverity=log(leaf2_dat$totalSeverity+.005)

#lntotalSeverity Plot: Control vs. Treatment, by time point 
ggplot(leaf2_dat,aes(x=timepoint, y=lntotalSeverity, fill=TC)) + geom_boxplot(outlier.size=-1) +
  xlab("Time After Exposure") + ylab("Disease Severity (% Lesion Coverage)") + ggtitle("Disease Severity")

#Total Severity Linear Model: two way ANOVA for time and treatment (ie fixed effects), random effects of cooler 
Severitytotal.model=lmer(totalSeverity~timepoint*TC+(1|tankno),  data=leaf2_dat)

#Check model assumptions 
###Normality of residuals 
hist(resid(Severitytotal.model))
###Homogenous variance 
plot(predict(Severitytotal.model),resid(Severitytotal.model))

#Severity Model ANOVA
summary(Severitytotal.model)
print(Severitytotal.model, correlation=TRUE)
Anova(Severitytotal.model, type="III")
emmip(Severitytotal.model, TC~timepoint, type ="response",CIs=T)
emmip(Severitytotal.model, TC~timepoint, type ="response",CIs=T, xlab="Time Point", ylab="Severity Modelled Response")
emmeans(Severitytotal.model, pairwise~TC|timepoint, type ="response")
emmeans(Severitytotal.model, pairwise~timepoint|TC, type ="response")


#PHENOLS

#Phenols Plot: Control vs. Treatment, by time point
ggplot(leaf2_dat,aes(x=timepoint, y=log(PT.), fill=TC)) + geom_boxplot(outlier.size=-1) +
  xlab("Time After Exposure") + ylab("Phenolic Acids (PT%)") + ggtitle("Phenolic Acid Concentration")

#Phenol Linear Model: two way ANOVA for time and treatment (ie fixed effects), random effects of tank
#With time_point
PT.model <- lmer(log(PT.) ~ timepoint*TC + (1|tankno) ,  data=leaf2_dat)

#Check model assumptions 
###Normality of residuals 
hist(resid(PT.model))
###Homogenous variance 
plot(predict(PT.model),resid(PT.model))
#Phenol Model ANOVA
summary(PT.model)
print(PT.model, correlation=TRUE)
Anova(PT.model, type="III")
emmip(PT.model, TC~time_after_exposure, type ="response",CIs=T, xlab="Time After Exposure (hrs)", ylab="Phenol Predicted Response", at=list(time_after_exposure=c(5,10,20,30,40,60,80,100,300,500)))
emmeans(PT.model, pairwise~TC|timepoint, type ="response")
emmeans(PT.model, pairwise~timepoint|TC, type ="response")

#AREA NEW GROWTH
#New Growth Plot: Control vs. Treatment, by time point
ggplot(TS1_dat,aes(x=timepoint, y=areaNG, fill=TC)) + geom_boxplot(outlier.size=-1) +
  xlab("Time After Exposure") + ylab("Area of New Growth (mm^2)") + ggtitle("Area of New Growth")

#New Growth Linear Model: two way ANOVA for time and treatment (ie fixed effects), random effects of cooler 
##NewGrowth.model <- areaNG ~ timepoint*TC + (1|tankno) ,  data=TS1_dat)
NewGrowth.model=lmer(areaNG~timepoint*TC+(1|tankno),  data=TS1_dat)
summary(NewGrowth.model)
print(NewGrowth.model, correlation=TRUE)
Anova(NewGrowth.model, type="III")
emmip(NewGrowth.model, TC~timepoint, type ="response",CIs=T, xlab="Time Point", ylab="New Growth Predicted Response")
pairs(emmeans(NewGrowth.model, ~TC|timepoint, type ="response"))
emmeans(NewGrowth.model, pairwise~TC|timepoint, type ="response")

#TOTAL AREA NEW GROWTH
#Total New Growth Plot: Control vs. Treatment, by time point
ggplot(leaf2_dat,aes(x=timepoint, y=totalareaNG, fill=TC)) + geom_boxplot(outlier.size=-1) +
  xlab("Time After Exposure") + ylab("Area of New Growth (mm^2)") + ggtitle("Area of New Growth")

#New Growth Linear Model: two way ANOVA for time and treatment (ie fixed effects), random effects of cooler 
##NewGrowth.model <- areaNG ~ timepoint*TC + (1|tankno) ,  data=TS1_dat)
TotalNewGrowth.model=lmer(totalareaNG~timepoint*TC+(1|tankno),  data=TS1_dat)
summary(TotalNewGrowth.model)
print(TotalNewGrowth.model, correlation=TRUE)
Anova(TotalNewGrowth.model, type="III")
emmip(TotalNewGrowth.model, TC~timepoint, type ="response",CIs=T, xlab="Time Point", ylab="New Growth Predicted Response")
pairs(emmeans(TotalNewGrowth.model, ~TC|timepoint, type ="response"))
emmeans(TotalNewGrowth.model, pairwise~TC|timepoint, type ="response")

#Leaf 2 Severity Linear Model: two way ANOVA for time and treatment (ie fixed effects), random effects of cooler 
leaf2_dat$lnSeverity=log(leaf2_dat$Severity+.005)
#lnSeverity Plot: Control vs. Treatment, by time point 
ggplot(leaf2_dat,aes(x=timepoint, y=lnSeverity, fill=TC)) + geom_boxplot(outlier.size=-1) +
  xlab("Time After Exposure") + ylab("Disease Severity (% Lesion Coverage)") + ggtitle("Disease Severity")
#Model
Severityleaf2.model=lmer(Severity~timepoint*TC+(1|tankno),  data=leaf2_dat)
#Check model assumptions 
###Normality of residuals 
hist(resid(Severityleaf2.model))
###Homogenous variance 
plot(predict(Severityleaf2.model),resid(Severityleaf2.model))
#Severity Model ANOVA
summary(Severityleaf2.model)
print(Severityleaf2.model, correlation=TRUE)
Anova(Severityleaf2.model, type="III")
emmip(Severityleaf2.model, TC~timepoint, type ="response",CIs=T)
emmeans(Severityleaf2.model, pairwise~TC|timepoint, type ="response")
emmeans(Severityleaf2.model, pairwise~timepoint|TC, type ="response")

#Interaction between Phenolics and Severity? 
PT.Severitymodel <- lmer(log(PT.) ~ timepoint*TC + totalSeverity*timepoint + (1|tankno) ,  data=leaf2_dat)
summary(PT.Severitymodel)
print(PT.Severitymodel, correlation=TRUE)
Anova(PT.Severitymodel, type="III")
emmip(PT.Severitymodel, TC~timepoint, type ="response",CIs=T, xlab="Time Point", ylab="New Growth Predicted Response")
pairs(emmeans(PT.Severitymodel, ~TC|timepoint, type ="response"))
emmeans(PT.Severitymodel, pairwise~TC|timepoint, type ="response")

#NO SIGNIFICANT ASSOCIATION BETWEEN SEVERITY AND PHENOLICS, if you add severity into the model, it does not explain phenolics