## Sterre Witte January 2024

## SUBSTRATE ANALYSES

##---------------------------------
## Load packages
library(tidyverse)
library(dplyr)
library(lme4)
library(lmerTest)
library(car)
library(MASS)
library(emmeans)
library(multcomp)
library(multcompView)
library(vegan)
library(viridis)
library(ggpubr)
library(colorBlindness)

##---------------------------------
## Load data sheets
setwd("...~fonder containing data")
data <- read.csv("Benthic species.csv", sep=";")
mobile <- read.csv("Mobile species.csv", sep=";")
source('PairwiseAdonis.R')

##---------------------------------
## Reorganise dataframe
mobiles <- mobile %>%
  filter(Sampling=="aut20" | Sampling=="aut21")
  
mobiles <- mobiles %>% mutate(sampleyear=factor(ifelse(mobiles$Sampling=='aut20', 2020, 2021)),
         block=factor(substr(mobiles$Cage_no, 2, 3)),
         substrate=factor(substr(mobiles$Cage_no, 1, 1)))
mobiles <- mobiles %>% 
         mutate(Part_Wadden_Sea=factor(ifelse(mobiles$block==5 | mobiles$block==6 |mobiles$block==7,
                                   "East", "West")),
         Substrate=factor(case_when(substrate=="C"~"Granite",
                                    substrate=="B"~"BESE",
                                    substrate=="P"~"Pebbles",
                                    substrate=="R"~"Reefs",
                                    substrate=="S"~"Shells",
                                    substrate=="W"~"Wood")))%>%
        dplyr::rename("fieldlabel"="Cage_no")%>%
        mutate_if(is.character,as.factor)

names(mobiles) <- gsub(".", " ", names(mobiles), fixed = TRUE)

sand20middle <- c(145237, 145236, 145209, 145208, 145249, 145238, 145196, 145197,
                  145178, 145179, 145169, 145168, 145265, 145264, 145219, 145218)

sand21 <- data %>%
  filter(substrate=="X" & sampleyear==2021) %>%
  distinct(sample_id, .keep_all = TRUE) %>%
  pull(sample_id)

sand <- c(sand20middle, sand21)
  
## Only keep selected sand IDs and add the two sand samples together 
## And add same species within one year/substrate/block together
cleaned <-data %>%
  group_by(substrate)%>%
  filter(if(any(sample_id %in% sand)) sample_id %in% sand else TRUE)%>%
  ungroup%>%
  filter(!(sampleyear==2020 & fieldlabel=="X7" | sampleyear==2020 & fieldlabel=="X10" | sample_id==152337))
  #remove faulty resampling in wood by sample_id
  
data_grouped <- cleaned%>%
  group_by(sampleyear, species_name, 
           block, substrate, Substrate, fieldlabel, Part_Wadden_Sea)%>%
  summarise(abundance=sum(n_m2))%>%
  ungroup()

# Make wide dataframe
data_wide <- data_grouped %>% 
  dplyr::mutate_if(is.character,as.factor)%>%
  mutate(Substrate=fct_relevel(Substrate,c("BESE","Granite","Pebbles", "Reefs", "Shells", "Wood", "Sand")))%>%
  mutate(sampleyear=factor(sampleyear),
         block=factor(block))  %>%
  pivot_wider(names_from = species_name, values_from = abundance) %>% 
  replace(is.na(.),0)%>%
  group_by(fieldlabel)%>%
  mutate(c=cur_group_id()) %>%
  arrange(c)%>%
  ungroup()%>%
  dplyr::select(-c)

data_wide[,7:length(data_wide)]  <- ceiling(data_wide[,7:length(data_wide)])
data_wide <- as.data.frame(data_wide)
row.names(data_wide) <- paste(data_wide$fieldlabel, data_wide$sampleyear, sep="_")

data_all <- left_join(data_wide, mobiles)%>%
  dplyr::select(-c(Sampling))%>%
  replace(is.na(.),0)%>%
  dplyr::filter(substrate != "X")
