##############################################
#  title: "PeerJ Supplemental Material"  #
#  author: "Cori Speights & Michael McCoy"   #
##############################################


#Data Input.


stone = read.csv(("StoneCrabSponges.csv"), header =T)
stone


#Isolating only data taken at the 24 hour mark. After 48 hours crabs at highest treatment (22 oysters) had eaten most/all which would make us unable to use the functional response curves.


new <- droplevels(subset(stone,hr==24))


#Necessary Libraries.


library(ggplot2)
library(grid)
library("bbmle")
require("bbmle")
library("emdbook")
require("emdbook")


#Rogers random predation model.


rogers.pred <- function(N0,a,h,T) {
  N0 - lambertW(a*h*N0*exp(-a*(T-h*N0)))/(a*h)
}


#Necessary variable assignments for the Rogers random predation model.


N0 <- new$oysnum
a <- 1
h <- 0.2
T0 <- 1 

eaten <- new$eaten
sdat <- data.frame(N0,eaten,sponge=new$sponge) 


#Two stone crabs did not eat during experiment or intermediate days so they were removed from the trials and data analysis.


nonew<-new[-c(4,14),]
nosdat<-sdat[-c(4,14),]


#Models fit using methods of maximum likelihood using the bbmle package.


(m8 <- mle2(eaten~dbinom(size=N0,prob=rogers.pred(N0,a,h,T)/N0), 
            #method="SANN",
            parameters=list(a~sponge, h~sponge),
            start=list(a=c(1.5,1.5),h=c(0.5,0.5)),
            data=c(as.list(nosdat),T=T0)))

(m9 <- mle2(eaten~dbinom(size=N0,prob=rogers.pred(N0,a,h,T)/N0), 
            #method="SANN",
            parameters=list(h~sponge),
            start=list(a=c(1.5),h=c(0.5,0.5)),
            data=c(as.list(nosdat),T=T0)))

(m10 <- mle2(eaten~dbinom(size=N0,prob=rogers.pred(N0,a,h,T)/N0), 
             #method="SANN",
             parameters=list(a~sponge),
             start=list(a=c(1.5,1.5),h=c(0.5)),
             data=c(as.list(nosdat),T=T0)))

(m11 <- mle2(eaten~dbinom(size=N0,prob=rogers.pred(N0,a,h,T)/N0), 
             #method="SANN",
             start=list(a=c(1.5),h=c(0.5)),
             data=c(as.list(nosdat),T=T0)))


#Model comparisons using size-corrected Akaike Information Criterion (AICc).


AICctab(m8,m9,m10,m11,delta=T,weights=T,nobs=38)


#Calculating the confindence intervals for the attack rate and handling time model parameters.


confint(m8,method="quad")
confint(m9,method="quad")
confint(m10,method="quad")
confint(m11,method="quad")


#Calculate mean and standard error for the no sponge treatments.


newns <- droplevels(subset(nonew,sponge=="NS"))
N0 <- newns$oysnum
eaten <- newns$eaten

newns.mean=with(newns,aggregate(eaten,list(c(N0)),mean)) #vector of means
sterr=function(x){sd(x)/sqrt(length(x))} #st. error
newns.err=with(newns,aggregate(eaten,list(c(N0)),sterr)) #vecotr of st. error
xx=data.frame(N0=newns.mean[,1],eaten=newns.mean[,2],err=newns.err[,2]) #data frame
limits=aes(ymax = eaten + err, ymin=eaten - err)



#Calculate mean and standard error for the sponge treatments.

news <- droplevels(subset(nonew,sponge=="S"))
N0 <- news$oysnum
eaten <- news$eaten

news.mean=with(news,aggregate(eaten,list(c(N0)),mean)) #vector of means
sterr=function(x){sd(x)/sqrt(length(x))} #st. error
news.err=with(news,aggregate(eaten,list(c(N0)),sterr)) #vecotr of st. error
xy=data.frame(N0=news.mean[,1],eaten=news.mean[,2],err=news.err[,2]) #data frame
limits=aes(ymax = eaten + err, ymin=eaten - err) ###CHECK THIS OUTS####


#Figure code using parameter estimates from the m1 model.

g2=ggplot(xy,aes(N0,eaten))+
  geom_point(col="black",size=9) +  ylim(c(0,20))+
  labs(x = "Initial Prey Density",y="Prey Consumed") + geom_errorbar(limits,  width=0.25)
g3=g2+ geom_errorbar(limits,  width=0.25)+stat_function(fun = function(x,a,h,T) { x - lambertW(a*h*x*exp(-a*(T-h*x)))/(a*h)}, 
                                                        args=list(a=1.970 ,h=.039,T=1),colour = "black",size=2) + 
  theme(panel.grid.minor=element_blank(),panel.grid.major=element_blank())+
  theme(axis.title.x = element_text(face="bold", size=25),
        axis.text.x  = element_text( size=20),
        axis.title.y = element_text(face="bold", angle=90,size=25),
        axis.text.y  = element_text( size=20))
g4=g3+geom_point(data=xx,aes(y=eaten, x=N0),col="gray",size=9)+geom_errorbar(data=xx,limits,  width=0.25)+ stat_function(data=xx,fun = function(x,a,h,T) { x - lambertW(a*h*x*exp(-a*(T-h*x)))/(a*h)}, 
                                                                                                                         args=list(a=4.079 ,h=0.112,T=1),colour = "gray",size=2) 


#Plot figure.


g4