#Environment, dosage and isolate moderate pathogen virulence in eelgrass wasting disease
#P.D. Dawkins , M.E. Eisenlord, R.M. Yoshioka, E. Fiorenza, S. Fruchter, F. Giammona, M. Winningham, C.D. Harvell1



####Dosage Analysis---------------

####Load Required Packages
rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
library(ggplot2)
library(lme4) 
library(MASS)
library(tidyr)
library(betareg)
require(lsmeans)
library(MuMIn)
library(betareg)
####Load a standard error function
SE=function(x){
  sqrt(var(x)/length(x))
}
BetaTrans=function(dat,group_id,s=.5){
  
  zz_dat=data.frame(dat,group_id,zz=(rep(1,length(group_id))))
  
  N=aggregate(zz~group_id,data=zz_dat,FUN=sum)
  zz_dat$N=N[match(zz_dat$group_id, N$group_id),2]
  yn=(zz_dat$dat*(zz_dat$N-1)+s)/zz_dat$N
}

#Load Data---------------------------
Dose=data.frame(read.table('Dosage_Dataset_EF.csv',header=T,sep=','))
Dose$Disease=0
Dose$Disease[which(Dose$LA>0)]=1
#Data Summary----------
Means=aggregate(Disease~Temp+Dose,data=Dose,FUN='mean')
colnames(Means)=c('Temp','Dose','Prev')
Means$PrevSE=aggregate(Disease~Temp+Dose,data=Dose,FUN='SE')[,3]
Means$Sev=aggregate(LA/BA~Temp+Dose,data=Dose,FUN='mean')[,3]
Means$SevSE=aggregate(LA/BA~Temp+Dose,data=Dose,FUN='SE')[,3]


#Is prevalence related to temperature or dosage---------------------
Prev.Mod21=glm(Disease~Dose*Temp,data=Dose,family='binomial')
Prev.Mod22=glm(Disease~Dose+Temp,data=Dose,family='binomial')
Prev.Mod23=glm(Disease~Dose,data=Dose,family='binomial')
Prev.Mod24=glm(Disease~Temp,data=Dose,family='binomial')
Prev.Mod25=glm(Disease~1,data=Dose,family='binomial')
AIC(Prev.Mod21,Prev.Mod22,Prev.Mod23,Prev.Mod24,Prev.Mod25)

PrevalenceMods=model.sel(Prev.Mod21,Prev.Mod22,Prev.Mod23,Prev.Mod24,Prev.Mod25)
write.table(PrevalenceMods, "PrevalenceModelTable.csv", col.names = TRUE, row.names = FALSE, sep = ',')




#Is Severity related to dosage or temperature-----------------
Dose$BA=round(Dose$BA)
Dose$LA=round(Dose$LA)
Dose$trans=BetaTrans(Dose$LA/Dose$BA,paste(Dose$Temp,Dose$Dose))

Sev.Mod31=betareg(trans~I(10^Dose)*Temp,data=Dose)
Sev.Mod32=betareg(trans~I(10^Dose)+Temp,data=Dose)
Sev.Mod33=betareg(trans~I(10^Dose),data=Dose)
Sev.Mod34=betareg(trans~Temp,data=Dose)
Sev.Mod35=betareg(trans~1,data=Dose)
AIC(Sev.Mod31,Sev.Mod32,Sev.Mod33,Sev.Mod34,Sev.Mod35)
SeverityMods=model.sel(Sev.Mod31,Sev.Mod32,Sev.Mod33,Sev.Mod34,Sev.Mod35)
summary(model.avg(SeverityMods))
write.table(SeverityMods, "SeverityModelTable.csv", col.names = TRUE, row.names = FALSE, sep = ',')

#Plot the data-------
pd <- position_dodge(.9)
ggplot(aes(x=Dose, y=Sev),data=Means)+geom_col(aes(fill=as.factor(Temp)),position = 'dodge')+geom_errorbar(aes(ymin=Sev-SevSE,ymax=Sev+SevSE,color=as.factor(Temp)),position=pd,width=.25)+scale_color_manual(values=c('black','black'))
ggplot(aes(x=Dose, y=Prev),data=Means)+geom_col(aes(fill=as.factor(Temp)),position = 'dodge')+geom_errorbar(aes(ymin=Prev-PrevSE,ymax=Prev+PrevSE,color=as.factor(Temp)),position=pd,width=.25)+scale_color_manual(values=c('black','black'))