row.names(data_all) <- paste(data_all$fieldlabel, data_all$sampleyear, sep="_")

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

# Calculate species richness
data_wide$richness <- rowSums(data_wide[,7:length(data_wide)] > 0)
data_all$richness <- rowSums(data_all[,7:length(data_all)] > 0)

# Benthic species
srich <- glmer(data=data_wide, richness~substrate+sampleyear+(1|block), family=poisson)
srich1 <- glm(data=data_wide, richness~substrate+sampleyear, family=poisson)
AIC(srich)-AIC(srich1)
Anova(srich)
cldrich <- emmeans(object=srich, spec="substrate")
pairs(cldrich, adjust='fdr')
# summary(srich)

cldsp <- data.frame(cld(object=cldrich, adjust = "fdr", alpha=0.05, Letters = LETTERS))
labelsp <- data_wide %>%
  group_by(Substrate, substrate)%>%
  summarize()%>%
  ungroup()%>%
  right_join(., cldsp) %>%
  rename(group = .group)

A <- ggplot(data_wide, aes(x=Substrate, y=richness, fill=Substrate))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(aes(pch=sampleyear), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16), name="Year")+
  scale_fill_manual(values=c("#01665E","#01665E","#01665E","#01665E","#01665E","#01665E","#FFAD65"), guide="none")+
  geom_text(data = labelsp, aes(label = group, y = 0.1), 
              position = position_dodge(width = .5), 
              show.legend = FALSE, size = 7)+
  labs(y="Species Richness")+
  theme_classic()+
  theme(text = element_text(size = 20), plot.margin=margin(5.5, 5.5, 5.5, 16, "points"))

# All species
sricha <- glmer(data=data_all, richness~substrate+sampleyear+(1|block), family=poisson)
sricha1 <- glm(data=data_all, richness~substrate+sampleyear, family=poisson)
AIC(sricha)-AIC(sricha1)
Anova(sricha)
#summary(sricha)

cldricha <- emmeans(object=sricha, spec="substrate")
pairs(cldricha, adjust='fdr')
# summary(srich)

# Compact letters
cldspall <- data.frame(cld(object=cldricha, adjust='fdr', alpha=0.05, Letters = LETTERS))
labelspall <- data_all %>%
  group_by(Substrate, substrate)%>%
  summarize()%>%
  ungroup()%>%
  # summarize(ypos = mean(richness)+ sd(richness) + .30 )%>%
  right_join(., cldspall)%>%
  rename(group = .group)

Aa <- ggplot(data_all, aes(x=Substrate, y=richness))+
  geom_boxplot(outlier.shape=NA, fill="lightgrey")+
  geom_point(aes(pch=sampleyear), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16), name="Year")+
  geom_text(data = labelspall, aes(label = group, y = 0.1), 
            position = position_dodge(width = .5), 
            show.legend = FALSE,  size = 7)+
  labs(y="Species Richness")+
  theme_classic()+
  theme(text = element_text(size = 20), plot.margin=margin(5.5, 5.5, 5.5, 16, "points"))

##---------------------------------
## Shannon's index
# Benthic
data_wide$shannid <- paste(data_wide$fieldlabel, data_wide$sampleyear, sep="_")
shann <- diversity(data_wide[,c(7:(length(data_wide)-2))]) 
data_wide$shannon <- shann[match(names(shann), data_wide$shannid)]

shannmod <- lmer(data=data_wide, shannon~substrate+sampleyear+(1|block))
shannmod1 <- lm(data=data_wide, shannon~substrate+sampleyear)
AIC(shannmod)-AIC(shannmod1)
Anova(shannmod)

labelshann <- data.frame(Substrate=c("BESE","Granite","Pebbles", "Reefs", "Shells", "Wood", "Sand"),
                       group=c("A","A","A","A","A","A", "A"))

