setwd("C:/Users/harperl/Documents/CRCP Files")
Recruits<-read.csv("Recruits.csv")
Live<-read.csv("LiveScans.csv")
MapCoords<-read.csv("Sites List Edit.csv") #coordinates nudged manually (as necessary and as little as possible) to eliminate point overlap
nums<-read.csv("Site Numbers.csv") #Site numbers in order of descending latitude, corresponding with Table 2

devtools::install_github("jslefche/piecewiseSEM@devel")

library(nlme)
library(rlang)
library(plyr)
library(tidyverse)
library(reshape2) 
library(viridis)

remotes::install_github("coolbutuseless/ggpattern")
library(ggpattern)

se<-function(x)sqrt(var(x)/length(x))

isSig <- function(p) {
  
  ifelse(p > 0.01 & p < 0.05, "*",
         ifelse(p > 0.001 & p <= 0.01, "**",
                ifelse(p <= 0.001, "***", "")))
  
}



#correct for errors in total corals
Recruits <- Recruits %>% mutate(TotalSto = AGAR+PORI+SIDE+FAVI+UNKS+Other)


#LongForm
Long <-melt(Recruits, id.vars =c("Year","Region","Habitat","Depth",
                                 "Site","Tile","Loc"), measure.vars = 
              c("AGAR", "SIDE","PORI","FAVI",
                "Other","UNKS","TotalSto","TotalOct"))

Long <- Long %>% rename("Taxon" = "variable", "Count" = "value")

Long$Taxon <- recode(Long$Taxon, "AGAR" = "Agariciidae",
                     "SIDE" = "Siderastreidae",
                     "PORI" = "Poritidae",
                     "FAVI" = "Faviidae",
                     "UNKS" = "Unknown Scleractinian", 
                     "TotalOct" = "Octocoral")

Long <- Long %>%
  mutate(Area = case_when(
    Loc == "Top" ~ 0.0225,
    Loc == "Bottom" ~ 0.0225,
    Loc == "Side" ~  0.012)) %>%
  mutate(Density = Count/Area) %>%
  mutate(Method = "Corallite")


Long <- Long %>%
  mutate(Loc = fct_relevel(Loc, 
                           "Top", "Side", "Bottom")) 

Reps <- Long %>% group_by(Year) %>% summarize(n = length(unique(Tile)))

Live <- Live %>% mutate(Method = "Live") %>%
  rename("Density" = "LiveArea",
         "Loc" = "Top.Bottom",
         "Tile" = "Tile.ID",
         "Site" = "SiteName",
         "Count" = "value",
         "Taxon" = "variable")

LiveRec <- rbind.fill(Long, Live)

LiveRec$Taxon <- recode(LiveRec$Taxon, "TotalSto" = "Stony",
                        "Octo" = "Octocoral")


mapdat <- MapCoords %>%
  mutate(Depth = Depth*0.3)

Taxa <- Long %>% subset(Taxon == "Agariciidae" | Taxon == "Poritidae" | Taxon == "Siderastreidae"|
                          Taxon == "Faviidae") %>% droplevels()
Taxa <- left_join(Taxa, mapdat, by = c("Site"="Site.Name")) %>%
  rename("Depth_m" = "Depth.y")


Octos <- Long %>% subset(Taxon == "Octocoral") %>% droplevels()
Octos <- left_join(Octos, mapdat, by = c("Site"="Site.Name")) %>%
  rename("Depth_m" = "Depth.y")

#proportion of corals on bottom of tile vs top or sides: Stony Corals
Df <- Taxa %>% dplyr::select("Year":"Habitat", "Site", "Depth_m", "Tile", "Taxon", "Loc", "Count") %>%
  mutate(row = row_number())

PivotLoc <- Df %>% pivot_wider(names_from = Loc, values_from = Count, values_fill = 0)

TaxTileDepth <- PivotLoc %>% group_by(Depth_m, Taxon, Year, Tile) %>%
  summarize(Top = sum(Top), Bottom = sum(Bottom), Side = sum(Side))