####Isolate Virulence Analysis----------------------
###Load the required packages
rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
library(lme4)
library(lsmeans)
library(ggplot2)
library(pscl)
library(betareg)
library(MuMIn)
#bring in disease data
df1=data.frame(read.table('2016V.2.csv', header=TRUE, sep=','))
df1=df1[which(df1$treat!='8.16 ConD D'),]
SE=function(x)
{
  sqrt(var(x,na.rm=T)/length(x))
}
BetaTrans=function(dat,group_id,s=.5){
  
  zz_dat=data.frame(dat,group_id,zz=(rep(1,length(group_id))))
  
  N=aggregate(zz~group_id,data=zz_dat,FUN=sum)
  zz_dat$N=N[match(zz_dat$group_id, N$group_id),2]
  yn=(zz_dat$dat*(zz_dat$N-1)+s)/zz_dat$N
}

#Data calculations----------
df1$blade.area=round(df1$ba,0)
df1$lesion.area=round(df1$la,0)
head(df1)
#Determine if there are differences in prevalence between treatments--------
mod1=glm(prev~treat,family=binomial,data=df1)
mod1null=glm(prev~1,family=binomial,data=df1)
AIC(mod1)
AIC(mod1null)
summary(mod1)
anova(mod1null,mod1,test='LRT')
lsmeans(mod1,spec=pairwise~treat,adj="bonferroni")

#See if there are differences in severity between treatments-----------
df1$propla=df1$lesion.area/df1$blade.area
df1$trans=BetaTrans(df1$propla,group_id = df1$treat)

mod2=betareg(trans~treat,data=df1)
mod2null=betareg(trans~1,data=df1)
AIC(mod2,mod2null)

#Plot the data as prevalence versus severity-------
clustering.prev=aggregate(prev~treat,data=df1,FUN='mean')
clustering.sev=aggregate(lesion.area/blade.area~treat,data=df1[which(df1$prev==1),],FUN='mean')
clustering=merge(clustering.prev,clustering.sev)
colnames(clustering)=c('Treat','Prev','Sev')
clustering$PrevSE=aggregate(prev~treat,data=df1,FUN='SE')[,2]
clustering$SevSE=aggregate(lesion.area/blade.area~treat,data=df1[which(df1$prev==1),],FUN='SE')[,2]

ggplot(aes(x=Prev,y=Sev),data=clustering)+geom_point(aes(color=Treat),size=3)+geom_errorbarh(aes(xmin=Prev-PrevSE,xmax=Prev+PrevSE))+geom_errorbar(aes(ymin=Sev-SevSE,ymax=Sev+SevSE))


####In-Vivo Light and Temperature Experiment-----
rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
# Make a function for standard error
#  ( SE = standard deviation (=square root of the variance) / sqrt(n) )----
SE=function(x)
{
  sqrt(var(x)/length(x))
}
BetaTrans=function(dat,group_id,s=.5){
  
  zz_dat=data.frame(dat,group_id,zz=(rep(1,length(group_id))))
  
  N=aggregate(zz~group_id,data=zz_dat,FUN=sum)
  zz_dat$N=N[match(zz_dat$group_id, N$group_id),2]
  yn=(zz_dat$dat*(zz_dat$N-1)+s)/zz_dat$N
}

data=read.csv("LTexp2016.csv")
data[is.na(data)]=0

# Data Calculations----
data$CA=as.numeric(as.character(data$CA))
data$LA=as.numeric(as.character(data$LA))
data$percLA=as.numeric(as.character(data$percLA))
data$ApercLA=as.numeric(as.character(data$ApercLA))
data["prev"]=1
data$prev=ifelse(data$ApercLA==0,0,1)
data2=data[which(data$Treatment=='T'),]
data2$Temperature=as.numeric(data2$Temperature)

#Determine if light or temperature is influential on EWD prevalence
Prevalence.T.x.L=glm(prev~Temperature*Light,data = data2,family='binomial')
Prevalence.T.add.L=glm(prev~Temperature+Light,data = data2,family='binomial')
Prevalence.T=glm(prev~Temperature,data = data2,family='binomial')
Prevalence.L=glm(prev~Light,data = data2,family='binomial')
Prevalence.null=glm(prev~1,data = data2,family='binomial')
AIC(Prevalence.T.x.L,Prevalence.T.add.L,Prevalence.T,Prevalence.L,Prevalence.null)

#Determine if light or temperature is influential on severity 
data2$ApercLAt=BetaTrans(data2$ApercLA/100,data2$treat)

Severity.T.x.L=betareg(ApercLAt~as.factor(Temperature)*as.factor(Light),data = data2)
Severity.T.add.L=betareg(ApercLAt~as.factor(Temperature)+as.factor(Light),data = data2)
Severity.T=betareg(ApercLAt~as.factor(Temperature),data = data2)
Severity.L=betareg(ApercLAt~as.factor(Light),data = data2)
Severity.null=betareg(ApercLAt~1,data = data2)
AIC(Severity.T.x.L,Severity.T.add.L,Severity.T,Severity.L,Severity.null)
summary(Severity.T.x.L)
sev.mods=model.sel(Severity.T.x.L,Severity.T.add.L,Severity.T,Severity.L,Severity.null)
Avg.Sev=model.avg(sev.mods)
summary(Avg.Sev)
sense=aggregate(LA/CA~Temperature+Light,data=data2,FUN=mean)
sense$SE=aggregate(LA/CA~Temperature+Light,data=data2,FUN=SE)[,3]