B  <- ggplot(data_wide, aes(x=Substrate, y=shannon, fill=Substrate))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(aes(pch=sampleyear), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16))+
  scale_fill_manual(values=c("#01665E","#01665E","#01665E","#01665E","#01665E","#01665E","#FFAD65"), guide="none")+
  geom_text(data = labelshann, aes(label = group, y = 0.1), 
              position = position_dodge(width = .5), 
              show.legend = FALSE, size=7)+ 
  labs(y="Shannon Diversity Index")+
  theme_classic()+
  theme(text = element_text(size = 20), plot.margin=margin(5.5, 5.5, 5.5, 30, "points") )

# All
data_all$shannid <- paste(data_all$fieldlabel, data_all$sampleyear, sep="_")
shanna <- diversity(data_all[,c(7:(length(data_all)-2))]) 
data_all$shannon <- shanna[match(names(shanna), data_all$shannid)]

shannmoda <- lmer(data=data_all, shannon~substrate+sampleyear+(1|block))
shannmoda1 <- lm(data=data_all, shannon~substrate+sampleyear)
AIC(shannmoda)-AIC(shannmoda1)
Anova(shannmoda)
# summary(shannmoda)

labelshanna <- data.frame(Substrate=c("BESE","Granite","Pebbles", "Reefs", "Shells", "Wood"),
                          group=c("A","A","A","A","A","A"))

Ba   <- ggplot(data_all, aes(x=Substrate, y=shannon))+
  geom_boxplot(outlier.shape=NA, fill="lightgrey")+
  geom_point(aes(pch=sampleyear), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16), name="Year")+
  geom_text(data = labelshanna, aes(label = group, y = 0.1), 
            position = position_dodge(width = .5), 
            show.legend = FALSE, size=7)+
  labs(y="Shannon Diversity Index")+
  theme_classic()+
  theme(text = element_text(size = 20), plot.margin=margin(5.5, 5.5, 5.5, 30, "points"))

ggarrange(A,B, labels = c("A", "B"), font.label = list(size = 30), ncol = 2, nrow = 1, common.legend = T, legend = "right")
ggarrange(Aa,Ba, labels = c("A", "B"), font.label = list(size = 30),  ncol = 2, nrow = 1, common.legend = T, legend = "right")

##---------------------------------
## Multivariate Analyses
  
# species <- data_wide[,c(12:140)
data_wide$substrate <- as.factor(data_wide$substrate)
data_wide$block <- as.factor(data_wide$block)
data_wide$Part_Wadden_Sea <- as.factor(data_wide$Part_Wadden_Sea)
rdadata <- data_wide[,c(1:(length(data_wide)-3))]
rdadata[,c(7:length(rdadata))] <- decostand(rdadata[,c(7:length(rdadata))], "hellinger")

rdalab <- rda(rdadata[,c(7:length(rdadata))] ~ substrate + block + sampleyear, data=rdadata)
anova.cca(rdalab, by="term")
summary(rdalab)
pairwise.adonis(rdadata[,c(7:length(rdadata))], rdadata$substrate)
pairwise.adonis(rdadata[,c(7:length(rdadata))], rdadata$block)

colvec <- c("#0B4151", "#007574", "#1D907F", "#48AA88", "#C7E5BE", "#D4F3A3","#FFAD65")

## Species composition substrate types
par(mar = c(5.1, 5.1, 2.1, 2.1))
plot(rdalab, type = "n", xlab="RDA1", ylab="RDA2", cex.lab=1.5, cex.axis=1.5, cex.main=1, cex.sub=2)
ordibar(rdalab, rdadata$substrate, kind="se", conf=0.95, lwd=3, col="black")
ordiellipse(rdalab, rdadata$substrate, col=colvec,kind="se", draw="polygon", conf=0.95, label=F, alpha=0.5)
ordiellipse(rdalab, rdadata$substrate, col="black",kind="se",  conf=0.95, lwd=2, label=F)
with(rdadata, points(rdalab, display = "sites", pch=c(21,24)[sampleyear], col="black",bg=colvec[substrate], lwd=0.75))
# box(); axis(1); axis(2); title(xlab="RDA1  31.5%", ylab="RDA2  27.8%")
with(rdadata, legend("topright", legend = c(levels(Substrate), levels(sampleyear)), pch = c(21,21,21,21,21,21,21,1,2), 
                     bty = "o", col= "black", pt.bg = c(colvec, "black", "black"), cex=1.45,                 
                     title=expression(bold("Substrate type")), title.adj = 0.25))