TaxSiteDepth <- PivotLoc %>% group_by(Depth_m, Taxon, Year, Site) %>%
  summarize(Top = sum(Top), Bottom = sum(Bottom), Side = sum(Side))

TaxLocDepth <- PivotLoc %>% group_by(Depth_m, Taxon) %>%
  summarize(Top = sum(Top), Bottom = sum(Bottom), Side = sum(Side))

#proportion of corals on bottom of tile vs top or sides: Octocorals
DfO <- Octos %>% select("Year":"Habitat", "Site", "Depth_m", "Tile", "Taxon", "Loc", "Count") %>%
  mutate(row = row_number())

PivotLocO <- DfO %>% pivot_wider(names_from = Loc, values_from = Count, values_fill = 0) %>%
  select(-c(Side))

OctoTileDepth <- PivotLocO %>% group_by(Depth_m, Taxon, Year, Tile) %>%
  summarize(Top = sum(Top), Bottom = sum(Bottom))

OctoSiteDepth <- PivotLocO %>% group_by(Depth_m, Taxon, Year, Site) %>%
  summarize(Top = sum(Top), Bottom = sum(Bottom))

OctoLocDepth <- PivotLocO %>% group_by(Depth_m, Taxon) %>%
  summarize(Top = sum(Top), Bottom = sum(Bottom))


library(splitstackshape)
library(factoextra)
library(gridExtra)
library(lemon)
library(multcompView)
library(rcompanion)
library(cowplot)
library(FSA)
library(scales)
library(gg.gap)
library(effects)
library(emmeans)
library(lme4)
library(lmtest)
library(piecewiseSEM)
library(MASS)
library(ggmap)
library(MuMIn)
library(car)
library(ggrepel)

###############################
#Table 1: Descriptive Stats from Our Study

stopercent <- Long %>% subset(Taxon == "TotalSto") %>%
  group_by(Loc) %>% summarize(TotalLoc = sum(Count)) %>%
  mutate(Total = sum(TotalLoc)) %>%
  mutate(Percent = (TotalLoc/Total)*100)

octpercent <- OctoTileDepth %>%
  group_by(Taxon) %>%
  summarize(TotalTop = sum(Top), TotalBot = sum(Bottom)) %>%
  mutate(Total = TotalTop + TotalBot) %>%
  mutate(Percent = (TotalTop/Total)*100)

###############################
#Figure 1
###############################
mapdat <- left_join(mapdat, nums, by = "Site.Name")

FLMap <- get_stamenmap(bbox = c(left = -82.1, bottom = 24.3, right = -79.9, top = 26.3), maptype = "terrain")

plot(FLMap)

tiff("SiteMapDepth.tif",width = 7, height = 7, units = "in", res = 800)
ggmap(FLMap)  +
  geom_point(data=mapdat, aes(x=Longitude, y = Latitude, fill = Depth),
             color = "white", size = 3.7, pch = 21, stroke = 1.1, alpha = 0.8) +
  geom_text_repel(data=mapdat, aes(x=Longitude, y=Latitude, label=No), 
                  size = 3.5, fontface = "bold", color="black", hjust=-0.4) +
  labs(y = "Latitude", x = "Longitude") +
  scale_fill_gradient("Depth (m)", low = "lightgoldenrod1", high = "black") +
  theme(legend.title = element_text(size=12), 
        legend.text = element_text(size=10),
        legend.key=element_rect(fill='skyblue4'),
        plot.margin = unit(c(0,0,0,0), "cm")) +
  guides(fill = guide_legend(override.aes = list(size=4, stroke=1.5))) 
dev.off()

#####################
#Figure 2: Images Assembled in Powerpoint
#####################

#####################
#Figure 3
#####################

LiveRec <- LiveRec %>%
  mutate(Method = fct_relevel(Method, 
                              "Live", "Corallite")) 

LiveRec <- LiveRec %>% mutate(Loc = fct_relevel(Loc, "Top", "Bottom"))


