### ANALYSES STRUCTURAL COMPLEXITY DATA
### July 2024

#-------------------------------------------------------------------
## Packages and setup
library(dplyr)
library(ggplot2)
library(tidyverse)
library(emmeans)
library(multcomp)
library(multcompView)
library(glmmTMB)
library(car)
library(vegan)

library(nlme)
library(MASS)
library(lmerTest)
library(pscl)
library(ggpubr)
library(patchwork)

#-------------------------------------------------------------------
## Data import and cleaning

## Load data
setwd("...~folder containing data")
data <- read.csv("Abundances.csv", sep=",", header=T, check.names=T)
mobile <- read.csv("Cover.csv", sep=",", header=T, check.names=T)

## Correct abundances and cover for missing pieces and subsamples taken
## Filter for only site Vlakte van Kerken
abundances_corrected <- abundances %>%
  separate(Site , c("Location", "Zone", "Platform"), sep= c(1, 2, 3)) %>%
  mutate(Barnacles=Barnacles/Subsample.barnacles)%>%
  mutate(across(c(Barnacles:Ensis.sp.), ~.*(1/(1-Missing)))) %>%
  mutate_if(is.character, as.factor)%>%
  mutate_at(c('Level', 'Vertical', 'Horizontal'), as.factor)%>%
  filter(Location=="K")

cover_corrected <- cover %>%
  separate(Site , c("Location", "Zone", "Platform"), sep= c(1, 2, 3)) %>%
  mutate(across(c(Missing:Small.anemones), ~./16)) %>%
  mutate(across(c(Algae:Small.anemones), ~.*(1/(1-Missing)))) %>%
  mutate_if(is.character, as.factor)%>%
  mutate_at(c('Level', 'Vertical', 'Horizontal'), as.factor)%>%
  filter(Location=="K")

#-------------------------------------------------------------------
## Group observations per structure

## Count the number of sides included in sample
no_sides <- abundances_corrected %>%
  filter(Side != "Residue")%>%
  group_by(Location,Zone, Level, Platform, Vertical, Horizontal)%>%
  count()
  
## Remove covering species in abundance data
## Add 0.25*residue for LV2 & LV3 and 1*residue for LV1
## Group per structure
## Round the numbers for poisson
pyramid_abundance <- merge(abundances_corrected, no_sides)%>%
  dplyr::select(-c("Seaweeds", "Hydroid.polyps"))%>%
  rowwise()%>%
  mutate(across(c(Barnacles),~ifelse(Side=="Residue", 0, .)))%>%
  mutate(across(c(Mussels:Ensis.sp.),~case_when(Level != "LV1" & Side=="Residue" & n==2 ~ .*0.5,
                                                 Level != "LV1" & Side=="Residue" & n==1 ~.*0.25, TRUE~.)))%>%
  group_by(Zone, Level, Platform)%>%
  summarise(across(Barnacles:Ensis.sp., sum))%>%
  mutate(across(c(Barnacles:Ensis.sp.),ceiling))%>%
  ungroup()

## Find abundances of each taxon
# mostabundant <- data.frame(colsums=colSums(pyramid_abundance[,c(5:23)]))

## Add number of sides assessed
## Group per structure & divide by 10 to get fraction coverage
pyramid_cover <- cover_corrected %>%
  filter(Side!="Residue") %>%
  group_by(Zone, Level, Platform, Vertical, Horizontal) %>%
  summarise(number=n(), across(c(Algae:Small.anemones), ~sum(.)/number))%>%
  group_by(Zone, Level, Platform) %>%
  summarise(across(c(Algae:Small.anemones), ~sum(.)/10))%>% ## transform to %? Then replace /10 by *10
  ungroup()

## Find coverage of each taxon
# mostcovered <- data.frame(colsums=colSums(pyramid_cover[,c(5:9)]))

#-------------------------------------------------------------------
## Analyse abundances of different taxa per treatment