## All species
data_all$substrate <- as.factor(data_all$substrate)
data_all$block <- as.factor(data_all$block)
data_all$Part_Wadden_Sea <- as.factor(data_all$Part_Wadden_Sea)
rdadata1 <- data_all[,c(1:(length(data_all)-3))]
rdadata1[,c(7:length(rdadata1))] <- decostand(rdadata1[,c(7:length(rdadata1))], "hellinger")

rdaall<- rda(rdadata1[,c(7:length(rdadata1))] ~ substrate + block + sampleyear, data=rdadata1)
summary(rdaall)
anova.cca(rdaall, by="term")

colveca <- c("#0B4151", "#007574", "#1D907F", "#48AA88", "#C7E5BE", "#D4F3A3")

## Species composition substrate types
par(mar = c(5.1, 5.1, 2.1, 2.1))
plot(rdaall, type = "n", xlab="RDA1", ylab="RDA2", cex.lab=1.5, cex.axis=1.5, cex.main=1, cex.sub=1.2)
ordibar(rdaall, rdadata1$substrate, kind="se", conf=0.95, lwd=3, col="black")
ordiellipse(rdaall, rdadata1$substrate, col=colveca, draw="polygon",kind="se",conf=0.95, lwd=0.75, label=F, alpha = 0.75)
ordiellipse(rdaall, rdadata1$substrate, col="black", kind="se", conf=0.95, lwd=2, label=F)
with(rdadata1, points(rdaall, display = "sites", pch=c(21,24)[sampleyear], col="black",bg=colveca[substrate], lwd=0.75))
with(rdadata1, legend("bottomright", legend = c(levels(Substrate)[1:6], levels(sampleyear)), pch = c(21,21,21,21,21,21,1,2), 
                     bty="O", col= "black", pt.bg = c(colveca, "black", "black"), cex=1.75, 
                  title=expression(bold("Substrate type")), title.adj = 0.25))

##---------------------------------
# Variation no. mobile species between substrates

mobiles$richness <- rowSums(mobiles[,3:25] > 0)

mobrich <- glmer(data=mobiles, richness~substrate+sampleyear+(1|block), family=poisson)
mobrich1 <- glm(data=mobiles, richness~substrate+sampleyear, family=poisson)
AIC(mobrich)-AIC(mobrich1)
Anova(mobrich)

meansmob <- emmeans(object=mobrich, spec="substrate")
pairs(meansmob, adjust='fdr')
# summary(srich)

cldmob <- data.frame(multcomp::cld(object=meansmob, adjust='fdr', Letters = LETTERS, alpha=0.05))
labelmob <- data_wide %>%
  group_by(Substrate, substrate)%>%
  summarize()%>%
  right_join(., cldmob)%>%
  rename(group = .group)
  
Am <- ggplot(mobiles, aes(x=Substrate, y=richness))+
  geom_boxplot(outlier.shape=NA, fill="lightgrey")+
  geom_point(aes(pch=sampleyear), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16), name="Year")+
  geom_text(data = labelmob, aes(label = group, y = -0.75), 
            position = position_dodge(width = .5), 
            show.legend = FALSE, size=7)+
  labs(y="Mobile Species Richness")+
  theme_classic()+
  theme(text = element_text(size = 20), plot.margin=margin(5.5, 5.5, 5.5, 10, "points"))