dotsum <- LiveRec %>% subset(Taxon == "Stony") %>% subset(Loc != "Side") %>% droplevels() %>%
  group_by(Year, Loc, Method, Site) %>%
  summarize(SiteDens = mean(Density))

dotsum2 <- LiveRec %>% subset(Taxon == "Stony") %>% subset(Loc != "Side") %>% droplevels() %>%
  group_by(Year, Loc, Method) %>%
  summarize(total = sum(Count))

dotsum3 <- LiveRec %>% subset(Taxon == "Stony") %>% subset(Loc != "Side") %>% droplevels() %>%
  subset(Year == "2018" & Loc == "Bottom") %>%
  group_by(Year, Loc, Method, Site) %>%
  summarize(n = n()) %>%
  arrange(Site, Loc)


Fig3 <- ggplot(dotsum, aes(x = Loc, y = SiteDens, fill = Method)) +
  geom_boxplot(color = "black") +
  facet_grid(~Year, switch = "x") +
  scale_y_continuous(expand = c(0,0), limits = c(0,1000)) +
  geom_vline(xintercept=c(2.6), linetype="dashed") +
  scale_x_discrete("") +
  scale_fill_manual("Scan Method",values=c("gray40","white")) +
  labs(y=expression(Recruits/m^2)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line.x = element_line(colour = "black"),
        axis.line.y = element_line(colour = "black"),
        axis.text = element_text(colour = "black", size = 12), axis.title.y = element_text(colour = "black", size = 13),
        axis.ticks.length.x = unit(0, 'cm'),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 13),
        strip.placement = "bottom", strip.background = element_rect(colour = "black", fill = "gray40"),
        strip.text = element_text(colour = "white", size = 13)) 

tiff("Fig3.tif",width = 9, height = 6, units = "in", res = 300)
Fig3
dev.off()

outliers3 <- dotsum %>% arrange(-(SiteDens))

#####################
##Table 2/Supporting analysis for Figure 3###
#####################
LiveRecStony <- LiveRec %>% subset(Taxon == "Stony") %>% subset(Loc != "Side") %>% droplevels()


norm.df <- LiveRecStony %>%
  group_by(Method) %>%
  summarise(statistic = shapiro.test(log(Density+1))$statistic, 
            p.value = shapiro.test(log(Density+1))$p.value) #not normal



onedf <- LiveRecStony %>% subset(Year == "2016")
twodf <- LiveRecStony %>% subset(Year == "2017")
threedf <- LiveRecStony %>% subset(Year == "2018")

totals1 <- onedf %>% group_by(Method, Loc) %>% summarize(sum = sum(Count))
totals2 <- twodf %>% group_by(Method, Loc) %>% summarize(sum = sum(Count))
totals3 <- threedf %>% group_by(Method, Loc) %>% summarize(sum = sum(Count))


#conduct separate linear model for each year
##2016

Fig1mod16 <- glm.nb(Count ~ Loc*Method + Site,
                    link = "log",
                    data = onedf)

# Examine assumptions
hist(resid(Fig1mod16))

# Get model results
Anova(Fig1mod16, type = "II")

summary(Fig1mod16)
rsquared(Fig1mod16)

pairwise16 <- as.data.frame(pairs(emmeans(Fig1mod16, ~ Loc*Method)))

pairwise16$sig <- isSig(pairwise16$p.value)

cldList(p.value ~ contrast,
        data = pairwise16,
        threshold = 0.05)


##2017
Fig1mod17 <- glm.nb(Count ~ Loc*Method + Site,
                    link = "log",
                    data = twodf)
# Examine assumptions
hist(resid(Fig1mod17))


# Get model results
Anova(Fig1mod17, type = "II")


rsquared(Fig1mod17)


pairwise17 <- as.data.frame(pairs(emmeans(Fig1mod17, ~ Loc*Method)))


pairwise17$sig <- isSig(pairwise17$p.value)

cldList(p.value ~ contrast,
        data = pairwise17,
        threshold = 0.05)