## Barnacles
nb.barnaclemodel <- glmmTMB(data=pyramid_abundance, Barnacles~Level*Zone+(1|Platform), family="nbinom1") # fit model
# summary(nb.barnaclemodel) # model summary
car::Anova(nb.barnaclemodel) # find sigificance of predictors
barnaclemeans <- emmeans(nb.barnaclemodel, pairwise~Level*Zone) # post-hoc comparison
summary(contrast(barnaclemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr") # fdr corrected planned-paired comparison

pb <- data.frame(summary(contrast(barnaclemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr")) # save post-hoc groups for plotting
pb$labels <- symnum(pb$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s.")) # translate into significance-symbol

cldb <- data.frame(cld(object=barnaclemeans, adjust = "fdr", alpha=0.05, Letters = LETTERS))%>%
  rename(group = .group) # find statistical groups of treatments

## Plot the panel for taxon
A<-  ggplot(data=pyramid_abundance, aes(x=Level, y=Barnacles, fill=Zone))+
  geom_boxplot(outlier.shape = NA, position='dodge')+
  geom_point(position=position_jitterdodge())+
  labs(x= "", y="Barnacle Abundance")+ 
  scale_fill_manual(values=c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  ylim(-100, 3000)+
  theme_classic()+
  theme(text=element_text(size=15), legend.position = "none")+
  geom_signif(y_position = c(1400, 1800, 2800), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = pb$labels[c(7,8,9)])+
  geom_text(data=cldb, aes(y = c(-90, -90, -90, 250, 550,  1250), label = group), position = position_dodge(width = .75))

#--------------------------
## Mussels
nb.musselmodel <- glmmTMB(data=pyramid_abundance, Mussels~Level*Zone+(1|Platform), family="nbinom2") #fit is worse
car::Anova(nb.musselmodel)
musselmeans <- emmeans(nb.musselmodel, pairwise~Level*Zone)
summary(contrast(musselmeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr")

pm <- data.frame(summary(contrast(musselmeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr"))
pm$labels <- symnum(pm$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

cldm <- data.frame(cld(object=musselmeans, adjust = "fdr", alpha=0.04, Letters = LETTERS))%>%
  rename(group = .group)

B <- ggplot(data=pyramid_abundance, aes(x=Level, y=round(Mussels), fill=Zone))+
  geom_boxplot(outlier.shape = NA)+
  labs(x= "", y="Mussel Abundance")+ 
  geom_point(position=position_jitter())+
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  theme_classic()+
  theme(text=element_text(size=15), legend.position = "none")+
  geom_signif(y_position = c(8, 12, 54), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = pm$labels[c(7,8,9)])+
  geom_text(data=cldm, aes(y = c(-2, -2, -2, -2, 1,  2), label = group), position = position_dodge(width = .75))

#--------------------------
## Tunicates
zi.tunicatemodel2 <- glmmTMB(data=pyramid_abundance, Sea.squirt~Level*Zone+(1|Platform), family="nbinom2", ziformula=~1)
car::Anova(zi.tunicatemodel2)
tunicatemeans <- emmeans(zi.tunicatemodel2, pairwise~Level*Zone)
summary(contrast(tunicatemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr")

pt <- data.frame(summary(contrast(tunicatemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr"))
pt$labels <- symnum(pt$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

cldt <- data.frame(cld(object=tunicatemeans, adjust = "fdr", alpha=0.05, Letters = LETTERS))%>%
  rename(group = .group)

C <-  ggplot(data=pyramid_abundance, aes(x=Level, y=Sea.squirt, fill=Zone))+
  geom_boxplot(outlier.shape = NA)+
  geom_point(position=position_jitterdodge())+
  ylim(-2,33)+
  labs(x= "", y="Tunicate Abundance")+ 
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  theme_classic()+
  theme(text=element_text(size=15), legend.position = "none")+
  geom_signif(y_position = c(14, 25, 32), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = pt$labels[c(7,8,9)])+
  geom_text(data=cldt, aes(y = c(-1.5, -1.5, -1.5, -1.5, 4,  4), label = group), position = position_dodge(width = .75))

#--------------------------
## Anemones
nb.anemonemodel <- glmmTMB(data=pyramid_abundance, Anemones~Level*Zone+(1|Platform), family="nbinom2")
car::Anova(nb.anemonemodel)
anemonemeans <- emmeans(nb.anemonemodel, pairwise~Level*Zone)
summary(contrast(anemonemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr") 

pa <- data.frame(summary(contrast(anemonemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr"))
pa$labels <- symnum(pa$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

clda <- data.frame(cld(object=anemonemeans, adjust = "fdr", alpha=0.05, Letters = LETTERS))%>%
  rename(group = .group)

D<- ggplot(data=pyramid_abundance, aes(x=Level, y=Anemones, fill=Zone))+
  geom_boxplot(outlier.shape = NA)+
  labs(x= "", y="Anemone Abundance")+ 
  geom_point(position=position_jitterdodge())+
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  theme_classic()+
  theme(text=element_text(size=15),legend.position = "none")+
  geom_signif(y_position = c(15, 10, 8), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = pa$labels[c(7,8,9)])+
  geom_text(data=clda, aes(y = c(-1, -1, -1, -1, -1, -1), label = group), position = position_dodge(width = .75))

#--------------------------
## Effect sizes

## Calculate average abundances for each treatment
abundanceszl<- pyramid_abundance %>%
  group_by(Zone, Level) %>%
  summarise(Barnacles=mean(Barnacles, na.rm=T), Mussels=mean(Mussels, na.rm=T),
            Anemones=mean(Anemones, na.rm=T), Tunicates=mean(Sea.squirt, na.rm=T))

## Compare abundances among treatments
round(mean(abundanceszl[abundanceszl$Zone=="I",]$Barnacles)/mean(abundanceszl[abundanceszl$Zone=="S",]$Barnacles), 2)
round(abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV3",]$Barnacles/abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV1",]$Barnacles, 2)
round(abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV3",]$Barnacles/abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV2",]$Barnacles, 2)

round(abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV3",]$Mussels/abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV1",]$Mussels, 2)
round(abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV3",]$Mussels/abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV1",]$Mussels, 2)
round(abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV3",]$Mussels/abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV2",]$Mussels, 2)
round(abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV3",]$Mussels/abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV2",]$Mussels, 2)

round(abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV3",]$Tunicates/abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV1",]$Tunicates, 2)
round(abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV2",]$Tunicates/abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV2",]$Tunicates, 2)
round(abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV3",]$Tunicates/abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV3",]$Tunicates, 2)

round(abundanceszl[abundanceszl$Zone=="S"&abundanceszl$Level=="LV2",]$Anemones/abundanceszl[abundanceszl$Zone=="I"&abundanceszl$Level=="LV2",]$Anemones, 2)

#--------------------------
## Analyse algae cover per treatment

algaemodel <- glmmTMB(data=pyramid_cover, Algae~Level*Zone+(1|Platform), family="gaussian")
# hist(resid(algaemodel))
summary(algaemodel)
car::Anova(algaemodel)
algaemeans <- emmeans(algaemodel, pairwise~Level*Zone)
summary(contrast(algaemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr")

palg <- data.frame(summary(contrast(algaemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr"))
palg$labels <- symnum(palg$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

cldalg <- data.frame(cld(object=algaemeans, adjust = "fdr", alpha=0.05, Letters = LETTERS))%>%
  rename(group = .group)

E <- ggplot(data=pyramid_cover, aes(x=Level, y=Algae, fill=Zone))+
  geom_boxplot(outlier.shape = NA)+
  labs(x= "", y="Algae coverage")+ 
  geom_point(position=position_jitterdodge())+
  ylim(-0.02, 0.7)+
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  theme_classic()+
  theme(text=element_text(size=15))+
  geom_signif(y_position = c(0.67, 0.43, 0.35), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = palg$labels[c(7,8,9)])+
  geom_text(data=cldalg, aes(y = c(-0.01, -0.01, -0.01, -0.01, -0.01, -0.01), label = group), position = position_dodge(width = .75))

#--------------------------
## Effect sizes

## Calculate cover per treatment
coverzl<- pyramid_cover %>%
  group_by(Zone, Level) %>%
  summarise(Algae=mean(Algae, na.rm=T))

## Compare cover among treatments
round(coverzl[coverzl$Zone=="S"&coverzl$Level=="LV1",]$Algae/coverzl[coverzl$Zone=="S"&coverzl$Level=="LV2",]$Algae, 2)
round(coverzl[coverzl$Zone=="S"&coverzl$Level=="LV2",]$Algae/coverzl[coverzl$Zone=="S"&coverzl$Level=="LV3",]$Algae, 2)
round(coverzl[coverzl$Zone=="S"&coverzl$Level=="LV1",]$Algae/coverzl[coverzl$Zone=="I"&coverzl$Level=="LV1",]$Algae, 2)
round(coverzl[coverzl$Zone=="I"&coverzl$Level=="LV3",]$Algae/coverzl[coverzl$Zone=="S"&coverzl$Level=="LV3",]$Algae, 2)

#--------------------------
## Generate the plot for abundances & cover from separate panels

(A+plot_layout(guides="collect")& theme(legend.position='left'))+ (B+theme(legend.position='left')) +C +D +E+guide_area()+plot_layout(ncol=2, axis_title="collect")+
  plot_annotation(tag_levels = 'A')
#-------------------------------------------------------------------
## Accumulation curves

## Merge the abundances and cover to get all species per sample
## Group per sample
## Round observations for analyses
total <- merge(abundances_corrected%>%
               dplyr::select(-c(Hydroid.polyps))%>%
               group_by(Location,Zone, Level, Platform, Vertical, Horizontal)%>%  
               summarise(across(Barnacles:Ensis.sp., sum)), 
               cover_corrected%>%
               dplyr::select(-c(Small.anemones, Tube.structures))%>%
               group_by(Location,Zone, Level, Platform, Vertical, Horizontal)%>%  
               summarise(across(c(Algae:Hydroid.polyps), sum)),
               by=c("Location", "Zone", "Level", "Platform", "Vertical", "Horizontal"))%>%
               mutate(across(c(Barnacles:Hydroid.polyps),ceiling))
               # mutate(across(c(Barnacles:Hydroid.polyps), ~ ifelse(. > 0, 1, 0)))

## Accumulation for Low complexity intertidal
total1i <- total[total$Location=="K" & total$Zone=="I" & total$Level=="LV1",] %>%
  dplyr::select(-c("Location", "Zone", "Level", "Platform", "Vertical", "Horizontal")) # subset data
sp1i <- specaccum(total1i, method="exact") # accumulate species
ac1i <- data.frame(ID="LV1_I", Sites=sp1i$sites, Richness=sp1i$richness, SD=sp1i$sd) # save as data frame
m1im <- nls(Richness~Vm*Sites/(K+Sites), data=ac1i, start=list(K=0.01, Vm=20)) # fit Michaelis-Menten accumulation curve
a11 <- data.frame(Zone="I", LV="LV1", half=coef(m1im)[[1]]) # save values of fit for plot

## Accumulation for Low complexity subtidal
total1s <- total[total$Location=="K" & total$Zone=="S" & total$Level=="LV1",] %>%
  dplyr::select(-c("Location", "Zone", "Level", "Platform", "Vertical", "Horizontal"))
sp1s <- specaccum(total1s, method="exact")
ac1s <- data.frame(ID="LV1_S", Sites=sp1s$sites, Richness=sp1s$richness, SD=sp1s$sd)
m1sm <- nls(Richness~Vm*Sites/(K+Sites), data=ac1s, start=list(K=0.01, Vm=20))
a1s <- data.frame(Zone="S", LV="LV1", half=coef(m1sm)[[1]])

## Accumulation for Mid complexity intertidal
total2i <- total[total$Location=="K" & total$Zone=="I" & total$Level=="LV2",] %>%
  dplyr::select(-c("Location", "Zone", "Level", "Platform", "Vertical", "Horizontal"))
sp2i <- specaccum(total2i, method="exact")
ac2i <- data.frame(ID="LV2_I", Sites=sp2i$sites, Richness=sp2i$richness, SD=sp2i$sd)
m2im <- nls(Richness~Vm*Sites/(K+Sites), data=ac2i, start=list(K=0.01, Vm=20)) 
a2i <- data.frame(Zone="I", LV="LV2", half=coef(m2im)[[1]])

## Accumulation for Mid complexity subtidal
total2s <- total[total$Location=="K" & total$Zone=="S" & total$Level=="LV2",] %>%
  dplyr::select(-c("Location", "Zone", "Level", "Platform", "Vertical", "Horizontal"))
sp2s <- specaccum(total2s, method="exact")
ac2s <- data.frame(ID="LV2_S", Sites=sp2s$sites, Richness=sp2s$richness, SD=sp2s$sd)
m2sm <- nls(Richness~Vm*Sites/(K+Sites), data=ac2s, start=list(K=0.01, Vm=20))
a2s <- data.frame(Zone="S", LV="LV2", half=coef(m2sm)[[1]])

## Accumulation for High complexity intertidal
total3i <- total[total$Location=="K" & total$Zone=="I" & total$Level=="LV3",] %>%
  dplyr::select(-c("Location", "Zone", "Level", "Platform", "Vertical", "Horizontal"))
sp3i <- specaccum(total3i, method="exact")
ac3i <- data.frame(ID="LV3_I", Sites=sp3i$sites, Richness=sp3i$richness, SD=sp3i$sd)
m3im <- nls(Richness~Vm*Sites/(K+Sites), data=ac3i, start=list(K=0.01, Vm=20)) 
a3i <- data.frame(Zone="I", LV="LV3", half=coef(m3im)[[1]])

## Accumulation for High complexity subtidal
total3s <- total[total$Location=="K" & total$Zone=="S" & total$Level=="LV3",] %>%
  dplyr::select(-c("Location", "Zone", "Level", "Platform", "Vertical", "Horizontal"))
sp3s <- specaccum(total3s, method="exact")
ac3s <- data.frame(ID="LV3_S", Sites=sp3s$sites, Richness=sp3s$richness, SD=sp3s$sd)
m3sm <- nls(Richness~Vm*Sites/(K+Sites), data=ac3s, start=list(K=0.01, Vm=20)) 
a3s <- data.frame(Zone="S", LV="LV3", half=coef(m3sm)[[1]])

## combune all the stored fitted values for the different treatments
accum <- rbind(ac1i, ac2i, ac3i, ac1s, ac2s, ac3s)%>%
  separate(ID, c("Level", "Zone"), remove = FALSE)%>%
  mutate(Label=case_when(ID=="LV1_S" | ID=="LV1_I" ~ "Low",
                         ID=="LV2_S" | ID=="LV2_I" ~ "Mid",
                         ID=="LV3_S" | ID=="LV3_I" ~ "High"))

## Find number of sites where half the max number of species is reached
halves <- data.frame(half=c(coef(m1im)[[1]], coef(m1sm)[[1]], coef(m2im)[[1]], coef(m2sm)[[1]], coef(m3im)[[1]], coef(m3sm)[[1]]),
                     Zone=c("I", "S", "I", "S","I", "S"), Level=c("LV1", "LV1", "LV2", "LV2", "LV3", "LV3"))

## Compare different accumulation curves
# m1 <- nls(Richness~Vm*Sites/(K+Sites), data=accum[accum$ID=="LV3_I",], start=list(K=0.01, Vm=20)) # Michaelis-Menten; K=value at which half of assymptote is reached
# m2 <- nls(Richness~Vm-Vm*exp(-K*Sites), data=accum[accum$ID=="LV3_I",], start=list(K=0.01, Vm=20)) # Bertalanffy
# m3 <- nls(Richness~Vm*Sites^K, data=accum[accum$ID=="LV3_I",], start=list(K=0.01, Vm=20))# Preston (power-law)
# m4 <- nls(Richness~log(K)+Vm*log(Sites), data=accum[accum$ID=="LV3_I",], start=list(K=0.01, Vm=20)) # Gleason semi-log
# m5 <- nls(Richness~Vm/(1+exp((V0-Sites)/K)), data=accum[accum$ID=="LV3_I",], start=list(K=0.5, Vm=20, V0=5)) # Logistic
# m6 <- nls(Richness~Vm-V.5*exp(-exp(K)*Sites^1), data=accum[accum$ID=="LV3_I",], start=list(K=-2, Vm=20, V.5=5)) # Weibull; Vm is assymptote
# AIC(m1, m2, m3, m4, m5, m6) #best fits:(L1I:m4, L1S:m4,L2I:m4,L2S:m1,L3I:m4, L3S:m1, m6 2e keus) so fit Gleason semilog for now

## Plot for accumulation curves
acplot <- ggplot(data=accum, aes(x=Sites, y=Richness, color=interaction(Zone, Level), linetype=interaction(Zone, Level), fill=interaction(Zone, Level)))+
  geom_ribbon(aes(ymin=(Richness-2*SD),ymax=(Richness+2*SD)),alpha=0.2, colour=NA , show.legend = FALSE)+
  geom_smooth(method="nls", method.args=list(formula=y~Vm*x/(K+x), start=list(K=2, Vm=20)), data=accum[accum$ID=="LV1_I",], se=F)+
  geom_smooth(method="nls", method.args=list(formula=y~Vm*x/(K+x), start=list(K=2, Vm=20)), data=accum[accum$ID=="LV1_S",], se=F)+
  geom_smooth(method="nls", method.args=list(formula=y~Vm*x/(K+x), start=list(K=2, Vm=20)), data=accum[accum$ID=="LV2_I",], se=F)+
  geom_smooth(method="nls", method.args=list(formula=y~Vm*x/(K+x), start=list(K=2, Vm=20)), data=accum[accum$ID=="LV2_S",], se=F)+
  geom_smooth(method="nls", method.args=list(formula=y~Vm*x/(K+x), start=list(K=2, Vm=20)), data=accum[accum$ID=="LV3_I",], se=F)+
  geom_smooth(method="nls", method.args=list(formula=y~Vm*x/(K+x), start=list(K=2, Vm=20)), data=accum[accum$ID=="LV3_S",], se=F)+
  geom_vline(data=halves, size=1, aes(xintercept=half, color = interaction(factor(Zone), factor(Level)), lty=interaction(factor(Zone), factor(Level))), show.legend = FALSE)+
  scale_fill_manual(values = c(rep(c("#F2AF4AFF","#3B7C70FF"), times=3),"#F2AF4AFF","#3B7C70FF") , guide="none")+
  scale_colour_manual(values = c(rep(c("#F2AF4AFF","#3B7C70FF"), times=3),"#F2AF4AFF","#3B7C70FF"),
                      labels = c("Low Intertidal", "Low Subtidal", "Mid Intertidal", "Mid Subtidal", "High Intertidal", "High Subtidal"))+
  scale_linetype_manual(values=c("solid", "solid", "dashed", "dashed", "dotted", "dotted", "solid", "dashed", "dotted"),
                        labels = c("Low Intertidal", "Low Subtidal", "Mid Intertidal", "Mid Subtidal", "High Intertidal", "High Subtidal"))+
  geom_label(data = accum %>% filter(Sites == last(Sites) & ID!= "LV2_S"), aes(label = Label,
              x = Sites+2, y = Richness, fill=interaction(Zone, Level)), colour="black", fontface="bold", size=6)+
  geom_label(data = accum %>% filter(Sites == last(Sites) & ID== "LV2_S"), aes(label = Label,
              x = Sites+2,y = Richness+1, fill=interaction(Zone, Level)), colour="black", fontface="bold", size=6)+
  xlab("Samples")+
  xlim(c(0,52.5))+
  labs(colour="Treatment", linetype="Treatment")+
  theme_classic()+
  theme(text=element_text(size=30),legend.key.size=unit(2,"lines"), legend.title = element_text(size = 12), legend.text = element_text(size = 10), legend.position="none")
       
#-------------------------------------------------------------------
## Taxonomic richness

## Combine abundance and cover dataframes
## Change abundances to presence/absence
## Sum per structure
## Calculate taxonomic richness
richness <- merge(pyramid_abundance%>%
                    mutate(across(c(Barnacles:Ensis.sp.), ~ ifelse(. > 0, 1, 0))),
                  pyramid_cover%>%
                    dplyr::select(-c("Small.anemones", "Tube.structures"))%>%
                    mutate(across(c(Algae:Hydroid.polyps), ~ ifelse(. > 0, 1, 0))),
                  by=c("Zone", "Level", "Platform"))%>%
  mutate(Richness =rowSums(across(c(Barnacles:Hydroid.polyps)))) 

## Compare richness among treatments
richnessmodel <- glmmTMB(data=richness, Richness~Level*Zone+(1|Platform), family="poisson")
car::Anova(richnessmodel)
richnessmeans <- emmeans(richnessmodel, pairwise~Level*Zone)
summary(contrast(richnessmeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr")

prich <- data.frame(summary(contrast(richnessmeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr"))
prich$labels <- symnum(prich$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

cldrich <- as.data.frame(cld(object=richnessmeans, adjust = "fdr", alpha=0.05, Letters = LETTERS))%>%
  rename(group = .group)

## Plot for richness
richplot <- ggplot(richness, aes(x=Level, y=Richness, fill=Zone))+
  geom_boxplot()+
  geom_point(position=position_jitterdodge())+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  labs(x= "Complexity", y="Taxonomic Richness")+ 
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  theme_classic()+
  geom_signif(y_position = c(15, 15, 15), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = prich$labels[c(7,8,9)], textsize=8)+
  geom_text(data=cldrich, aes(y = c(2, 5, 7, 9, 10, 10), label = group), position = position_dodge(width = .75), size=6)+
  theme(text=element_text(size=30))

#--------------------------
## Effect sizes

## Calculate richness per treatment
richnesszl<- richness %>%
  group_by(Zone, Level) %>%
  summarise(richness=mean(Richness, na.rm=T))

## Compare richness among treatments
round(richnesszl[richnesszl$Zone=="I"&richnesszl$Level=="LV2",]$richness/richnesszl[richnesszl$Zone=="I"&richnesszl$Level=="LV1",]$richness, 2)
round(richnesszl[richnesszl$Zone=="I"&richnesszl$Level=="LV3",]$richness/richnesszl[richnesszl$Zone=="I"&richnesszl$Level=="LV1",]$richness, 2)
round(richnesszl[richnesszl$Zone=="S"&richnesszl$Level=="LV1",]$richness/richnesszl[richnesszl$Zone=="I"&richnesszl$Level=="LV1",]$richness, 2)

#--------------------------
## Generate the plot for accumulation and richness from separate panels

(richplot + plot_layout(guides = "collect") & theme(legend.position = "bottom"))+acplot+plot_annotation(tag_levels = 'A')

#-------------------------------------------------------------------
## Calculate densities per structure

## Surface area per complexity level
# LV3 --> 144 sides (36 per side)
# Area 21892.78
# 21892.78/144=152.03cm2 per side (16 'sides' sampled)
# sampled surface: 2432.531cm2
# LV2 --> 144 sides (36 per side)
# Area 9730.13cm2
# 9730.13/144=67.57cm2 per side (16 'sides' sampled)
# sampled surface: 1081.126cm2
# LV1 --> 64 sides (16 per side)
# Area 4324.5	cm2
# 4324.5/64=67.57cm2 per side (10 'sides' sampled)
# sampled surface: 675.7cm2

surface_abundance <- abundances_corrected%>% 
  dplyr::select(-c("Seaweeds", "Hydroid.polyps"))%>%
  dplyr::filter(Side != "Residue")%>%
  group_by(Zone, Level, Platform)%>%
  summarise(across(Barnacles:Ensis.sp., sum))%>%
  rowwise()%>%
  mutate(across(c(Barnacles:Ensis.sp.), ~case_when((Level== "LV3")~ ./2432.5,
                                                   (Level== "LV2")~ ./1081.1,
                                                   (Level== "LV1")~ ./675.7, TRUE~.)))%>%
  ungroup()

#--------------------------
## Compare densities among treatments

## Barnacles
barnaclemodel <- glmmTMB(data=surface_abundance, Barnacles~Level*Zone+(1|Platform), family="gaussian")
# hist(resid(barnaclemodel))
# qqnorm(resid(barnaclemodel))
# qqline(resid(barnaclemodel))
car::Anova(barnaclemodel)
barnaclemeans <- emmeans(barnaclemodel, pairwise~Level*Zone)
summary(contrast(barnaclemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr")

pb <- data.frame(summary(contrast(barnaclemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr"))
pb$labels <- symnum(pb$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

cldb <- data.frame(cld(object=barnaclemeans, adjust = "fdr", alpha=0.05, Letters = LETTERS))%>%
  rename(group = .group)

a <- ggplot(data=surface_abundance, aes(x=Level, y=Barnacles, fill=Zone))+
  geom_boxplot(outlier.shape = NA, position='dodge')+
  geom_point(position=position_jitterdodge())+
  labs(x= "", y=expression(paste("Barnacle Density (n/cm"^2,")")))+ 
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  theme_classic()+
  theme(text=element_text(size=15))+
  geom_signif(y_position = c(1.8, 1.6, 1.2), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = pb$labels[c(7,8,9)])+
  geom_text(data=cldb, aes(y = c(-0.1, -0.1, -0.1, -0.1, -0.1,  -0.1), label = group), position = position_dodge(width = .75))

#--------------------------
## Mussels

musselmodel <- glmmTMB(data=surface_abundance, Mussels~Level*Zone+(1|Platform), family="gaussian")
car::Anova(musselmodel)
musselmeans <- emmeans(musselmodel, pairwise~Level*Zone)
summary(contrast(musselmeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr") 

pm <- data.frame(summary(contrast(musselmeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="none"))
pm$labels <- symnum(pm$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

cldm <- data.frame(cld(object=musselmeans, adjust = "none", alpha=0.05, Letters = LETTERS))%>%
  rename(group = .group)

b <-  ggplot(data=surface_abundance, aes(x=Level, y=Mussels, fill=Zone))+
  geom_boxplot(outlier.shape = NA)+
  labs(x= "", y=expression(paste("Mussel Density (n/cm"^2,")")))+ 
  geom_point(position=position_jitter())+
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  theme_classic()+
  theme(text=element_text(size=15), legend.position = "none")+
  geom_signif(y_position = c(0.006, 0.006, 0.016), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = pm$labels[c(7,8,9)])+
  geom_text(data=cldm, aes(y = c(-0.002, -0.002, -0.002, -0.002, -0.002,  -0.002), label = group), position = position_dodge(width = .75))

#--------------------------
## Tunicates
tunicatemodel <- glmmTMB(data=surface_abundance, Sea.squirt~Level*Zone+(1|Platform), family="gaussian")
car::Anova(tunicatemodel) 
tunicatemeans <- emmeans(tunicatemodel, pairwise~Level*Zone)
summary(contrast(tunicatemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr")

pt <- data.frame(summary(contrast(tunicatemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr"))
pt$labels <- symnum(pt$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

cldt <- as.data.frame(cld(object=tunicatemeans, adjust = "fdr", alpha=0.05, Letters = LETTERS))%>%
  rename(group = .group)

c <-   ggplot(data=surface_abundance, aes(x=Level, y=Sea.squirt, fill=Zone))+
  geom_boxplot(outlier.shape = NA)+
  geom_point(position=position_jitterdodge())+
  labs(x= "", y=expression(paste("Tunicate Density (n/cm"^2,")")))+ 
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  theme_classic()+
  theme(text=element_text(size=15), legend.position = "none")+
  geom_signif(y_position = c(0.017, 0.021, 0.012), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = pt$labels[c(7,8,9)])+
  geom_text(data=cldt, aes(y = c(-0.002, -0.002, -0.002, -0.002, -0.002,  -0.002), label = group), position = position_dodge(width = .75))

#--------------------------
## Anemones
anemonemodel <- glmmTMB(data=surface_abundance, Anemones~Level*Zone+(1|Platform), family="gaussian")
car::Anova(anemonemodel)
anemonemeans <- emmeans(anemonemodel, pairwise~Level*Zone)
summary(contrast(anemonemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="fdr")
summary(contrast(anemonemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="none") # more lenient correction, otherwise ns

pa <- data.frame(summary(contrast(anemonemeans, "pairwise")[c(1,2,6,13,14,15,3,8,12)], adjust="none"))
pa$labels <- symnum(pa$p.value, corr = FALSE, cutpoints = c(0,  .001,.01,.05, 1), symbols = c("***","**","*","n.s."))

clda <- data.frame(cld(object=anemonemeans, adjust = "none", alpha=0.01, Letters = LETTERS))%>%
  rename(group = .group)

d <-  ggplot(data=surface_abundance, aes(x=Level, y=Anemones, fill=Zone))+
  geom_boxplot(outlier.shape = NA)+
  labs(x= "", y=expression(paste("Anemone Density (n/cm"^2,")")))+ 
  geom_point(position=position_jitterdodge())+
  scale_fill_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  scale_x_discrete(labels=c("LV1" = "Low", "LV2" = "Mid", "LV3" = "High"))+
  theme_classic()+
  theme(text=element_text(size=15),legend.position = "none")+
  geom_signif(y_position = c(0.018, 0.006, 0.003), xmin = unique(as.numeric(pyramid_abundance$Level))-.2, xmax = unique(as.numeric(pyramid_abundance$Level))+.2, 
              annotations = pa$labels[c(7,8,9)])+
  geom_text(data=clda, aes(y = c(-0.001, -0.001, -0.001, -0.001, -0.001, -0.001), label = group), position = position_dodge(width = .75))

#--------------------------
## Effect sizes

## Calculate densities per treatment
surfacezl<- surface_abundance %>%
  group_by(Zone, Level) %>%
  summarise(Barnacles=mean(Barnacles, na.rm=T), Mussels=mean(Mussels, na.rm=T),
            Tunicates=mean(Sea.squirt, na.rm=T), Anemones=mean(Anemones, na.rm=T))

## Compare densities among treatments
round(surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV1",]$Barnacles/surfacezl[surfacezl$Zone=="S"&surfacezl$Level=="LV1",]$Barnacles, 2)
round(surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV2",]$Barnacles/surfacezl[surfacezl$Zone=="S"&surfacezl$Level=="LV2",]$Barnacles, 2)
round(surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV3",]$Barnacles/surfacezl[surfacezl$Zone=="S"&surfacezl$Level=="LV3",]$Barnacles, 2)

round(surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV3",]$Mussels/surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV2",]$Mussels, 2)
round(surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV3",]$Mussels/surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV1",]$Mussels, 2)

# round(surfacezl[surfacezl$Zone=="S"&surfacezl$Level=="LV1",]$Tunicates/surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV1",]$Tunicates, 2)
round(surfacezl[surfacezl$Zone=="S"&surfacezl$Level=="LV2",]$Tunicates/surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV2",]$Tunicates, 2)
round(surfacezl[surfacezl$Zone=="S"&surfacezl$Level=="LV3",]$Tunicates/surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV3",]$Tunicates, 2)

round(surfacezl[surfacezl$Zone=="S"&surfacezl$Level=="LV2",]$Anemones/surfacezl[surfacezl$Zone=="I"&surfacezl$Level=="LV2",]$Anemones, 2)

#--------------------------
## Generate the plot for densities from separate panels

a + b + c + d + plot_annotation(tag_levels = 'A') + plot_layout(ncol=2, guides = 'collect') & theme(legend.position = "bottom") 

#-------------------------------------------------------------------
## Assess positions on structures

## Only select edge of structure
## Group per structure
## Round counts for analyses
## Add plotting variables
vertical_abundance <- abundances_corrected %>%
  dplyr::select(-c("Seaweeds", "Hydroid.polyps"))%>%
  filter(Horizontal=="H1")%>%
  group_by(Location, Zone, Level, Vertical, Platform) %>%
  summarise(across(c(Barnacles:Ensis.sp.), sum))%>%
  mutate(across(c(Barnacles:Ensis.sp.),round))%>%
  ungroup()%>%
  mutate(ID=paste(Zone, Level,Platform, sep="_"))%>%
  mutate(Vertical.plot=rev(as.numeric(Vertical)))%>%
  mutate(Level.plot=factor(case_when(Level=="LV1"~ "Low", 
                              Level=="LV2" ~ "Mid",
                              Level=="LV3" ~ "High"), 
         levels=c("Low", "Mid", "High")))
  
## Only select middle of structure
## Group per structure
## Round counts for analyses
## Add plotting variables
horizontal_abundance <- abundances_corrected %>%
  dplyr::select(-c("Seaweeds", "Hydroid.polyps"))%>%
  filter(Vertical=="V3")%>%
  group_by(Location, Zone, Level, Horizontal, Platform) %>%
  summarise(across(c(Barnacles:Ensis.sp.), sum))%>%
  mutate(across(c(Barnacles:Ensis.sp.),round))%>%
  ungroup()%>%
  mutate(ID=paste(Zone, Level, Platform, sep="_"))%>%
  mutate(Horizontal.plot=as.numeric(Horizontal))%>%
  mutate(Level.plot=factor(case_when(Level=="LV1"~ "Low", 
                                     Level=="LV2" ~ "Mid",
                                     Level=="LV3" ~ "High"), 
                           levels=c("Low", "Mid", "High")))

## Only select middle of structure
## Group per structure
## Add plotting variables
horizontal_cover <- cover_corrected %>%
  dplyr::select(-c("Small.anemones", "Tube.structures"))%>%
  filter( Vertical=="V3") %>%
  group_by(Location, Zone, Level, Horizontal, Platform) %>%
  summarise(across(c(Algae:Hydroid.polyps), mean))%>%
  ungroup()%>%
  mutate(ID=paste(Zone, Level,Platform, sep="_"))%>%
  mutate(Horizontal.plot=as.numeric(Horizontal))%>%
  mutate(Level.plot=factor(case_when(Level=="LV1"~ "Low", 
                                     Level=="LV2" ~ "Mid",
                                     Level=="LV3" ~ "High"), 
                           levels=c("Low", "Mid", "High")))

## Only select edge of structure
## Group per structure
## Add plotting variables
vertical_cover <- cover_corrected %>%
  dplyr::select(-c("Small.anemones", "Tube.structures"))%>%
  filter(Horizontal=="H1") %>%
  group_by(Location, Zone, Level, Vertical, Platform) %>%
  summarise(across(c(Algae:Hydroid.polyps), mean))%>%
  ungroup()%>%
  mutate(ID=paste(Zone, Level,Platform, sep="_"))%>%
  mutate(Vertical.plot=as.numeric(Vertical))%>%
  mutate(Level.plot=factor(case_when(Level=="LV1"~ "Low", 
                                     Level=="LV2" ~ "Mid",
                                     Level=="LV3" ~ "High"), 
                           levels=c("Low", "Mid", "High")))

#--------------------------
## Barnacles

## Vertical positions
i <-  ggplot(data=vertical_abundance, aes(x=Vertical, y=Barnacles, col=Zone))+
  geom_point(position=position_jitter(width = 0.15, height = 0))+
  geom_smooth()+
  coord_flip()+
  scale_colour_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  theme_linedraw()+
  theme(text=element_text(size=15))+
  labs(x= "Vertical position", y="Barnacle abundance")+ 
  facet_grid(~Level.plot)

barnacle.v <- glmmTMB(data=vertical_abundance, Barnacles~Vertical+(1|Zone/Platform/Level),family="nbinom1")
summary(barnacle.v)
car::Anova(barnacle.v)

barnacle.v.means <- emmeans(barnacle.v, ~Vertical)
contrast(barnacle.v.means, "pairwise", adjust="fdr",  simple = "each", combine = TRUE)

## Horizontal positions
hi <- ggplot(data=horizontal_abundance, aes(x=Horizontal, y=Barnacles, col=Zone))+
  geom_point(position=position_jitter(width = 0.15, height = 0))+
  geom_smooth()+
  labs(x="Horizontal position", y="Barnacle abundance")+
  scale_colour_manual(values = c("#F2AF4AFF","#3B7C70FF"), labels=c("intertidal", "subtidal"))+
  theme_linedraw()+
  theme(text=element_text(size=15))+
  scale_x_continuous(breaks = c(1, 2, 3))+
  facet_grid(~Level.plot)

barnacle.h <- glmmTMB(data=horizontal_abundance, Barnacles~Horizontal+(1|Zone/Platform/Level),family="nbinom1")
summary(barnacle.h)
car::Anova(barnacle.h)
barnacle.h.means <- emmeans(barnacle.h, ~Horizontal)
contrast(barnacle.h.means, "pairwise", adjust="fdr",  simple = "each", combine = TRUE)

#--------------------------
## Mussels

## Vertical positions
ii <-   ggplot(data=vertical_abundance, aes(x=Vertical, y=Mussels, col=Zone))+
  geom_point(position=position_jitter(width = 0.15, height = 0))+
  geom_smooth(lty="dashed")+
  theme_linedraw()+
  theme(text=element_text(size=15))+
  scale_colour_manual(values = c("#F2AF4AFF","#3B7C70FF"), guide="none")+
  coord_flip()+
  labs(x= "Vertical position", y="Mussel abundance")+ 
  facet_grid(~Level.plot)

mussel.v <- glmmTMB(data=vertical_abundance, Mussels~Vertical+(1|Zone/Platform/Level),family="nbinom1")
summary(mussel.v)
car::Anova(mussel.v)
# mussel.v.means <- emmeans(mussel.v, ~Level*Zone*Vertical)

## Horizontal positions
hii <- ggplot(data=horizontal_abundance, aes(x=Horizontal, y=Mussels, col=Zone))+
  geom_point(position=position_jitter(width = 0.15, height = 0))+
  geom_smooth(linetype="dashed")+
  labs(x="Horizontal position", y="Mussel abundance")+
  scale_colour_manual(values = c("#F2AF4AFF","#3B7C70FF"), guide="none")+
  scale_x_continuous(breaks = c(1, 2, 3))+
  theme_linedraw()+
  theme(text=element_text(size=15))+
  facet_grid(~Level.plot)

mussel.h <- glmmTMB(data=horizontal_abundance, Mussels~Horizontal+(1|Zone/Platform/Level),family="nbinom1")
summary(mussel.h)
car::Anova(mussel.h)
mussel.h.means <- emmeans(mussel.h, ~Horizontal)
contrast(mussel.h.means, "pairwise", adjust="fdr",  simple = "each", combine = TRUE)

#--------------------------
## Algae

## Vertical positions
iii <- ggplot(data=vertical_cover, aes(x=Vertical, y=Algae, col=Zone))+
  geom_point(position=position_jitter(width = 0.15, height = 0))+
  geom_smooth(linetype="dashed")+
  coord_flip()+
  scale_x_reverse()+
  scale_colour_manual(values = c("#F2AF4AFF","#3B7C70FF"), guide="none")+
  theme_linedraw()+
  theme(text=element_text(size=15))+
  labs(x="Vertical position", y="Algae cover")+
  facet_grid(~Level.plot)

algae.v<- glmmTMB(data=vertical_cover, logit(Algae, percents = FALSE)~Vertical+(1|Level/Zone/Platform),family="gaussian")
summary(algae.v)
car::Anova(algae.v)
# algae.v.means <- emmeans(algae.v, ~Vertical)
# contrast(algae.v.means, "pairwise", adjust="fdr",  simple = "each", combine = TRUE)

## Horizontal positions
hiii <- ggplot(data=horizontal_cover, aes(x=Horizontal, y=Algae, col=Zone))+
  geom_point(position=position_jitter(width = 0.15, height = 0))+
  geom_smooth(linetype="dashed")+
  labs(x="Horizontal position", y="Algae cover")+
  theme_linedraw()+
  scale_colour_manual(values = c("#F2AF4AFF","#3B7C70FF"), guide="none")+
  theme(text=element_text(size=15))+
  scale_x_continuous(breaks = c(1, 2, 3))+
  facet_grid(~Level.plot)

algae.h<- glmmTMB(data=horizontal_cover, logit(Algae, percents = FALSE)~Horizontal+(1|Level/Zone/Platform),family="gaussian")
summary(algae.h)
car::Anova(algae.h) # Level, level:zone, level:horizontal, zone:horizontal significant
# algae.h.means <- emmeans(algae.h, ~Horizontal|Level*Zone)
# contrast(algae.h.means, "pairwise", adjust="none",  simple = "each", combine = TRUE)

#--------------------------
## Effect sizes

## Calculate abundances per vertical position
abundancesvzl<- vertical_abundance %>%
  group_by(Zone, Level, Vertical) %>%
  summarise(Barnacles=mean(Barnacles, na.rm=T), Mussels=mean(Mussels, na.rm=T))

# Calculate abundances per horizontal position
abundanceshzl<- horizontal_abundance %>%
  group_by(Zone, Level, Horizontal) %>%
  summarise(Barnacles=mean(Barnacles, na.rm=T), Mussels=mean(Mussels, na.rm=T))

round(mean(abundancesvzl[abundancesvzl$Zone=="I"&abundancesvzl$Vertical=="V2"|abundancesvzl$Zone=="I"&abundancesvzl$Vertical=="V3",]$Barnacles)/
        mean(abundancesvzl[abundancesvzl$Zone=="I"&abundancesvzl$Vertical=="V4"|abundancesvzl$Zone=="I"&abundancesvzl$Vertical=="V1",]$Barnacles), 2)

round(abundanceshzl[abundanceshzl$Zone=="I"&abundanceshzl$Level=="LV3"&abundanceshzl$Horizontal=="H1",]$Barnacles/
        abundanceshzl[abundanceshzl$Zone=="I"&abundanceshzl$Level=="LV3"&abundanceshzl$Horizontal=="H3",]$Barnacles, 2)
round(abundanceshzl[abundanceshzl$Zone=="I"&abundanceshzl$Level=="LV3"&abundanceshzl$Horizontal=="H2",]$Barnacles/
        abundanceshzl[abundanceshzl$Zone=="I"&abundanceshzl$Level=="LV3"&abundanceshzl$Horizontal=="H3",]$Barnacles, 2)

#--------------------------
## Generate the plots from the panels

## Vertical positions
i / ii / iii + plot_annotation(tag_levels = 'A') + plot_layout(guides = 'collect') & theme(legend.position = "bottom") 

## Horizontal positions
hi / hii / hiii + plot_annotation(tag_levels = 'A') + plot_layout(guides = 'collect') & theme(legend.position = "bottom") 

## YOU MADE IT TO THE END OF THIS SCRIPT!! ###