## Sterre Witte January 2024

##---------------------------------
## Packages etc.
setwd("...~fonder containing data")

library(tidyr)
library(dplyr)
library(ggplot2)
library(igraph)
library(NetIndices)
library(viridis)
library(ggnetwork)
library(forcats)
library(lme4)
library(lmerTest)
library(multcomp)
library(multcompView)
library(emmeans)
library(car)
library(ggpubr)
source("Foodwebmetrics_functions.R")

##---------------------------------
## Benthos and mobile species data

## Field
mobilesall <- read.csv("Mobile species.csv", head=T, sep=";", check.names=FALSE, fileEncoding="UTF-8-BOM")
mobilesall$sampleyear <- ifelse(grepl("20", mobilesall$Sampling), 2020, 2021)
mobiles <- mobilesall %>%
  dplyr::filter(Sampling=="aut20" | Sampling=="aut21")

uniquefield <- mobiles %>%
  gather(species_name, abundance, "Asterias rubens":"Syngnathus sp", factor_key=TRUE) %>%
  filter(abundance > 0)%>%
  mutate(mobile=1, benthic=0, Samping=NULL, species_name=as.character(species_name), abundance=as.numeric(abundance))%>%
  rename(fieldlabel=Cage_no)%>%
  dplyr::select(sampleyear, fieldlabel, species_name, abundance, mobile, benthic)

## Lab
lab <- read.csv("Benthic species.csv", sep=";", head=T, fileEncoding="UTF-8-BOM")

## Select only 2 sand samples from the middle of each block to correspond with the surface of the cage samples
sand20middle <- c(145237, 145236, 145209, 145208, 145249, 145238, 145196, 145197,
                  145178, 145179, 145169, 145168, 145265, 145264, 145219, 145218)

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

sand <- c(sand20middle, sand21)

## Find the unique species in all the different cages
uniquelab<-lab %>%
  group_by(substrate)%>%
  filter(if(any(sample_id %in% sand)) sample_id %in% sand else TRUE)%>%  #Ik begrijp deze filter nog niet helemaal eerlijk gezegd, maar hij werkt
  ungroup %>%
  filter(!(sampleyear==2020 & fieldlabel=="X7" | sampleyear==2020 & fieldlabel=="X10" | sample_id==152337)) %>%
  group_by(sampleyear, fieldlabel, species_name)%>%
  dplyr::summarise(abundance=sum(n_m2)) %>%
  mutate(mobile=0, benthic=1)

speciesobserved <- rbind(uniquelab, uniquefield)
speciesobserved$mobile <- 1

##---------------------------------
## Food web data

## Interaction matrix from the NIOZ dataverse DOI: 10.25850/nioz/7b.b.gg

interactions <- read.csv("~...folder containing file/Interaction Matrix.csv", head=T, sep= ";",  check.names = FALSE)
rownames(interactions) <-interactions[,1]
interactions[,1] <- NULL     #matrix1 <- subset(matrix1, select=-c(X)), dit is het voor een matrix
interactions1=as.matrix(interactions)

##---------------------------------
## Make a foodweb for each sample

fwmetrics <- NULL

