## Code for species richness and species abundances of the predation and 
# substrate experiment
# Sterre Witte September 2024

##----------------------------------------------------------------------
## Prepare workspace
library(ggplot2)
library(dplyr)
library(readxl)
library("tidyverse")
library(car)
library(ggpattern)
library(glmmTMB)
library(ggprism)
library(ggbreak)
library(ggridges)
library(ggpubr)

##----------------------------------------------------------------------
## Prepare the data

# Set working directory
setwd("...~folder containing data")

# Read in data 
species <- read.csv("Exclosures_Species.csv", sep=";")
samples <- read.csv("Exclosures_Samples.csv", sep=";")

# Create unique sample_id for each exclosure
samples$Treatment <- gsub(" ", "_", samples$Treatment)
samples$sample_ID <- paste(samples$Location, samples$Plot_number, samples$Treatment)

# merge the sample information and species information
merged <- merge(samples, species, by = "Poskey", all = TRUE)

# merge with habitat requirements
merged <- merge(merged, habitat, by = "Species", all.x = TRUE, all.y = TRUE)

##----------------------------------------------------------------------
## Shellfish abundances per species

# sum abundances per species
species_sum <- merged[!is.na(merged$Species),] %>% 
   group_by(Species) %>%
   summarise(total=length(Species))%>%
   mutate(Species=ifelse(Species == 'Limecola balthica', 'Macoma balthica', Species))%>%
   mutate(Species = as.factor(Species),
          Species= fct_reorder(Species, -total),
          in_out=ifelse(Species=="Mytilus edulis"|Species=="Cerastoderma edule"|
                           Species=="Mya arenaria"|Species=="Veneridae"|
                           Species=="Ensis leei" | Species=="Macoma balthica",
                        "in", "out")) %>%
   filter(Species != "Asterias rubens" & Species != "Carcinus maenas"
          & Species != "Hemigrapsus sanguineus") %>%
   ungroup()

sum(species_sum[species_sum$in_out=="in",]$total)/sum(species_sum$total) #fraction of total individuals that is a common species
textcol <- rep(c("black", "grey"), times=c(6, 12))
 
# plot the total abundances
ggplot(species_sum, aes(y=fct_rev(Species), x=total, fill=in_out))+
   geom_bar(stat="identity", col="black") + 
   # labs(x= "", x="Total abundance") +
   scale_fill_manual(values = colorRampPalette(c("#003C30","#8C510A"))(2), guide="none")+
   theme_classic() +
   labs(x= "Total abundance", x="")+
   scale_x_continuous(limits=c(0, 12550))+
   scale_y_discrete(guide = "prism_bracket",
                    labels = c(expression(bold(Species), colour=)))+
   theme(text=element_text(size=15),
         axis.text.y.right = element_blank(), axis.line.y.right = element_blank(),
         axis.ticks.y.right = element_blank(), axis.title.y.left = element_blank(),
         axis.text.y.left=element_text( colour=c("grey50",  "grey50", "grey50", "grey50", "grey50", "grey50", "grey50", "grey50", "grey50", "grey50",
                  "grey50", "grey50", "black", "black", "black", "black","black","black"), face="italic", family="sans"))
                                       
##----------------------------------------------------------------------
## Gather data per sample

# calculate number of individuals per species per exclosure
per_sample <- merged %>%
  mutate(abundance=ifelse(is.na(Species), 0, 1))%>% # add abundance 0 for empty samples
  mutate(per_m2=abundance/surf_area)%>% # calculate densities
  group_by(Location, Plot_number, Treatment, Species, sample_ID)%>%
  summarise(abundance_m2=sum(per_m2))

# Expand dataset to include if a species was NOT found in a sample
# all unique species
Species <- unique(merged$Species[!is.na(merged$Species)])
# all unique sample_IDs
sample_ID <- unique(samples$sample_ID)
# all possible sample/species combinations
allpossible <- tidyr::crossing(sample_ID, Species)

# combine all data
data <- left_join(allpossible, per_sample, by=c("sample_ID", "Species")) %>%
  separate(sample_ID, c("Location", "Plot_number", "Treatment"), sep=" ", remove=F) %>%
  dplyr::filter(Treatment!="mystery")%>% #LEAVE OUT THE MYSTERY STARFISH FOR NOW
  mutate_if(is.character, as.factor)%>%
  mutate(abundance_m2=replace_na(abundance_m2, 0))%>%
  mutate(abundance_m2=round(abundance_m2))%>%
  mutate(Treatment=factor(Treatment, levels = c("control_sediment", "control_substrate",
                                                "open_sediment", "open_substrate", 
                                                "closed_sediment", "closed_substrate")))%>%
  separate(Treatment, c("Cage", "Sedsub"), sep="_", remove=F)%>%
  mutate(Cage=factor(Cage, levels=c("control", "open", "closed")))%>%
  mutate_if(is.character,as.factor)