check <- twodf %>% group_by(Loc) %>% summarize(mean = mean(Count), se = se(Count))
##2018

Fig1mod18 <- glm.nb(Count ~ Loc*Method + Site,
                    link = "log",
                    data = threedf)
# Examine assumptions
hist(resid(Fig1mod18))


# Get model results
Anova(Fig1mod18, type = "II")


rsquared(Fig1mod18)


pairwise18 <- as.data.frame(pairs(emmeans(Fig1mod18, ~ Loc*Method)))

pairwise18$sig <- isSig(pairwise18$p.value)

cldList(p.value ~ contrast,
        data = pairwise18,
        threshold = 0.05)


##########################
##Table 3
##########################
##Model with Year as Fixed Factor
Table3Mod <- glm.nb(Count ~ Loc*Method + as.factor(Year) + Site,
                    link = "log", maxit = 100,
                    data = LiveRecStony)

# Examine assumptions
hist(resid(Table3Mod))


# Get model results
Anova(Table3Mod, type = "II")


rsquared(Table3Mod)


#descriptive stats: percent found in live of total found in corallite
sum1 <- LiveRecStony %>% 
  group_by(Method) %>% summarize(sum = sum(Count), dens = mean(Density))

(11059/2915)*100
(88.9/23.5)*100

###########################
##Figure 4
###########################

FamilyLoc <- Long %>% subset(Taxon == "Agariciidae"|Taxon == "Poritidae"|
                                 Taxon == "Siderastreidae"|Taxon == "Faviidae") %>% droplevels() %>%
  group_by(Loc, Taxon, Site) %>%
  summarize(SiteDens = mean(Density))

#Get percentage of four primary families
TotIDSto <- Long %>% subset(Taxon != "TotalSto" & Taxon != "Octocoral" & Taxon != "Unknown Scleractinian") %>%
  summarize(Sto = sum(Count))
MainFam <- Long %>% subset(Taxon == "Agariciidae"|Taxon == "Poritidae"|
                         Taxon == "Siderastreidae"|Taxon == "Faviidae") %>%
  summarize(MF = sum(Count))

MainFam/TotIDSto*100

FamilyLoc$Loc <- fct_relevel(FamilyLoc$Loc,"Top", "Side", "Bottom")

FamilyLoc <- FamilyLoc %>%
  mutate(Taxon = fct_relevel(Taxon, "Siderastreidae", "Poritidae", "Agariciidae", "Faviidae"))

##Make one with facet labels
Fig4 <- ggplot(FamilyLoc, aes(x = Loc, y = SiteDens, fill = Loc)) +
  geom_boxplot() +
  facet_grid(~Taxon, switch = "x") +
  scale_y_continuous(expand = c(0,0), limits = c(0,100)) +
  scale_x_discrete("", labels = c("", "", "")) +
  scale_fill_manual("Tile Surface", values = c( Top = "gray90", Side = "yellow", Bottom = "blue")) +
  labs(y=expression(Recruits/m^2)) +
  geom_vline(xintercept=c(3.6), linetype="dashed") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line.x = element_line(colour = "black"),
        axis.line.y = element_line(colour = "black"),
        axis.text = element_text(colour = "black", size = 12), axis.title.y = element_text(colour = "black", size = 13),
        axis.ticks.length.x = unit(0, 'cm'),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 13),
        strip.placement = "bottom", strip.background = element_rect(colour = "black", fill = "gray40"),
        strip.text = element_text(colour = "white", size = 13)) 


tiff("Fig4.tif",width = 9, height = 6, units = "in", res = 300)
plot(Fig4)
dev.off()


outliers4 <- FamilyLoc %>% arrange(-(SiteDens))

#######################
##Table 4
#######################

#testing normality and homogeneity of variances

norm.df <- Long %>%
  group_by(Taxon, Loc) %>%
  summarise(statistic = shapiro.test(Count)$statistic, 
            p.value = shapiro.test(Count)$p.value)