ggplot(aes(x=Light,y=`LA/CA`),data=sense)+geom_col(aes(fill=as.factor(Temperature)),position = position_dodge(width=.9))+geom_errorbar(aes(ymin=`LA/CA`-SE,ymax=`LA/CA`+SE,color=as.factor(Temperature)),position = position_dodge(width=.9),width=.5)





####In-Vitro Light-Temperature experiment-------
# Clear Workspace
#  (gets rid of data floating around from other runs of R)
rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
# Make a function for standard error
#  ( SE = standard deviation (=square root of the variance) / sqrt(n) )
SE=function(x)
{
  sqrt(var(x)/length(x))
}

#Load the data----
Growth=data.frame(read.table('Cell_growth.csv',header=T,sep=','))
Growth=Growth[which(Growth$Date<3.1),]
Growth$Light2=ifelse(Growth$Light=='D',1,0)

#Run the full model set to determine which factors influence growth rate------------
T.x.L.x.D=glmer.nb(Number~Date*I(scale(Temperature))*Light2+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
T.x.D.add.D.x.L=glmer.nb(Number~Date*I(scale(Temperature))+Date*Light2+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
T.x.D.add.L=glmer.nb(Number~Date*I(scale(Temperature))+Light2+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
T.add.D.x.L=glmer.nb(Number~I(scale(Temperature))+Date*Light2+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
T.x.L.add.D=glmer.nb(Number~I(scale(Temperature))*Light2+Date+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
T.add.D.add.L=glmer.nb(Number~Date+I(scale(Temperature))+Light2+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
T.add.L=glmer.nb(Number~Temperature+Light2+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
T.add.D=glmer.nb(Number~Date+I(scale(Temperature))+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
D.add.L=glmer.nb(Number~Date+Light2+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
Temp=glmer.nb(Number~Temperature+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
Date=glmer.nb(Number~Date+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
Light=glmer.nb(Number~Light2+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))
Nul=glmer.nb(Number~1+(1|Vial),data=Growth, control = glmerControl(optimizer ="bobyqa"))

AIC(T.x.L.x.D,T.x.D.add.D.x.L,T.x.D.add.L,T.add.D.x.L,T.x.L.add.D,T.add.D.add.L,T.add.L,T.add.D,D.add.L,Temp,Date,Light,Nul)

summary(T.x.D.add.L)

#Plot the data and model results----
ggplot(aes(x=Date,y=Number),data=Growth)+geom_point(aes(color=as.factor(Temperature),shape=Light))
Growth.pred1=data.frame(Date=c(1,2,3),Temp=-0.9953596,Lig=0,Temperature=11,Light='L')
Growth.pred2=data.frame(Date=c(1,2,3),Temp= 0.9953596,Lig=0,Temperature=18,Light='L')
Growth.pred3=data.frame(Date=c(1,2,3),Temp=-0.9953596,Lig=1,Temperature=11,Light='D')
Growth.pred4=data.frame(Date=c(1,2,3),Temp= 0.9953596,Lig=1,Temperature=18,Light='D')
Growth.pred=rbind(Growth.pred1,Growth.pred2,Growth.pred3,Growth.pred4)
Growth.pred$Number=exp(0.1295+.1619*Growth.pred$Temp+1.3228*Growth.pred$Date+.6249*Growth.pred$Lig+.2178*Growth.pred$Temp*Growth.pred$Date)
Growth.pred$LL=exp(0.1295+.1619*Growth.pred$Temp+1.3228*Growth.pred$Date+.6249*Growth.pred$Lig+.2178*Growth.pred$Temp*Growth.pred$Date-sqrt(.2486^2+.1018^2+.2265^2+.1632^2+.1015^2))
Growth.pred$UL=exp(0.1295+.1619*Growth.pred$Temp+1.3228*Growth.pred$Date+.6249*Growth.pred$Lig+.2178*Growth.pred$Temp*Growth.pred$Date+sqrt(.2486^2+.1018^2+.2265^2+.1632^2+.1015^2))


ggplot(aes(x=Date,y=Number),data=Growth)+geom_jitter(height=0,width=.1,aes(color=as.factor(Temperature),shape=Light))+geom_line(aes(color=as.factor(Temperature),linetype=Light),data=Growth.pred)+geom_ribbon(aes(ymin=LL,ymax=UL,color=as.factor(Temperature),fill=Light),data=Growth.pred)+xlab('Day')+ylab('Number of Cells')