# pdf("~...empty folder directory for plot pdfs")
for(s in unique(speciesobserved$sampleyear)){
  for (i in unique(speciesobserved$fieldlabel)){
    for (c in c("Benthic", "All")){
      # subset the matrix per cage
      data <- speciesobserved%>%
        filter(sampleyear==s & fieldlabel==i)
      
      if(nrow(data)==0) {
        print(paste("skipped", s,i,c,sep=",")); next}
      
      Macroalgae <- data.frame(sampleyear=s, fieldlabel= i, species_name= "Macroalgae", abundance = 1, mobile= 1, benthic= 1)
      Phytoplankton <- data.frame(sampleyear=s, fieldlabel= i, species_name= "Phytoplankton", abundance = 1, mobile= 1, benthic= 1)
      Zooplankton <- data.frame(sampleyear=s, fieldlabel= i, species_name= "Zooplankton", abundance = 1, mobile= 1, benthic= 1)
      wPOM <- data.frame(sampleyear=s, fieldlabel= i, species_name= "wPOM", abundance = 1, mobile= 1, benthic= 1)
      sPOM <- data.frame(sampleyear=s, fieldlabel= i, species_name= "sPOM", abundance = 1, mobile= 1, benthic= 1)
      
      data <- rbind(data, Macroalgae, Phytoplankton, Zooplankton, wPOM, sPOM)
      
      if(c=="Benthic"){species=unique(data[data$benthic==1,]$species_name) } else{
          species=unique(data$species_name)}
      
      interact<- interactions%>%
        filter(row.names(.) %in% species)%>%
        dplyr::select(all_of(species))%>%
        relocate(any_of(names(interactions)))
     
      matrix=as.matrix(interact)
      mode(matrix) = "numeric"
      
      network <- graph.adjacency(matrix, mode="undirected")
      
      #all metrics
      metrics <- Foodwebmetrics(s,i,c, interact, network)
      
      # Using rbind() to append the output of one iteration to the dataframe
      fwmetrics = rbind(fwmetrics, metrics)

      #calculate Trophic level
      TL <- TrophInd(matrix) #trophic level (with cannibals included)
      TL$round <- as.factor(findInterval(TL$TL, c(0, 1.9, 2.2, 2.5, 3.2, 4), all.inside=T))
      V(network)$TL<- TL$TL[match(V(network)$name, row.names(TL))]
      V(network)$round <- TL$round[match(V(network)$name, row.names(TL))]
      
      coordsp <- cbind(runif(length(V(network))), TL$TL)
      pal=viridis(5)
# 
#       figurepath <- file.path("~/figure folder", paste("Foodweb", s, i, c, ".pdf", sep = ""))
# 
#       pdf(file=figurepath)
# 
#       #Save plot as pdf
#       print(
#         ggplot(ggnetwork(network, layout=coordsp, scale = FALSE, cell.jitter = 0.25), aes(x, y, xend = xend, yend = yend))+
#         geom_edges(size = 0.3, color = "black", arrow = arrow(length = unit(6,"pt"), type = "closed"), curvature = 0.1) +
#         geom_nodes(aes(color = round), size=8)+
#         # geom_nodelabel_repel(aes(label=name, color= sg), size=2, fontface = "bold", box.padding = unit(0.1, "lines"))+
#         # scale_colour_viridis()+
#         scale_colour_manual(values=pal)+
#         theme_classic()+
#         ylim(c(0.9,3.9))+
#         # guides(colour = guide_colourbar(order = 1,  barwidth = 0.75, nbin = 6, ticks.colour='black', frame.colour='black'),
#         # alpha = guide_legend(order = 2), size = guide_legend(order = 3))+
#         coord_cartesian(clip = "off") +
#         theme(legend.position="none", text = element_text(size = 20),axis.line.x = element_blank(),
#               axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank())+
#           labs(y="Trophic level")
#         # labs(title = paste("Foodweb", s, i, c, sep = " "))
#         )
#       dev.off()

      print(paste(s,i,c,sep=","))
      
    }
  }
}

dev.off()

##---------------------------------
## Reorganise metrics

metrics <- as.data.frame(fwmetrics)
  
metrics <- metrics %>% 
  mutate_at(c(4:21), as.numeric)%>%
  mutate(Substrate=factor(substr(metrics$Field.ID, 1, 1)),
         Block=factor(substr(metrics$Field.ID, 2, 3)))

metrics$Part_Wadden_Sea <- ifelse(metrics$Block==5 | metrics$Block==6 |metrics$Block==7,
                                  "East", "West")

metrics$substrate <- factor(ifelse(metrics$Substrate=="B", "BESE", 
                                   ifelse(metrics$Substrate=="C", "Granite", 
                                          ifelse(metrics$Substrate=="P", "Pebbles",
                                                 ifelse(metrics$Substrate=="R", "Reefs",
                                                        ifelse(metrics$Substrate=="S", "Shells",
                                                               ifelse(metrics$Substrate=="W" , "Wood",
                                                                      ifelse(metrics$Substrate=="X", "Sand", ""))))))))

#Filter out all for sand (should not be analysed)
metrics <- metrics %>%
  filter(!grepl('X', Substrate) | !Community=="All") %>%
  mutate_if(is.character,as.factor) %>%
  mutate(Year=factor(Year))%>%
  mutate(substrate=fct_relevel(substrate,c("BESE","Granite","Pebbles", "Reefs", "Shells", "Wood", "Sand")))

# Benthic only
benthic <- metrics %>%
  filter(Community=="Benthic")

# All, including mobile species
all <- metrics %>%
  filter(Community=="All")

##---------------------------------
## Analyse the food web; Link density and connectance

## Link density

# benthic species
lnLDb <- lmer(data=benthic, log1p(Linkdensity)~Substrate+Year+(1|Block))
lnLDb1 <- lm(data=benthic, log1p(Linkdensity)~Substrate+Year)

AIC(lnLDb)-AIC(lnLDb1)
Anova(lnLDb)
# summary(lnLDb)
meansld <- emmeans(object=lnLDb, spec="Substrate")

pairs(meansld, adjust='fdr')
pairs(meansld, adjust="none")