#none of the groups are normally distributed
homo.df <- Long %>%
  group_by(Taxon) %>%
  summarise(statistic = bartlett.test(Count~Loc)$statistic, 
            p.value = bartlett.test(Count~Loc)$p.value)
#none of the groups are homogeneous

#try again sqrt-transformed
LongAreaSqrt <- Long %>% mutate(Count = sqrt(Count))
norm.df <- LongAreaSqrt %>%
  group_by(Taxon, Loc) %>%
  summarise(statistic = shapiro.test(Count)$statistic, 
            p.value = shapiro.test(Count)$p.value)
#none of the groups are normally distributed
homo.df <- LongAreaSqrt %>%
  group_by(Taxon) %>%
  summarise(statistic = bartlett.test(Count~Loc)$statistic, 
            p.value = bartlett.test(Count~Loc)$p.value)


#Linear Model: top vs bottom by taxon
FamLong <- Long %>% subset(Taxon == "Agariciidae"|Taxon == "Poritidae"|
                             Taxon == "Siderastreidae"|Taxon == "Faviidae") %>% droplevels()

mapdat <- mapdat %>% rename("Depth_m" = "Depth")

####
#print mapdat file for data submission with depth corrected
mapdat <- mapdat %>% dplyr::select("No", "Site.Name":"Longitude")
write.csv(mapdat, file = "Sites List Edit_V2.csv")

FamLongDepth <- left_join(FamLong, mapdat, by = c("Site" = "Site.Name")) %>%
  dplyr::select("Year":"Method", "Depth_m")


Table4Mod <- glm.nb(Count ~ offset(Area) + Taxon*Loc*Depth_m + as.factor(Year),
                   link = "log", maxit = 500,
                   data = FamLongDepth)
# Examine assumptions
hist(resid(Table4Mod))

# Get model results
Anova(Table4Mod, type = "II")

rsquared(Table4Mod)

pairwise <- as.data.frame(pairs(emmeans(Table4Mod, ~ Taxon*Loc)))

pairwise$sig <- isSig(pairwise$p.value)

pairwise2 <- do.call(rbind, lapply(1:nrow(pairwise), function(i) {
  
  a <- strsplit(pairwise[i, "contrast"], " ")[[1]]
  
  if(sub(" .*", "", a[3]) == sub(" .*", "", a[7])) data.frame(Taxon = sub(" .*", "", a[3]) , pairwise[i, ]) else data.frame()
  
} ) )


#get letters for each taxon
ag <- pairwise2 %>% filter(Taxon == "Agariciidae")
cldList(p.value ~ contrast,
        data = ag,
        threshold = 0.05)

por <- pairwise2 %>% filter(Taxon == "Poritidae")
cldList(p.value ~ contrast,
        data = por,
        threshold = 0.05)

fav <- pairwise2 %>% filter(Taxon == "Faviidae")
cldList(p.value ~ contrast,
        data = fav,
        threshold = 0.05)

sid <- pairwise2 %>% filter(Taxon == "Siderastreidae")
cldList(p.value ~ contrast,
        data = sid,
        threshold = 0.05)


########################
#Figure 5
########################
TaxSiteDepth <- TaxSiteDepth %>% mutate(Total = Top + Side + Bottom) %>%
  filter(Total != 0) %>%
  mutate(Prop_Bot = Bottom/Total,Prop_Top = Top/Total, Prop_Side = Side/Total) 

TopProp_p <- ggplot(TaxSiteDepth, aes(x = Depth_m, y = Prop_Top, color = Taxon)) +
  geom_point(size = 1, position = "jitter", alpha = 0.3) +
  geom_smooth(method = "lm", se = F) +
  scale_y_continuous("Prop. on Top") +
  scale_x_continuous("") +
  scale_color_viridis("Taxon", discrete = TRUE) +
  theme(plot.title = element_text(size = 12,hjust = 0.5),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        axis.text.y = element_text(colour = "black", hjust = 1, size = 11),
        axis.text.x = element_blank(),
        axis.title = element_text(size = 12), legend.position = "none")