##----------------------------------------------------------------------
## Shellfish densities

# M. edulis
mussel <- data[data$Species=="Mytilus edulis",]

# Data for plotting
musselplot <- mussel%>%
  group_by(Treatment, Location, Cage, Sedsub)%>%
  summarise(abundance=mean(abundance_m2), sd=sd(abundance_m2))

# Make the plot
A <- ggplot(musselplot, aes( y=abundance, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+  
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = abundance+sd, ymin = abundance-sd),
                position = position_dodge(.9), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone") +
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "M. edulis", x= "Treatment", y=bquote("Density" ~(ind/m^-2)))+ 
  scale_alpha_manual(values = c(control=0.3, open=1, closed=1), guide="none")+
  scale_y_continuous(limits=c(-50,30000))+
  scale_y_break(c(770, 6000), scales="free", space=0.4)+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# Test for cage effect
musselcontrol <- mussel %>%
  filter(str_detect(Treatment, "open|control"))
nb.musselcontrol <- glmmTMB(data=musselcontrol, abundance_m2~Treatment+(1|Plot_number), family="nbinom1")
car::Anova(nb.musselcontrol)

# remove the control treatment for testing because of the cage-effect
mussel2 <- mussel %>%
  filter(!str_detect(Treatment, "control"))

# A mussel model
nb.musselmodel <- glmmTMB(data=mussel2, abundance_m2~Location*Sedsub*Cage+(1|Plot_number), family="nbinom1")
# check_model(nb.musselmodel)
summary(nb.musselmodel)
car::Anova(nb.musselmodel) 

# Effect sizes
effects_me <- mussel %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(abundance=mean(abundance_m2))
round(mean(effects_me[effects_me$Cage=="closed",]$abundance)/mean(effects_me[effects_me$Cage=="open",]$abundance), 2)
round(mean(effects_me[effects_me$Sedsub=="substrate",]$abundance)/mean(effects_me[effects_me$Sedsub=="sediment",]$abundance), 2)
round(mean(effects_me[effects_me$Location=="subtidal",]$abundance)/mean(effects_me[effects_me$Location=="intertidal",]$abundance), 2)

# C. edule
cockle <- data[data$Species=="Cerastoderma edule",]

# Data for plotting
cockleplot <- cockle%>%
  group_by(Treatment, Location, Cage, Sedsub)%>%
  summarise(abundance=mean(abundance_m2), sd=sd(abundance_m2))

# Make the plot
B <- ggplot(cockleplot, aes( y=abundance, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+  
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = abundance+sd, ymin = abundance-sd),
                position = position_dodge(.9), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone") +
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "C. edule", x= "Treatment", y=bquote("Density" ~(ind/m^-2)))+
  scale_alpha_manual(values = c(control=0.3, open=1, closed=1), guide="none")+
  scale_y_continuous(limits=c(-25,5500))+
  scale_y_break(c(200, 700), scales="free", space=0.4)+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# Test for cage effect
cocklecontrol <- cockle %>%
  filter(str_detect(Treatment, "open|control"))
nb.cocklecontrol <- glmmTMB(data=cocklecontrol, abundance_m2~Treatment+(1|Plot_number), family="nbinom1")
car::Anova(nb.cocklecontrol)

# remove the control treatment for testing because of the cage-effect
cockle2 <- cockle %>%
  filter(!str_detect(Treatment, "control"))

# A cockle model
nb.cocklemodel <- glmmTMB(data=cockle2, abundance_m2~Location*Sedsub*Cage+(1|Plot_number), family="nbinom1")
summary(nb.cocklemodel)
car::Anova(nb.cocklemodel)

# Effect sizes
effects_ce <- cockle %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(abundance=mean(abundance_m2))
round(mean(effects_ce[effects_ce$Cage=="closed",]$abundance)/mean(effects_ce[effects_ce$Cage=="open",]$abundance), 2)
round(mean(effects_ce[effects_ce$Sedsub=="sediment",]$abundance)/mean(effects_ce[effects_ce$Sedsub=="substrate",]$abundance), 2)
round(mean(effects_ce[effects_ce$Location=="intertidal",]$abundance)/mean(effects_ce[effects_ce$Location=="subtidal",]$abundance), 2)

