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)

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

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

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

glm.douse2<-glm(Prevalence~Species,family="binomial",data=douse.dict)
anova(glm.douse2,test="Chi")
summary(glm.douse2) #aic=236.65#

glm.douse3<-glm(Prevalence~Treatment*Species,family="binomial",data=douse.dict)
anova(glm.douse3,test="Chi")
summary(glm.douse3) #aic=244.71#

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

glm.douse5<-glm(Prevalence~Species+Treatment+Time,family="binomial",data=douse.dict)
anova(glm.douse5,test="Chi")
summary(glm.douse5) #aic=111.34#

glm.douse6<-glm(Prevalence~Species*Treatment*Time,family="binomial",data=douse.dict)
anova(glm.douse6,test="Chi")
summary(glm.douse6) #aic=124.79#


##the glm model with the lowest AIC is the best fit.#
#interaction with all 3 variables is most accurate- douse5#
summary(glm.douse5)
TukeyHSD(glm.douse5)

#graph probability of predation based on best fit model#
visreg(glm.douse,"Species",by="Treatment",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")