plot(TopProp_p)


BotProp_p <- ggplot(TaxSiteDepth, aes(x = Depth_m, y = Prop_Bot, color = Taxon)) +
  geom_point(size = 1, position = "jitter", alpha = 0.3) +
  geom_smooth(method = "lm", se = F) +
  scale_y_continuous("Prop. on Bottom") +
  scale_x_continuous("Depth (m)") +
  scale_color_viridis("Taxon", discrete = TRUE) +
  theme(plot.title = element_text(size = 12,hjust = 0.5),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        axis.text.y = element_text(colour = "black", hjust = 1, size = 11),
        axis.text.x = element_text(colour = "black", hjust = 1, size = 11),
        axis.title = element_text(size = 12), legend.position = "none")
plot(BotProp_p)


legend <- ggplot(TaxSiteDepth, aes(x = Depth_m, y = Prop_Side, color = Taxon)) +
  geom_point(size = 2, position = "jitter") +
  geom_smooth(method = "lm", se = F) +
  scale_y_continuous("Proportion Settled on Side") +
  scale_x_continuous("Depth (m)") +
  scale_color_viridis("Taxon", discrete = TRUE) +
  theme(plot.title = element_text(size = 12,hjust = 0.5),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        axis.text.y = element_text(colour = "black", hjust = 1, size = 11),
        axis.text.x = element_blank(),
        axis.title = element_text(size = 12))
tiff("TaxDepthLocLegend.tif",width = 7, height = 7, units = "in", res = 800)
plot(legend)
dev.off()

SideProp_p <- ggplot(TaxSiteDepth, aes(x = Depth_m, y = Prop_Side, color = Taxon)) +
  geom_point(size = 1, position = "jitter", alpha = 0.3) +
  geom_smooth(method = "lm", se = F) +
  scale_y_continuous("Prop. on Side") +
  scale_x_continuous("") +
  scale_color_viridis("Taxon", discrete = TRUE) +
  theme(plot.title = element_text(size = 12,hjust = 0.5),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        axis.text.y = element_text(colour = "black", hjust = 1, size = 11),
        axis.text.x = element_blank(),
        axis.title = element_text(size = 12), legend.position = "none")
plot(SideProp_p)

tiff("TaxDepthLoc.tif",width = 7, height = 7, units = "in", res = 800)
grid.arrange(TopProp_p, SideProp_p, BotProp_p, nrow = 3)
dev.off()


#check effect of depth on proportion of each taxon on each suface
#p-values and r coefficients for Fig5


#Top

models <- TaxSiteDepth %>% 
  split(.$Taxon) %>% 
  map(~lm(Prop_Top ~ Depth_m, data = .))

summary(models$Agariciidae)
summary(models$Poritidae)
summary(models$Siderastreidae)
summary(models$Faviidae)

by(TaxSiteDepth, TaxSiteDepth$Taxon,function(x){cor((x$Prop_Top),(x$Depth_m), method = "pearson")})



models <- TaxSiteDepth %>% 
  split(.$Taxon) %>% 
  map(~lm(Prop_Side ~ Depth_m, data = .))

summary(models$Agariciidae)
summary(models$Poritidae)
summary(models$Siderastreidae)
summary(models$Faviidae)

by(TaxSiteDepth, TaxSiteDepth$Taxon,function(x){cor((x$Prop_Side),(x$Depth_m), method = "pearson")})


models <- TaxSiteDepth %>% 
  split(.$Taxon) %>% 
  map(~lm(Prop_Bot ~ Depth_m, data = .))

summary(models$Agariciidae)
summary(models$Poritidae)
summary(models$Siderastreidae)
summary(models$Faviidae)

by(TaxSiteDepth, TaxSiteDepth$Taxon,function(x){cor((x$Prop_Bot),(x$Depth_m), method = "pearson")})
#########################
#Figure 6
#########################

OctoSum <- Octos %>% subset(Loc != "Side") %>% droplevels() %>%
  group_by(Year, Loc, Site) %>%
  summarize(SiteDens = mean(Density), Count = sum(Count))