##---------------------------------
## Fish richness and abundance
fish <- mobiles[,c(2, 26:28, 30, 10:15, 21, 22, 25)]
fish$richness <- rowSums(fish[, 6:14] > 0)
fish$abundance <- rowSums(fish[, 6:14])

fishrich <- glmer(data=fish, richness~substrate+sampleyear+(1|block), family=poisson)
fishrich1 <- glm(data=fish, richness~substrate+sampleyear, family=poisson) 
AIC(fishrich)-AIC(fishrich1)
Anova(fishrich)

meansfish <- emmeans(object=fishrich, spec="substrate")
pairs(meansfish, adjust='fdr')
pairs(meansfish, adjust="none")

cldfish <- data.frame(cld(object=meansfish, adjust = "none", Letters = LETTERS, alpha=0.05))
labelfish <- data_wide %>%
  group_by(Substrate, substrate)%>%
  summarize()%>%
  ungroup()%>%
  right_join(., cldfish) %>%
  rename(group = .group)

Bf <-  ggplot(fish, aes(x=Substrate, y=richness))+
  geom_boxplot(outlier.shape=NA, fill="lightgrey")+
  geom_point(aes(pch=sampleyear), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16), name="Year")+
  geom_text(data = labelfish, aes(label = group, y = -0.75), 
            position = position_dodge(width = .5), 
            show.legend = FALSE, size=7)+
  labs(y="Fish Species Richness")+  
  theme_classic()+
  theme(text = element_text(size = 20), plot.margin=margin(5.5, 5.5, 5.5, 20, "points"))

fishabundance <- glmer(data=fish, abundance~substrate+sampleyear+(1|block), family=poisson)
fishabundance1 <- glm(data=fish, abundance~substrate+sampleyear, family=poisson)
AIC(fishabundance)-AIC(fishabundance1)
Anova(fishabundance)

meansfish <- emmeans(object=fishrich, spec="substrate")
pairs(meansfish, adjust='fdr')
pairs(meansfish, adjust="none")

cldfish <- data.frame(cld(object=meansfish, adjust = "none", Letters = LETTERS, alpha=0.05))
labelfish <- data_wide %>%
  group_by(Substrate, substrate)%>%
  summarize()%>%
  ungroup()%>%
  right_join(., cldfish) %>%
  rename(group = .group)

meansfisha <- emmeans(object=fishabundance, spec="substrate")
pairs(meansfisha, adjust="fdr")

cldfisha <- data.frame(cld(object=meansfisha, adjust = "fdr", Letters = LETTERS, alpha=0.05))
labelfisha <- data_wide %>%
  group_by(Substrate, substrate)%>%
  summarize()%>%
  ungroup()%>%
  right_join(., cldfisha) %>%
  rename(group = .group)

Cf <- ggplot(fish, aes(x=Substrate, y=abundance))+
  geom_boxplot(outlier.shape=NA, fill="lightgrey")+
  geom_point(aes(pch=sampleyear), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16), name="Year")+
  geom_text(data = labelfisha, aes(label = group, y = -2), 
            position = position_dodge(width = .5), 
            show.legend = FALSE, size=7)+
  labs(y="Fish Densities")+
  theme_classic()+
  theme(text = element_text(size = 20), plot.margin=margin(5.5, 5.5, 5.5, 10, "points"))


ggarrange(Am, labels = c("A"), font.label = list(size = 30) ,ncol = 1, nrow = 1)
ggarrange(Bf, labels = c("B"), font.label = list(size = 30), ncol = 1, nrow = 1)
ggarrange(Cf, labels = c("C"), font.label = list(size = 30), ncol = 1, nrow = 1)

##---------------------------------
## Shellfish abundance

