library(mrds)
library(lubridate)
library(tidyverse)
library(readxl)

#this contains functions to format the data and run the analyses
## This section runs the density analyses.
source('HelperFuncs.R') ##assumes all files are in the currect working directory


#run all estimates on Lake Burgan
LB.distance.est   <- distance.dens.est(curr.lake="Lake Burgan", curr.lake2="Burgan") #distance estimate
LB.quadrat.est    <- quadrat.dens.est(curr.lake="Lake Burgan", curr.lake2="Burgan")#quadrat estimate
LB.removal.est     <- removal.dens.est(curr.lake="Lake Burgan", curr.lake2="Burgan")#removal estimate


#run all estimates on Little Birch Lake
LBL.distance.est   <- distance.dens.est(curr.lake="Little Birch Lake", curr.lake2="Little Birch Lake")
LBL.quadrat.est    <- quadrat.dens.est(curr.lake="Little Birch Lake", curr.lake2="Little Birch Lake")
LBL.removal.est    <- removal.dens.est(curr.lake="Little Birch Lake", curr.lake2="Little Birch Lake")


#run all estimates on Lake Florida
LF.distance.est   <- distance.dens.est(curr.lake="Lake Florida", curr.lake2="Florida")
LF.quadrat.est    <- quadrat.dens.est(curr.lake="Lake Florida", curr.lake2="Florida")
LF.removal.est     <- removal.dens.est(curr.lake="Lake Florida", curr.lake2="Florida")


#combine relevant output into a dataframe
dens.df <- data.frame(Dhat= c(LF.distance.est$Dhat, LB.distance.est$Dhat, LBL.distance.est$Dhat, LF.removal.est$Dhat, LB.removal.est$Dhat, LBL.removal.est$Dhat, LF.quadrat.est$Dhat, LB.quadrat.est$Dhat, LBL.quadrat.est$Dhat), 
           Dhat.se= c(LF.distance.est$Dhat.se, LB.distance.est$Dhat.se, LBL.distance.est$Dhat.se, LF.removal.est$Dhat.se, LB.removal.est$Dhat.se, LBL.removal.est$Dhat.se, LF.quadrat.est$Dhat.se, LB.quadrat.est$Dhat.se, LBL.quadrat.est$Dhat.se), 
           Lake=rep(c("Lake Florida", "Lake Burgan", "Little Birch Lake"), 3), Method=c(rep("Distance", 3), rep("Removal", 3), rep("Quadrat", 3)) )


dens.df$Lake   <- factor(dens.df$Lake, levels=c("Lake Florida", "Lake Burgan", "Little Birch Lake"))
dens.df$Lake   <- relevel(dens.df$Lake, "Lake Florida")
dens.df$Method <- factor(dens.df$Method, levels=c("Quadrat", "Removal", "Distance"))


### plot the estimated densities
ggplot(dens.df, aes(x=Method, y=Dhat, colour=Method)) + 
  geom_point(size=4) + 
  facet_wrap(~Lake, ncol=3, scales="free") + 
  geom_errorbar(aes(ymin=Dhat - 2*Dhat.se, ymax=Dhat + 2*Dhat.se, width=0)) + 
  labs(x="Survey design", y=expression(paste("Estimated density (", hat(D),")" ))) + 
  #guides(fill=FALSE) +
  theme(legend.position="none")


### plot the estimated detection probabilities
phat.vec  <- c(LF.distance.est$phat, LF.removal.est$phat, LB.distance.est$phat, LB.removal.est$phat, LBL.distance.est$phat, LBL.removal.est$phat)

phat.se.vec <- c(LF.distance.est$phat.se, LF.removal.est$phat.se, LB.distance.est$phat.se, LB.removal.est$phat.se, LBL.distance.est$phat.se, LBL.removal.est$phat.se)

phat.df <- data.frame(phat =phat.vec, se=phat.se.vec, Lake=c("Lake Florida", "Lake Florida", "Lake Burgan", "Lake Burgan", "Little Birch Lake", "Little Birch Lake"), Design=rep(c("Distance", "Removal"), 3))

phat.df$Lake <- factor(phat.df$Lake, levels = c("Lake Florida", "Lake Burgan", "Little Birch Lake"))
ggplot(phat.df, aes(x=Lake, y=phat,  color=Design)) + 
  geom_point(size=4, position=position_dodge(width = .5)) +
  geom_errorbar(aes(ymin=phat-2*se, ymax=phat+2*se), width=0,  position=position_dodge(width = .5)) +
  labs(x="", y="Estimated probability of detection")


## This section conducts analysis of survey time by lake and method
time.df <- rbind(LB.distance.est$df, LB.removal.est$df, LB.quadrat.est$df, LBL.distance.est$df ,LBL.removal.est$df, LBL.quadrat.est$df, LF.distance.est$df, LF.removal.est$df, LF.quadrat.est$df)

time.df$Lake <- as.factor(time.df$Lake)
time.df$Lake <- relevel(time.df$Lake, ref="Lake Florida")

levels(time.df$Type)[levels(time.df$Type)=="removal"] <- "Removal" 

summary(aov(lm(log(Time/Length) ~ Type*Lake , data=time.df)))

summary(lm(log(Time) ~ Type*Lake - Type - Lake -1, data=time.df))


ggplot(time.df, aes(x=Lake, y=(Time/60^2), fill=Type)) +  
  geom_boxplot() + 
  labs(x="Lake", y = "Transect survey time (hours)") + 
  geom_point(size=1, alpha=0.5, aes(fill=Type), position = position_jitterdodge()) +
  theme(plot.title = element_text(hjust = 0.5)) + 
  scale_y_continuous(trans='log10')