# M. arenaria
mya <- data[data$Species=="Mya arenaria",]

# Data for plotting
myaplot <- mya%>%
  group_by(Treatment, Location, Cage, Sedsub)%>%
  summarise(abundance=mean(abundance_m2), sd=sd(abundance_m2))

# Make the plot
C <- ggplot(myaplot, aes( y=abundance, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+  
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = abundance+sd, ymin = abundance-sd),
                position = position_dodge(.9), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone") +
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "M. arenaria", x= "Treatment", y=bquote("Density" ~(ind/m^-2)))+
  scale_alpha_manual(values = c(control=0.3, open=1, closed=1), guide="none")+
  scale_y_continuous(limits=c(-7,1000))+
  scale_y_break(c(57, 230), scales="free", space=0.4)+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# Test for cage effect
myacontrol <- mya %>%
  filter(str_detect(Treatment, "open|control"))
nb.myacontrol <- glmmTMB(data=myacontrol, abundance_m2~Treatment+(1|Plot_number), family="nbinom1")
car::Anova(nb.myacontrol)

# remove the control treatment for testing because of the cage-effect
mya2 <- mya %>%
  filter(!str_detect(Treatment, "control"))

# A mya model
nb.myamodel <- glmmTMB(data=mya2, abundance_m2~Location*Sedsub*Cage+(1|Plot_number), family="nbinom1")
car::Anova(nb.myamodel)

# Effect sizes
effects_ma <-mya %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(abundance=mean(abundance_m2))
round(mean(effects_ma[effects_ma$Cage=="closed",]$abundance)/mean(effects_ma[effects_ma$Cage=="open",]$abundance), 2)
round(mean(effects_ma[effects_ma$Sedsub=="substrate",]$abundance)/mean(effects_ma[effects_ma$Sedsub=="sediment",]$abundance), 2)

# Veneridae
veneridae <- data[data$Species=="Veneridae",]

# Data for plotting
veneridaeplot <- veneridae%>%
  group_by(Treatment, Location, Cage, Sedsub)%>%
  summarise(abundance=mean(abundance_m2), sd=sd(abundance_m2))

# Make the plot
D <- ggplot(veneridaeplot, aes( y=abundance, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+  
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = abundance+sd, ymin = abundance-sd),
                position = position_dodge(.9), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone") +
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "Veneridae", x= "Treatment", y=bquote("Density" ~(ind/m^-2)))+
  scale_alpha_manual(values = c(control=1, open=0.3, closed=1), guide="none")+
  scale_y_continuous(limits=c(-3.5,460))+
  scale_y_break(c(15, 150), scales="free", space=0.4)+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# Test for cage effect
veneridaecontrol <- veneridae %>%
  filter(str_detect(Treatment, "open")) # because open are only 0; only keep open and look for significant intercept
nb.veneridaecontrol <- glmmTMB(data=veneridaecontrol, abundance_m2~(1|Plot_number), family="nbinom1")
summary(nb.veneridaecontrol)

# remove the open treatment
veneridae2 <- veneridae%>%
  filter(!str_detect(Treatment, "open"))

# A venereidae model
veneridaemodel <- glmmTMB(data=veneridae2, abundance_m2~Location*Sedsub*Cage+(1|Plot_number), family="poisson")
car::Anova(veneridaemodel)

# Effect sizes
effects_ve <- veneridae %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(abundance=mean(abundance_m2))
round(mean(effects_ve[effects_ce$Cage=="closed",]$abundance), 2)
round(mean(effects_ve[effects_ve$Sedsub=="substrate",]$abundance)/mean(effects_ve[effects_ve$Sedsub=="sediment",]$abundance), 2)
round(mean(effects_ve[effects_ve$Location=="subtidal",]$abundance)/mean(effects_ve[effects_ve$Location=="intertidal",]$abundance), 2)

# E.leei
ensis <- data[data$Species=="Ensis leei",]

# Data for plotting
ensisplot <- ensis%>%
  group_by(Treatment, Location, Cage, Sedsub)%>%
  summarise(abundance=mean(abundance_m2), sd=sd(abundance_m2))

# Make the plot
E <- ggplot(ensisplot, aes( y=abundance, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+  
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = abundance+sd, ymin = abundance-sd),
                position = position_dodge(.9), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone") +
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "E. leei", x= "Treatment", y=bquote("Density" ~(ind/m^-2)))+
  scale_alpha_manual(values = c(control=1, open=0.3, closed=1), guide="none")+
  scale_y_continuous(limits=c(-15,600))+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# Test for cage effect