OctoStat <- Octos %>% subset(Loc != "Side") %>% droplevels()

OctoSum <- OctoSum %>%
  mutate(Loc = fct_relevel(Loc, 
                           "Top","Bottom")) 


Fig6 <- ggplot(OctoSum, aes(x = Loc,y = SiteDens, fill = Loc)) +
  facet_grid(~Year, switch = "x") +
  geom_boxplot() +
  scale_y_continuous(expand = c(0,0), limits = c(0,110)) +
  scale_fill_manual("Location on Tile",values=c("gray90", "blue"),
                    labels=c("Top","Bottom")) +
  labs(y=expression(Recruits/m^2), x = "") +
  geom_vline(xintercept=c(2.6), linetype="dashed") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line.x = element_line(colour = "black"),
        axis.line.y = element_line(colour = "black"),
        axis.text = element_text(colour = "black", size = 12), axis.title.y = element_text(colour = "black", size = 13),
        axis.ticks.length.x = unit(0, 'cm'),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 13),
        strip.placement = "bottom", strip.background = element_rect(colour = "black", fill = "gray40"),
        strip.text = element_text(colour = "white", size = 13)) 

tiff("Fig6.tif",width = 9, height = 6, units = "in", res = 300)
Fig6
dev.off()

outliers6 <- OctoSum %>% arrange(-(SiteDens))


###################
##Table 5
###################

#Generalized Linear Model: top vs bottom for octos

Table5Mod <- glm(Count ~ offset(Area) + Loc*Depth_m + as.factor(Year),
                    maxit = 200,
                    family = "quasipoisson",
                    data = OctoStat)

# Examine assumptions
hist(resid(Table5Mod))

summary(Table5Mod)
# Get model results
Anova(Table5Mod, type = "II")

rsquared(Table5Mod)

pairwise <- as.data.frame(pairs(emmeans(Fig5mod, ~ as.factor(Year)*Loc)))

pairwise$sig <- isSig(pairwise$p.value)

pairwise2 <- do.call(rbind, lapply(1:nrow(pairwise), function(i) {
  
  a <- strsplit(pairwise[i, "contrast"], "- ")[[1]]
  
  if(sub(" .*", "", a[1]) == sub(" .*", "", a[2])) data.frame(Year = sub(" .*", "", a[1]) , pairwise[i, ]) else data.frame()
  
} ) )

#check for relationship between proportion and depth

OctoSiteDepth <- OctoSiteDepth %>% mutate(Total = Bottom + Top) %>%
                   filter(Total != 0) %>%
                   mutate(Prop_Bot = Bottom/Total,Prop_Top = Top/Total) 

#top
modelO <- lm(Prop_Top ~ Depth_m, data = OctoSiteDepth)

summary(modelO)

cor(OctoSiteDepth$Prop_Top, OctoSiteDepth$Depth_m)

#bottom
modelO <- lm(Prop_Bot ~ Depth_m, data = OctoSiteDepth)

summary(modelO)

cor(OctoSiteDepth$Prop_Bot, OctoSiteDepth$Depth_m)



#Descriptive stats for Octos
OctoTotals <- OctoStat %>% group_by(Year, Loc) %>% summarize(sum = sum(Count))

OctoSum <- OctoStat %>% group_by(Year, Site) %>% summarize(Area = sum(Area), Count = sum(Count)) %>% 
   mutate(SiteDens = Count/Area) %>% arrange(SiteDens)

MeanOct <- OctoSum %>% group_by(Year) %>% summarize(yrmean = mean(SiteDens), se = se(SiteDens))

OctoSumSite <- OctoStat %>% group_by(Site) %>% summarize(mean = mean(Density)) %>% arrange(mean)
range <- OctoSumSite %>% filter(mean > 0.6 & mean < 4.3)

OctoSum <- OctoStat %>% group_by(Year, Site) %>% summarize(mean = mean(Density)) %>% arrange(mean)