cldld <- data.frame(cld(object=meansld, adjust = "none", Letters = LETTERS, alpha=0.05))
labelld <- benthic %>%
  group_by(Substrate, substrate)%>%
  summarize()%>%
  ungroup()%>%
  right_join(., cldld) %>%
  rename(group = .group)

C <-ggplot(benthic, aes(x=substrate, y=log1p(Linkdensity), fill=substrate))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(aes(pch=Year), position=position_jitter(0.1), size=3)+
  geom_text(data = labelld, aes(label = group, y = 0.75), 
              position = position_dodge(width = .5), show.legend = FALSE, size=7 ) +
  scale_shape_manual(values=c(1, 16), name = "Year")+
  scale_fill_manual(values=c("#01665E","#01665E","#01665E","#01665E","#01665E","#01665E","#FFAD65"), guide="none")+
  labs(y=bquote(log(1 + "Link Density")))+
  theme_classic()+
  theme(text = element_text(size = 20))

# all species
lnLDa <- lmer(data=all, log1p(Linkdensity)~Substrate+Year+(1|Block))
lnLDa1 <- lm(data=all, log1p(Linkdensity)~Substrate+Year)
AIC(lnLDa)-AIC(lnLDa1)  
Anova(lnLDa)

labelLDa <- data.frame(substrate=c("BESE","Granite","Pebbles", "Reefs", "Shells", "Wood"),
                               group=c("A","A","A","A","A","A"))
  
Ca <- ggplot(all, aes(x=substrate, y=log1p(Linkdensity)))+
  geom_boxplot(outlier.shape=NA, fill="lightgrey")+
  geom_point(aes(pch=Year), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16))+
  geom_text(data = labelLDa, aes(label = group, y = 0.75), 
            position = position_dodge(width = .5), show.legend = FALSE, size=7)+
  labs(y=bquote(log(1 + "Link Density")))+
  theme_classic()+
  theme(text = element_text(size = 20))

## Connectance
# benthic species
COb <- lmer(data=benthic, Connectance~Substrate+Year+(1|Block))
COb1 <- lm(data=benthic, Connectance~Substrate+Year)
AIC(COb)-AIC(COb1)
Anova(COb)

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

D <- ggplot(benthic, aes(x=substrate, y=Connectance, fill=substrate))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(aes(pch=Year), position=position_jitter(0.1), size=3)+
  geom_text(data = labelCO, aes(label = group, y = 0.075), 
            position = position_dodge(width = .5), show.legend = FALSE,  size=7) +
  scale_shape_manual(values=c(1, 16), name = "Year")+
  scale_fill_manual(values=c("#01665E","#01665E","#01665E","#01665E","#01665E","#01665E","#FFAD65"), guide="none")+
  theme_classic()+
  theme(text = element_text(size = 20))

#all species
COa <- lmer(data=all, Connectance~Substrate+Year+(1|Block))
COa1 <- lm(data=all, Connectance~Substrate+Year)
AIC(COa)-AIC(COa1)
Anova(COa)

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

Da <- ggplot(all, aes(x=substrate, y=Connectance))+
  geom_boxplot(outlier.shape=NA, fill="lightgrey")+
  geom_point(aes(pch=Year), position=position_jitter(0.1), size=3)+
  scale_shape_manual(values=c(1, 16))+
  geom_text(data = labelCOa, aes(label = group, y = 0.09), 
            position = position_dodge(width = .5), show.legend = FALSE, size=7) +
  theme_classic()+
  theme(text = element_text(size = 20)) 

ggarrange(C,D,labels = c("C", "D"), font.label = list(size = 30), ncol = 2, nrow = 1, common.legend = T, legend = "right")
ggarrange(Ca,Da,labels = c("C", "D"), font.label = list(size = 30), ncol = 2, nrow = 1, common.legend = T, legend = "right")

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

year <- benthic %>%
  group_by(Year) %>%
  summarise(LD=mean(Linkdensity, na.rm=T), CO=mean(Connectance, na.rm=T))

(year[year$Year==2021,]$CO-year[year$Year==2020,]$CO)/
  year[year$Year==2020,]$CO*100

substrate <- benthic %>%
  group_by(Substrate) %>%
  summarise(LD=mean(Linkdensity, na.rm=T), CO=mean(Connectance, na.rm=T),
            median=median(Linkdensity, na.rm=T))
  
(mean(substrate[substrate$Substrate=="W" | substrate$Substrate=="S" | substrate$Substrate=="C",]$LD)- 
    substrate[substrate$Substrate=="X",]$LD)/substrate[substrate$Substrate=="X",]$LD*100
  
(substrate[substrate$Substrate=="B",]$median-substrate[substrate$Substrate=="X",]$median)/substrate[substrate$Substrate=="X",]$median*100