ensiscontrol <- ensis %>%
  filter(str_detect(Treatment, "open|control"))
nb2.ensiscontrol <- glmmTMB(data=ensiscontrol, abundance_m2~Treatment+(1|Plot_number), family="nbinom2")
car::Anova(nb2.ensiscontrol)

# remove the open treatment for testing because of the cage-effect
ensis2 <- ensis %>%
  filter(!str_detect(Treatment, "open"))

# An Ensis model
nb.ensismodel <- glmmTMB(data=ensis2, abundance_m2~Location*Sedsub*Cage+(1|Plot_number), family="nbinom1")
car::Anova(nb.ensismodel)

# Effect sizes
effects_el <- ensis %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(abundance=mean(abundance_m2))
round(mean(effects_el[effects_el$Cage=="closed",]$abundance)/mean(effects_el[effects_el$Cage=="open",]$abundance), 2)
round(mean(effects_el[effects_el$Sedsub=="sediment",]$abundance)/mean(effects_el[effects_el$Sedsub=="substrate",]$abundance), 2)
round(mean(effects_el[effects_el$Location=="subtidal",]$abundance)/mean(effects_el[effects_el$Location=="intertidal",]$abundance), 2)

# L. balthica
Macoma <- data[data$Species=="Limecola balthica",]

# Data for plotting
Macomaplot <- Macoma%>%
  group_by(Treatment, Location, Cage, Sedsub)%>%
  summarise(abundance=mean(abundance_m2), sd=sd(abundance_m2))

# Make the plot
G <- ggplot(Macomaplot, aes( y=abundance, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+  
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = abundance+sd, ymin = abundance-sd),
                position = position_dodge(.9), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone") +
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "M. balthica", x= "Treatment", y=bquote("Density" ~(ind/m^-2)))+
  scale_alpha_manual(values = c(control=1, open=0.3, closed=1), guide="none")+
  scale_y_continuous(limits=c(-9,400))+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position="none")

# Test for cage effect
Macomacontrol <- Macoma %>%
  filter(str_detect(Treatment, "open|control"))
nb.Macomacontrol <- glmmTMB(data=Macomacontrol, abundance_m2~Treatment+(1|Plot_number), family="nbinom1")
car::Anova(nb.Macomacontrol)

# remove the open treatment for testing because of the cage-effect
Macoma2 <- Macoma %>%
  filter(!str_detect(Treatment, "open"))

# A tellin model
nb.Macomamodel <- glmmTMB(data=Macoma2, abundance_m2~Location*Sedsub*Cage+(1|Plot_number), family="nbinom1")
summary(nb.Macomamodel)
car::Anova(nb.Macomamodel)

# Effect sizes
effects_lb <- Macoma %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(abundance=mean(abundance_m2))
round(mean(effects_lb[effects_lb$Cage=="closed",]$abundance)/mean(effects_lb[effects_lb$Cage=="open",]$abundance), 2)
round(mean(effects_lb[effects_lb$Location=="intertidal",]$abundance)/mean(effects_lb[effects_lb$Location=="subtidal",]$abundance), 2)

# Make combined plot
ggarrange(A,B,C,D,E,G, ncol=2, nrow=3, common.legend = T, legend="bottom")

##----------------------------------------------------------------------
## Species richness 

# Select only shellfish
allsamples <- data.frame(sample_ID)
species_richness <- merged %>%
  group_by(sample_ID, Location, Plot_number, Treatment)%>%
  filter(Species!="Asterias rubens" & Species!="Carcinus maenas" & Species!="Hemigrapsus sanguineus")%>%
  summarise(species_richness=length(unique(Species[!is.na(Species)]))) #nas are for empty samples
  
# Calculate species richness per sample
species_richness <- merge(allsamples, species_richness, all.x=T)%>%
  filter(sample_ID!="mystery mystery mystery")%>%
  separate(sample_ID, c("Location", "Plot_number", "Treatment"), sep=" ", remove=F) %>%
  mutate(species_richness=replace_na(species_richness,0))%>%
  mutate(Treatment=as.factor(Treatment))%>%
  dplyr::filter(Treatment!="mystery")%>% #LEAVE OUT THE MYSTERY STARFISH FOR NOW
  mutate(Treatment=factor(Treatment, levels = c("control_sediment", "control_substrate",
                                                "open_sediment", "open_substrate", 
                                                "closed_sediment", "closed_substrate")))%>%
  separate(Treatment, c("Cage", "Sedsub"), sep="_", remove=F)%>%
  mutate(Cage=factor(Cage, levels=c("control", "open", "closed")))
  