shellfish <- data_wide[c("fieldlabel", "sampleyear","block", "substrate", "Substrate", "Mytilidae", "Mytilus edulis", "Magallana gigas")]
shellfish <- dplyr::filter(shellfish, substrate != "X")
shellfish$abundance <- rowSums(shellfish[, 6:7])
shellfish$musselabundance <- sqrt(rowSums(shellfish[, 6:7]))
shellfish <- shellfish[shellfish$sampleyear!=2021 | shellfish$fieldlabel!="P9",]

# mean(shellfish$abundance, na.rm=T)
# sd(shellfish$abundance, na.rm=T)
mussels <- lmer(data=shellfish, musselabundance~sampleyear+substrate+(1|block))
mussels1 <- lm(data=shellfish, musselabundance~sampleyear+substrate)
AIC(mussels)-AIC(mussels1)
Anova(mussels)

musselmeans <- emmeans(object=mussels, spec="substrate")
pairs(musselmeans, adjust="fdr")

cldshellfish <- data.frame(cld(object=musselmeans, adjust='fdr', Letters = LETTERS, alpha=0.05))
labelshellfish <- shellfish %>%
  group_by(substrate, Substrate)%>%
  summarize()%>%
  right_join(., cldshellfish) %>%
  rename(group = .group)

# Plot
ggplot(shellfish, aes(x=Substrate, y=musselabundance))+
  geom_boxplot(outlier.shape = NA, fill="lightgrey")+
  geom_point(aes(pch=sampleyear), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16), name="Year")+
  geom_text(data = labelshellfish, aes(label = group, y = -0.75),
            position = position_dodge(width = .5),
            show.legend = FALSE, size=7)+
  labs(y=bquote(sqrt("Mussel Densities") ~ (ind/m^2)))+
  theme_classic()+
  theme(text = element_text(size = 20))

##---------------------------------
## Percentages increase 

allyear <- data_all %>%
  group_by(sampleyear) %>%
  summarise(richness=mean(richness, na.rm=T), shannon=mean(shannon, na.rm=T))

(allyear[allyear$sampleyear==2021,]$richness-allyear[allyear$sampleyear==2020,]$richness)/
  allyear[allyear$sampleyear==2020,]$richness*100

bsubstrate <- data_wide %>%
  group_by(Substrate, substrate) %>%
  summarise(richness=mean(richness, na.rm=T), shannon=mean(shannon, na.rm=T))  

(mean(bsubstrate[bsubstrate$substrate=="W" | bsubstrate$substrate=="S" | bsubstrate$substrate=="C",]$richness)-
    bsubstrate[bsubstrate$substrate=="X" ,]$richness)/bsubstrate[bsubstrate$substrate=="X" ,]$richness*100

allsub <- data_all %>%
  group_by(Substrate, substrate) %>%
  summarise(richness=mean(richness, na.rm=T), shannon=mean(shannon, na.rm=T))

(mean(allsub[allsub$substrate=="W" | allsub$substrate=="S" | allsub$substrate=="C",]$richness)-
    allsub[allsub$substrate=="B" ,]$richness)/allsub[allsub$substrate=="B" ,]$richness*100

(allsub[allsub$substrate=="W",]$richness-mean(allsub[allsub$substrate=="P" | allsub$substrate=="R",]$richness))/
  mean(allsub[allsub$substrate=="P" | allsub$substrate=="R",]$richness)*100

mobileyear <- mobiles %>%
  group_by(sampleyear) %>%
  summarise(richness=mean(richness, na.rm=T))

(mobileyear[mobileyear$sampleyear==2021,]$richness-mobileyear[mobileyear$sampleyear==2020,]$richness)/
  mobileyear[mobileyear$sampleyear==2020,]$richness*100  

mobilesub <- mobiles %>%
  group_by(Substrate, substrate) %>%
  summarise(richness=mean(richness, na.rm=T))

(mean(mobilesub[mobilesub$substrate=="S",]$richness)-mean(mobilesub[mobilesub$substrate!="S" | mobilesub$substrate!="P",]$richness))/
  mean(mobilesub[mobilesub$substrate!="S" | mobilesub$substrate!="P",]$richness)*100

