setwd("C:/Users/caili/Documents/UMiami/Internships")
#install and run these libraries in order to complete analysis as written#
library(ggplot2)
#The MASS library is useful for general and generalized linear models
library(MASS)
library(gridExtra)
library(abd)
library(mgcv)
#The mgcv library fits GAM models
library(abd)
library(dplyr)
library(MASS)
library(rgl)
library(survival)
library(knitr)
library(Matrix)
library(mgcv)
library(visreg)
library(lattice)
library(ggpubr)
library(FSA)

##DICTYOTA FEEDING TRIAL##

#read in data set-ours are species specific- and make column titles simple#
#OFAV#
ofav.dict<-read.csv("OFAV.dict.csv")
ofav.dict<-rename(ofav.dict,Treatment=ï..Treatment)

#run generalized linear models that utilize different combinations of the independent variables- Treatment, Time- to determine best fit with dependent variable- Predation
glm.ofav.dict<-glm(Prevalence~Treatment,family="binomial",data=ofav.dict)
anova(glm.ofav.dict,test="Chi")
summary(glm.ofav.dict) #aic=324.08#

glm.ofav.dict1<-glm(Prevalence~Treatment+Time,family="binomial", data=ofav.dict)
anova(glm.ofav.dict1,test="Chi")
summary(glm.ofav.dict1) #aic=318.66#

glm.ofav.dict2<-glm(Prevalence~Time+Treatment,family="binomial", data=ofav.dict)
anova(glm.ofav.dict2,test="Chi")
summary(glm.ofav.dict2) #aic=318.66#

glm.ofav.dict3<-glm(Prevalence~Time,family="binomial",data=ofav.dict)
anova(glm.ofav.dict3,test="Chi")
summary(glm.ofav.dict3) #aic=338.26#

glm.ofav.dict4<-glm(Prevalence~Treatment*Time,family="binomial", data=ofav.dict)
anova(glm.ofav.dict4,test="Chi")
summary(glm.ofav.dict4) #aic=308.55#

##the glm model with the lowest AIC is the best fit.#
#interaction models appear to be best fits#
summary(glm.ofav.dict1)
TukeyHSD(glm.ofav.dict1)

summary(glm.ofav.dict4)
TukeyHSD(glm.ofav.dict4)

#graph probability of predation based on best fit model#
visreg(glm.ofav.dict1,"Treatment",by="Time",gg=TRUE, scale="response",partial=FALSE,rug=FALSE,line=list(col="black"),fill="Treatment")+
  theme(axis.line = element_line(colour = "black"))+
  theme_bw(base_size=15)+
  theme(axis.text.x = element_text(angle = 0, vjust = 0.5, hjust=0.5))+
  ylab("Probability of Predation")+
  ylim(0.0,1.00)

##REPEAT CODE ABOVE FOR EACH ADDITIONAL SPECIES##


#MCAV PREVALENCE#
mcav.dict<-read.csv("mcav.csv")
mcav.dict<-rename(mcav.dict,Treatment=ï..Treatment)

glm.mcav.dict<-glm(Prevalence~Treatment,family="binomial",data=mcav.dict)
anova(glm.mcav.dict,test="Chi")
summary(glm.mcav.dict) #aic=87.594#

glm.mcav.dict1<-glm(Prevalence~Treatment+Time,family="binomial", data=mcav.dict)
anova(glm.mcav.dict1,test="Chi")
summary(glm.mcav.dict1) #aic=89.254#

glm.mcav.dict2<-glm(Prevalence~Time+Treatment,family="binomial", data=mcav.dict)
anova(glm.mcav.dict2,test="Chi")
summary(glm.mcav.dict2) #aic=89.254#

glm.mcav.dict3<-glm(Prevalence~Time,family="binomial",data=mcav.dict)
anova(glm.mcav.dict3,test="Chi")
summary(glm.mcav.dict3) #aic=91.583#

glm.mcav.dict4<-glm(Prevalence~Treatment*Time,family="binomial", data=mcav.dict)
anova(glm.mcav.dict4,test="Chi")
summary(glm.mcav.dict4) #aic=97.254#

#best fit: treatment as only factor#
summary(glm.mcav.dict)
TukeyHSD(glm.mcav.dict)

#graph probability of predation#
visreg(glm.mcav.dict1,"Treatment",by="Time",gg=TRUE, scale="response",partial=FALSE,rug=FALSE,line=list(col="black"),fill="Treatment")+
  theme(axis.line = element_line(colour = "black"))+
  theme_bw(base_size=15)+
  theme(axis.text.x = element_text(angle = 0, vjust = 0.5, hjust=0.5))+
  ylab("Probability of Predation")

#PSTR PREVALENCE#
pstr.dict<-read.csv("pstr.csv")
pstr.dict<-rename(pstr.dict,Treatment=ï..Treatment)

glm.pstr.dict<-glm(Prevalence~Treatment,family="binomial",data=pstr.dict)
anova(glm.pstr.dict,test="Chi")
summary(glm.pstr.dict) #aic=263.36#

glm.pstr.dict1<-glm(Prevalence~Treatment+Time,family="binomial", data=pstr.dict)
anova(glm.pstr.dict1,test="Chi")
summary(glm.pstr.dict1) #aic=55.646#

glm.pstr.dict2<-glm(Prevalence~Time+Treatment,family="binomial", data=pstr.dict)
anova(glm.pstr.dict2,test="Chi")
summary(glm.pstr.dict2) #aic=55.646#

glm.pstr.dict3<-glm(Prevalence~Time,family="binomial",data=pstr.dict)
anova(glm.pstr.dict3,test="Chi")
summary(glm.pstr.dict3) #aic=53.516#

glm.pstr.dict4<-glm(Prevalence~Treatment*Time,family="binomial", data=pstr.dict)
anova(glm.pstr.dict4,test="Chi")
summary(glm.pstr.dict4) #aic=63.646#

#best fits?#
summary(glm.pstr.dict2)
TukeyHSD(glm.pstr.dict2)

summary(glm.pstr.dict3)
TukeyHSD(glm.pstr.dict3)

#graph probability of predation#
visreg(glm.pstr.dict2,"Treatment",by="Time",gg=TRUE, scale="response",partial=FALSE,rug=FALSE,line=list(col="black"),fill="Treatment")+
  theme(axis.line = element_line(colour = "black"))+
  theme_bw(base_size=15)+
  theme(axis.text.x = element_text(angle = 0, vjust = 0.5, hjust=0.5))+
  ylab("Probability of Predation")