# Plot the species richness
ggplot(species_richness, aes(x=Cage, y=species_richness, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_boxplot_pattern(outlier.shape = NA, col="black", pattern_fill = "black",
                       pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+
  geom_point(aes(), size=3, shape=21, position=position_jitterdodge(0.15, 0.45))+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone")+
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  scale_alpha_manual(values = c(control=0.3, open=1, closed=1), guide="none")+
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme_classic() + 
  labs(y="Shellfish species", x="") +
  theme(text=element_text(size=12))

# Test for cage effect
richnesscontrol <- species_richness %>%
  filter(str_detect(Treatment, "open|control"))
p.richnesscontrol <- glmmTMB(data=richnesscontrol, species_richness~Treatment+(1|Plot_number), family="poisson")
car::Anova(p.richnesscontrol)

# remove the contolrol treatment for testing because of the cage-effect
richness2 <- species_richness %>%
  filter(!str_detect(Treatment, "control"))

# A Richness model
richnessmodel <- glmmTMB(data=richness2, species_richness~Location*Sedsub*Cage+(1|Plot_number), family="poisson")
car::Anova(richnessmodel)

# Effect sizes
effects_sr<- species_richness %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(Richness=mean(species_richness))%>%
  filter(Cage!="control")
round(mean(effects_sr[effects_sr$Cage=="closed",]$Richness)/mean(effects_sr[effects_sr$Cage=="open",]$Richness), 2)
round(mean(effects_sr[effects_sr$Cage=="open"&effects_sr$Location=="intertidal",]$Richness)/
        mean(effects_sr[effects_sr$Cage=="open"&effects_sr$Location=="subtidal",]$Richness), 2)
 
##----------------------------------------------------------------------
## Shellfish sizes

# reframe data
abundantspecies <- merged %>%
  mutate(Treatment=as.factor(Treatment), length=as.numeric(length))%>%
  filter(Treatment!="mystery" & !is.na(Species) & !is.na(length))%>%
  filter(Species=="Mytilus edulis"|
         Species=="Cerastoderma edule"|
         Species=="Ensis leei"|
         Species=="Mya arenaria"|
         Species=="Veneridae" |
         Species=="Limecola balthica")%>%
  mutate(cutoff=case_when(
    Species=="Mytilus edulis" ~ 20,
    Species=="Cerastoderma edule" ~ 17,
    Species=="Ensis leei" ~ 35,
    Species=="Mya arenaria" ~ 24,
    Species=="Limecola balthica" ~ 15,
    TRUE ~ NA))%>%
  mutate(Species=ifelse(Species == 'Limecola balthica', 'Macoma balthica', Species))

# Remove outliers
removed <- abundantspecies %>% 
  filter(Species=="Mytilus edulis" & length>20 |
           Species=="Cerastoderma edule" & length>17 |
           Species=="Ensis leei" & length>35 |
           Species=="Mya arenaria" & length>24 |
           Species=="Macoma balthica" & length>15)%>%
  group_by(Species) %>% 
  summarise(n = n()) %>%
  rbind(c(Species="Veneridae", n=0))

# Combine data
abundantspecies <- left_join(abundantspecies, removed, by="Species")
  
# Plot the outliers
p <- ggplot(abundantspecies, aes(x = length, y = Species)) + 
  geom_density_ridges(fill = "#00AFBB", rel_min_height = 0.01, scale = 0.8,
                      jittered_points = TRUE, position = position_raincloud(width=1, height=0.05), #position = position_points_jitter(width = 0.05, height = -0.3),
                      point_size = 2)
cuts <- abundantspecies %>% distinct(Species, cutoff) %>% mutate(number = row_number())
q <- ggplot_build(p)$data[[1]]
density_lines <- q %>%
  group_by(group) %>% 
  filter(density == max(density)) %>% 
  ungroup()
density_lines_complete <- left_join(density_lines, cuts, by = c("group" = "number"))

ggplot(abundantspecies, aes(x = length, y = Species))+
  geom_density_ridges(fill = "#00AFBB", rel_min_height = 0.01, scale = 0.8,
                      jittered_points = TRUE, position = position_raincloud(width=1, height=0.05), #position = position_points_jitter(width = 0.05, height = -0.3),
                      point_size = 2)+
  geom_text(aes(label=n, x=-3, y=Species), size=3, col='red')+
  geom_segment(data = density_lines_complete, aes(x = cutoff, xend = cutoff, y = ymin-0.2, yend = ymin+0.80), color = "red")+
  labs(x="Length (mm)")+
  theme_bw()+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.y = element_blank(), legend.position = "none")

# Calculate the average length per sample
lengths <- merged %>%
  mutate(Treatment=as.factor(Treatment), length=as.numeric(length))%>%
  filter(Treatment!="mystery" & !is.na(Species) & !is.na(length))%>%
  filter(Species=="Mytilus edulis" & length<=20 |
         Species=="Cerastoderma edule" & length<=17 |
         Species=="Ensis leei" & length<=35 |
         Species=="Mya arenaria" & length<=24 |
         Species=="Veneridae" |
         Species=="Limecola balthica" & length<=15)%>%
  separate(Treatment, c("Cage", "Sedsub"), sep="_", remove=F)%>%
  mutate(Cage=factor(Cage, levels=c("control", "open", "closed")))%>%
  mutate_if(is.character, as.factor)%>%
  group_by(Location, Cage, Sedsub, Plot_number, Species)%>%
  summarise(Average_length=mean(length))

# M. edulis
mussel_size <- lengths[lengths$Species=="Mytilus edulis",] 

# Data for plotting
musselsizeplot <- mussel_size%>%
  group_by(Location, Cage, Sedsub)%>%
  summarise(length=mean(Average_length), sd=sd(Average_length))

# Test for cage effect
mussel_sizecontrol <- mussel_size %>%
  filter(Cage=="open" | Cage=="control")
musselsizecontrol <- glmmTMB(data=mussel_sizecontrol, Average_length~Cage+(1|Plot_number), family="gaussian")
car::Anova(musselsizecontrol)

# remove the open treatment for testing
mussel_size2 <- mussel_size %>%
  filter(Cage!="open")

# Make the plot
a <- ggplot(musselsizeplot, aes(y=length, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(preserve="single"), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+  
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = length+sd, ymin = length -sd),
                 position = position_dodge(.9, preserve="single"), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone") +
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "M. edulis", x= "Treatment", y="Length (mm)")+
  scale_alpha_manual(values = c(control=1, open=0.3, closed=1), guide="none")+
  scale_y_continuous(limits=c(0,7))+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# Mussel size model
musselsize <- glmmTMB(data=mussel_size2, Average_length~Location*Sedsub*Cage+(1|Plot_number), family="gaussian")
car::Anova(musselsize)

# Effect sizes
effects_ms <- mussel_size %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(Length=mean(Average_length))
round(mean(effects_ms[effects_ms$Cage=="closed",]$Length),2)-round(mean(effects_ms[effects_ms$Cage=="control",]$Length), 2)
round(mean(effects_ms[effects_ms$Cage=="closed",]$Length)/mean(effects_ms[effects_ms$Cage=="control",]$Length), 2)

# C. edule 
cockle_size <- lengths[lengths$Species=="Cerastoderma edule",]

# Data for plotting
cocklesizeplot <- cockle_size%>%
  group_by(Location, Cage, Sedsub)%>%
  summarise(length=mean(Average_length), sd=sd(Average_length))

# Test for cage effect
cockle_sizecontrol <- cockle_size %>%
  filter(Cage=="open" | Cage=="control")
cocklesizecontrol <- glmmTMB(data=cockle_sizecontrol, Average_length~Cage+(1|Plot_number), family="gaussian")
car::Anova(cocklesizecontrol)

# remove the open treatment
cockle_size2 <- cockle_size %>%
  filter(Cage!="open")

# Make the plot
b <- ggplot(cocklesizeplot, aes(y=length, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(preserve="single"), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = length+sd, ymin = length -sd),
                position = position_dodge(.9, preserve="single"), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone")+
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white")))+
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "C. edule", x= "Treatment", y="Length (mm)")+
  scale_alpha_manual(values = c(control=1, open=0.3, closed=1), guide="none")+
  scale_y_continuous(limits=c(0,15))+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# A cockle size model
cocklesize <- glmmTMB(data=cockle_size2, Average_length~Cage*Sedsub*Location+(1|Plot_number), family="gaussian")
car::Anova(cocklesize)

# Effect sizes
effects_cs <- cockle_size %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(Length=mean(Average_length))
round(mean(effects_cs[effects_cs$Cage=="closed",]$Length),2)-round(mean(effects_cs[effects_cs$Cage=="control",]$Length), 2)

# E. leei  
ensis_size <- lengths[lengths$Species=="Ensis leei",]

# Data for plotting
ensissizeplot <- ensis_size%>%
  group_by(Location, Cage, Sedsub)%>%
  summarise(length=mean(Average_length), sd=sd(Average_length))

# Test for cage effect
ensis_sizecontrol <- ensis_size %>%
  filter(Cage=="open" | Cage=="control")
ensissizecontrol <- glmmTMB(data=ensis_sizecontrol, Average_length~Cage, family="gaussian")
car::Anova(ensissizecontrol) 

# remove the contolrol treatment for testing because of the cage-effect
ensis_size2 <- ensis_size %>%
  filter(Cage!="open")

# Make the plot
c <- ggplot(ensissizeplot, aes(y=length, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(preserve="single"), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = length+sd, ymin = length -sd),
                position = position_dodge(.9, preserve="single"), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone")+
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white")))+
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "E. leei", x= "Treatment", y="Length (mm)")+
  scale_alpha_manual(values = c(control=1, open=0.3, closed=1), guide="none")+
  scale_y_continuous(limits=c(0,35))+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# An ensis size model
ensissize <- glmmTMB(data=ensis_size, Average_length~Cage*Sedsub*Location+(1|Plot_number), family="gaussian")
summary(ensissize)
car::Anova(ensissize)

# Effect sizes
effects_es <- ensis_size %>%
  filter(Cage!="open")%>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(Length=mean(Average_length))
round(mean(effects_es[effects_es$Sedsub=="sediment",]$Length),2)-round(mean(effects_es[effects_es$Sedsub=="substrate",]$Length), 2)
round(mean(effects_es[effects_es$Location=="intertidal" & effects_es$Cage=="control",]$Length),2)-
  round(mean(effects_es[effects_es$Location=="subtidal" & effects_es$Cage=="control",]$Length), 2)

# M. arenaria
mya_size <- lengths[lengths$Species=="Mya arenaria",]

# Data for plotting
myasizeplot <- mya_size%>%
  group_by(Location, Cage, Sedsub)%>%
  summarise(length=mean(Average_length), sd=sd(Average_length))

# Test for cage effect
mya_sizecontrol <- mya_size %>%
  filter(Cage=="open" | Cage=="control")
myasizecontrol <- glmmTMB(data=mya_sizecontrol, Average_length~Cage+(1|Plot_number), family="gaussian")
car::Anova(myasizecontrol)

# remove the open treatment
mya_size2 <- mya_size %>%
  filter(Cage!="open")

# Make the plot
d <- ggplot(myasizeplot, aes(y=length, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(preserve="single"), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = length+sd, ymin = length -sd),
                position = position_dodge(.9, preserve="single"), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone")+
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white")))+
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "M. arenaria", x= "Treatment", y="Length (mm)")+
  scale_alpha_manual(values = c(control=1, open=0.3, closed=1), guide="none")+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# A mya size model
myasize <- glmmTMB(data=mya_size, Average_length~Cage*Sedsub*Location+(1|Plot_number), family="gaussian")
summary(myasize)
car::Anova(myasize)

# Effect sizes
effects_ms <- mya_size %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(Length=mean(Average_length))
round(mean(effects_ms[effects_ms$Cage=="control",]$Length),2)-round(mean(effects_ms[effects_ms$Cage=="closed",]$Length), 2)
round(mean(effects_ms[effects_ms$Sedsub=="sediment",]$Length),2)-round(mean(effects_ms[effects_ms$Sedsub=="substrate",]$Length), 2)
round(mean(effects_ms[effects_ms$Location=="intertidal" & effects_ms$Cage=="control",]$Length),2)-
  round(mean(effects_ms[effects_ms$Location=="subtidal" & effects_ms$Cage=="control",]$Length), 2)

# Veneridae
veneridae_size <- lengths[lengths$Species=="Veneridae",]

# Data for plotting
veneridaesizeplot <- veneridae_size%>%
  group_by(Location, Cage, Sedsub)%>%
  summarise(length=mean(Average_length), sd=sd(Average_length))

# Test for cage effect
veneridae_sizecontrol <- veneridae_size %>%
  filter(Cage=="open") # because open are only 0; only keep open and look for significant intercept
veneridaesizecontrol <- glmmTMB(data=veneridae_sizecontrol, Average_length~(1|Plot_number), family="gaussian")
summary(veneridaesizecontrol)

# remove the open treatment for testing because of the cage-effect
veneridae_size2 <- veneridae_size %>%
  filter(Cage!="control")

# Make the plot
e <- ggplot(veneridaesizeplot, aes(y=length, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(preserve="single"), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+  
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = length+sd, ymin = length-sd),
                position = position_dodge(.9, preserve="single"), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone") +
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white"))) +
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "Veneridae", x= "Treatment", y="Length (mm)")+
  scale_alpha_manual(values = c(control=0.3, open=1, closed=1), guide="none")+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket", drop="false")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# A venus size model
veneridaesize <- glmmTMB(data=veneridae_size2, Average_length~Cage*Sedsub*Location+(1|Plot_number), family="gaussian")
summary(veneridaesize)
car::Anova(veneridaesize)

# Effect sizes
effects_vs <- veneridae_size %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(Length=mean(Average_length))
round(mean(effects_vs[effects_vs$Cage=="closed",]$Length),2)-round(mean(effects_vs[effects_vs$Cage=="open",]$Length), 2)
mean(effects_vs[effects_vs$Location=="subtidal",]$Length)-mean(effects_vs[effects_vs$Location=="intertidal",]$Length)
mean(effects_vs[effects_vs$Sedsub=="substrate"&effects_vs$Cage=="closed",]$Length)-mean(effects_vs[effects_vs$Sedsub=="sediment"&effects_vs$Cage=="closed",]$Length)

# M. balthica
Macoma_size <- lengths[lengths$Species=="Limecola balthica",]

# Data for plotting
Macomasizeplot <- Macoma_size%>%
  group_by(Location, Cage, Sedsub)%>%
  summarise(length=mean(Average_length), sd=sd(Average_length))

# Test for cage effect
Macoma_sizecontrol <- Macoma_size %>%
  filter(Cage=="open" | Cage=="control")
Macomasizecontrol <- glmmTMB(data=Macoma_sizecontrol, Average_length~Cage+(1|Plot_number), family="gaussian")
car::Anova(Macomasizecontrol)

# Remove the open treatment
Macoma_size2 <- Macoma_size %>%
  filter(Cage!="control")

# Make the plot
f <- ggplot(Macomasizeplot, aes(y=length, x=Cage, fill=Location, pattern=Sedsub, alpha=Cage))+
  geom_bar_pattern(stat="identity", position=position_dodge(preserve="single"), col="black", pattern_fill = "black",
                   pattern_key_scale_factor = 0.25, pattern_angle = 45, pattern_spacing = 0.05)+
  geom_errorbar(aes(group = interaction(Sedsub, Location), ymax = length+sd, ymin = length -sd),
                position = position_dodge(.9, preserve="single"), width=0.3)+
  scale_fill_manual(values = colorRampPalette(c("#BF812D", "#01665E"))(2), name="Zone")+
  scale_pattern_manual(values = c("none", "stripe"), name="Substrate", labels=c("sediment", "granite"),
                       guide=guide_legend(override.aes = list(fill="white")))+
  guides(fill = guide_legend(override.aes = list(pattern = 'none')), y=guide_prism_offset_minor()) +
  theme_classic()+
  labs(title = "L. balthica", x= "Treatment", y="Length (mm)")+
  scale_alpha_manual(values = c(control=1, open=0.3, closed=1), guide="none")+
  scale_y_continuous(limits=c(0,12.5))+
  scale_x_discrete(labels=c("control", "open cage", "closed cage"), guide = "prism_bracket")+
  theme(text=element_text(size=15), axis.text.y.right = element_blank(),
        axis.line.y.right = element_blank(), plot.title = element_text(size = 15, face = "italic"),
        axis.ticks.y.right = element_blank(), axis.title.x = element_blank(), legend.position = "none")

# A macoma size model
Macomasize <- glmmTMB(data=Macoma_size2, Average_length~Cage*Sedsub*Location+(1|Plot_number), family="gaussian")
summary(Macomasize)
car::Anova(Macomasize)

# Effect sizes
effects_ls <- Macoma_size %>%
  group_by(Location, Sedsub, Cage) %>%
  summarise(Length=mean(Average_length))
round(mean(effects_ls[effects_ls$Sedsub=="sediment",]$Length),2)-round(mean(effects_ls[effects_ls$Sedsub=="substrate",]$Length), 2)

# Make the combined plot
ggarrange(a, b, c, d, e, f, ncol=2, nrow=3, common.legend = T, legend="bottom")