(mean(mobilesub[mobilesub$substrate=="P",]$richness)-mean(mobilesub[mobilesub$substrate=="C" | mobilesub$substrate=="W",]$richness))/
  mean(mobilesub[mobilesub$substrate=="C" | mobilesub$substrate=="W",]$richness)*100

fishyear <-  fish%>%
  group_by(sampleyear) %>%
  summarise(richness=mean(richness, na.rm=T), abundance=mean(abundance, na.rm=T))

(fishyear[fishyear$sampleyear==2021,]$abundance-fishyear[fishyear$sampleyear==2020,]$abundance)/
  fishyear[fishyear$sampleyear==2020,]$abundance*100

fishsub<-  fish%>%
  group_by(Substrate, substrate) %>%
  summarise(richness=mean(richness, na.rm=T), abundance=mean(abundance, na.rm=T))

(mean(fishsub[fishsub$substrate=="S" ,]$richness)-mean(fishsub[fishsub$substrate=="B" | fishsub$substrate=="C" | fishsub$substrate=="W",]$richness))/
  mean(fishsub[fishsub$substrate=="B" | fishsub$substrate=="C" | fishsub$substrate=="W",]$richness)*100

(mean(fishsub[fishsub$substrate=="B" | fishsub$substrate=="C" | fishsub$substrate=="W",]$abundance)-mean(fishsub[fishsub$substrate=="S" | fishsub$substrate=="P",]$abundance))/
  mean(fishsub[fishsub$substrate=="S" | fishsub$substrate=="P",]$abundance)*100

(mean(fishsub[fishsub$substrate=="R" ,]$abundance)-mean(fishsub[fishsub$substrate=="S" | fishsub$substrate=="P",]$abundance))/
  mean(fishsub[fishsub$substrate=="S" | fishsub$substrate=="P",]$abundance)*100


allshell <- shellfish%>%
  group_by(Substrate, substrate)%>%
  summarise(abundance = mean(abundance))

(mean(allshell[allshell$substrate=="W",]$abundance)-mean(allshell[allshell$substrate!="W"| allshell$substrate!="S",]$abundance))/
  mean(allshell[allshell$substrate!="W"| allshell$substrate!="S",]$abundance)*100

##---------------------------------
# Data for species tied to one unique substrate
data_substrates <- bind_rows(mobiles, data_wide)%>%
  dplyr::select("sampleyear", "fieldlabel", "block", "substrate", "Part_Wadden_Sea",
         "Asterias rubens":"Syngnathus sp", "Abra alba":"Teredo navalis")%>%
  group_by(sampleyear, fieldlabel, block, substrate, Part_Wadden_Sea) %>%
  summarise_all(funs(sum(., na.rm = TRUE)))%>%
  pivot_longer(
    cols = "Asterias rubens":"Teredo navalis", 
    names_to = "species_name",
    values_to = "abundance")%>%
  filter(abundance>0)

persub <- data_substrates%>%
  group_by(species_name)%>%
  mutate(count=length(species_name))%>%
  group_by(substrate, species_name)%>%
  dplyr::summarise(substrateabundance=sum(abundance),
                   count=mean(count))%>%
  ungroup()%>%
  pivot_wider(names_from = substrate, values_from = substrateabundance)

perspecies <- data_substrates%>%
  group_by(species_name)%>%
  dplyr::summarise(totalabundance=sum(abundance))

abundances <- right_join(persub, perspecies) %>%
  mutate_at(vars(B:X),list(fraction=~./totalabundance))%>%
  mutate(maxfraction=pmax(B_fraction,C_fraction,P_fraction,
                          R_fraction,S_fraction, W_fraction, X_fraction, na.rm=T))%>%
  filter(count>2)%>%
  arrange(desc(maxfraction))

# write.csv(abundances, 'species_substrate.csv', row.names = F)