---
title: "SEM_seagrass_data_formating"
author: "Jade Millot"
date: "2023-08-28"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE, eval = T, cache = TRUE)
```

```{r, echo = FALSE}
rm(list=ls())
graphics.off()
```

# Load library

```{r packages_analysis, echo = TRUE}

# Data Manipulation
library(dplyr)
library(tidyr)
library(plyr)
library(reshape2)
library(data.table)
library(stringr)
library(lubridate)
library(vegan)
library(BiodiversityR)
library(tidyverse)
library(magrittr)
library(ade4)     
library(factoextra)
library(klaR)
library(diverse)
library(piecewiseSEM)
library(MASS)
library(nlme)
library(lme4)
library(janitor)
library(DHARMa)
library(benthos)
library(VIM)

#Missing data
library(VIM)
library(naniar)

#Visualisation
library(ggplot2)
library(ggcorrplot)
library(kableExtra)
library(viridis)
library(visdat)
library(naniar)
library(corrplot)
library(ggrepel)
library(gridExtra)
library(ggforce)
library(ggpubr)
library(ggforce)
```

# Biometry data 

## Import data

```{r import_biometrie}
# Biometry before 2016
load("Input/Extract_marben_biom_bf16_30_octobre_2020.Rdata")

# Biometry after 2016
load("Input/Extract_marben_biom_af16_30_octobre_2020.Rdata")

# Biometry between 2019-2020
biom19.20 <- read.csv("Input/Biometrie_HI_2019_2020.csv", sep = ";",
    dec = ".")
```

## Data formating

```{r}
marben_biom_after_2016 <- marben_biom_after_2016 %>% mutate_if(is.character, as.factor)
summary(marben_biom_after_2016)

biom19.20 <- biom19.20 %>% mutate_if(is.character, as.factor)
summary(biom19.20)

marben_biom_before_2016 <- marben_biom_before_2016 %>% mutate_if(is.character, as.factor)
summary(marben_biom_before_2016)

# biom19.20 formating
#------------------------

biom19.20.cl <- biom19.20 %>%
  dplyr::rename(date = DATE, point = POINT, collection = PRELEVEMENT, pied = ECHANTILLON, value = RESULTAT, number = NUMERO_INDIVIDU, variable = PARAMETRE, value = RESULTAT, units = UNITE.1) %>%
  mutate(broken = case_when(
    RESULTAT_QUALITATIF_CODE_SANDRE == "Cassee" ~ TRUE, 
    TRUE ~ FALSE
  )) %>%  
  dplyr:: select(date, point, collection, pied, value, number, broken, variable, units)

# Dates formating
#------------------------

# Dates formating to remove hour
marben_biom_after_2016$date <- as.Date(marben_biom_after_2016$date,format="%d-%m-%y")
marben_biom_before_2016$date <- as.Date(marben_biom_before_2016$date,format="%d-%m-%y")
biom19.20.cl$date <- dmy(biom19.20.cl$date)

# Creation of a variable "year"
marben_biom_after_2016$annee <- as.factor(format.Date(marben_biom_after_2016$date,"%Y"))
marben_biom_before_2016$annee <- as.factor(format.Date(marben_biom_before_2016$date,"%Y"))
biom19.20.cl$annee <- as.factor(format.Date(biom19.20.cl$date,"%Y"))

# Creation of a variable "month"
marben_biom_after_2016$mois <- as.numeric(format.Date(marben_biom_after_2016$date,"%m"))
marben_biom_before_2016$mois <- as.numeric(format.Date(marben_biom_before_2016$date,"%m"))
biom19.20.cl$mois <- as.numeric(format.Date(biom19.20.cl$date,"%m"))

# Creation of a variable "season"
marben_biom_after_2016$saison[marben_biom_after_2016$mois<7] <- "Printemps"
marben_biom_after_2016$saison[marben_biom_after_2016$mois>=7] <- "Automne"
marben_biom_after_2016$saison<- as.factor(marben_biom_after_2016$saison)

marben_biom_before_2016$saison[marben_biom_before_2016$mois<7] <- "Printemps"
marben_biom_before_2016$saison[marben_biom_before_2016$mois>=7] <- "Automne"
marben_biom_before_2016$saison<- as.factor(marben_biom_before_2016$saison)

biom19.20.cl$saison[biom19.20.cl$mois<7] <- "Printemps"
biom19.20.cl$saison[biom19.20.cl$mois>=7] <- "Automne"
biom19.20.cl$saison<- as.factor(biom19.20.cl$saison)

# Site names formating
#---------------------

# Habitat mention removed
site.split1 <- str_split(marben_biom_after_2016$site," ")
for (i in 1:length(site.split1)){
  site.split1[[i]]<- site.split1[[i]][1:(length(site.split1[[i]])-1)] 
}
site1 <- sapply(X = site.split1,paste,collapse=" ")

# Special features changed
site1 <- str_replace_all(site1,"L'","")
site1 <- str_replace_all(site1,"â","a")
site1 <- str_replace_all(site1,"â","a")
site1 <- str_replace_all(site1,"é","e")
site1 <- str_replace_all(site1,"è","e")
site1 <- str_replace_all(site1,"ë","e")
site1 <- str_replace_all(site1,"î","i")

# # Habitat mention removed
site.split2 <- str_split(marben_biom_before_2016$site," ")
for (i in 1:length(site.split2)){
  site.split2[[i]]<- site.split2[[i]][1:(length(site.split2[[i]])-1)] 
}
site2 <- sapply(X = site.split2,paste,collapse=" ")

# Special features changed
site2 <- str_replace_all(site2,"L'","")
site2 <- str_replace_all(site2,"â","a")
site2 <- str_replace_all(site2,"â","a")
site2 <- str_replace_all(site2,"é","e")
site2 <- str_replace_all(site2,"è","e")
site2 <- str_replace_all(site2,"ë","e")
site2 <- str_replace_all(site2,"î","i")

# Biom19.20
# Site name extracted
site.split3 <- str_split(biom19.20.cl$point," ")
for (i in 1:length(site.split3)){
  site.split3[[i]] <- site.split3[[i]][1]
}
site3 <- sapply(X = site.split3,paste,collapse=" ")

# Special features changed
site3 <- str_replace_all(site3,"L'","")
site3 <- str_replace_all(site3,"l'","")
site3 <- str_replace_all(site3,"â","a")
site3 <- str_replace_all(site3,"â","a")
site3 <- str_replace_all(site3,"é","e")
site3 <- str_replace_all(site3,"è","e")
site3 <- str_replace_all(site3,"ë","e")
site3 <- str_replace_all(site3,"î","i")

# New site names saved
marben_biom_after_2016$site <- as.factor(site1)
marben_biom_before_2016$site <- as.factor(site2)
biom19.20.cl$site <- as.factor(site3)

# Visualize new site name
unique(marben_biom_after_2016$site)
unique(marben_biom_before_2016$site)
unique(biom19.20.cl$site)

# Point name cleaning
marben_biom_after_2016$point <- str_replace_all(marben_biom_after_2016$point,"Herbier ","") 
marben_biom_after_2016$point <- str_replace_all(marben_biom_after_2016$point,"l'","")

marben_biom_before_2016$point <- str_replace_all(marben_biom_before_2016$point,"Herbier ","")
marben_biom_before_2016$point <- str_replace_all(marben_biom_before_2016$point,"H","")
marben_biom_before_2016$point <- str_replace_all(marben_biom_before_2016$point,"l'","")

biom19.20.cl$point <- str_replace_all(biom19.20.cl$point,"Herbier ","")
biom19.20.cl$point <- str_replace_all(biom19.20.cl$point,"H","")
biom19.20.cl$point <- str_replace_all(biom19.20.cl$point,"l'","")

# Visualize point name
unique(marben_biom_after_2016$point) 
unique(marben_biom_before_2016$point)
unique(biom19.20.cl$point)

# Summary
marben_biom_after_2016 <- marben_biom_after_2016 %>% mutate_if(is.character, as.factor)
summary(marben_biom_after_2016)

marben_biom_before_2016 <- marben_biom_before_2016 %>% mutate_if(is.character, as.factor)
summary(marben_biom_before_2016)

biom19.20.cl <- biom19.20.cl %>% mutate_if(is.character, as.factor)
summary(biom19.20.cl)

# Data selection
#---------------

# Type of collection
unique(marben_biom_after_2016$collection)
unique(marben_biom_before_2016$collection)
unique(biom19.20.cl$collection)

# "Haveneau", "Carotte" and "Granulo" collections are removed
marben_biom_before_2016 <- marben_biom_before_2016 %>% 
  filter (!collection %in% c("Haveneau 1", "Haveneau 2", "Haveneau 3", "Carotte 1", "Carotte 2", "Carotte 3", "Granulo")) %>% 
  droplevels()

# Visualize variables
unique(marben_biom_after_2016$variable)
unique(marben_biom_before_2016$variable)
unique(biom19.20.cl$variable)

# Column "theme" added
marben_biom_after_2016 <- marben_biom_after_2016 %>% 
  mutate (theme = "Herbiers Intertidaux") %>% 
  mutate_if(is.character, as.factor)

marben_biom_before_2016 <- marben_biom_before_2016 %>% 
  mutate (theme = "Herbiers Intertidaux") %>% 
  mutate_if(is.character, as.factor)

biom19.20.cl <- biom19.20.cl %>% 
  mutate (theme = "Herbiers Intertidaux") %>% 
  mutate_if(is.character, as.factor)

# Variables selection
marben_biom_after_2016 <- marben_biom_after_2016 %>%
  dplyr::select(theme, site, point_name = point, date, annee, saison, collection, pied, feuille = number, variable, value, units, broken)

marben_biom_before_2016 <- marben_biom_before_2016 %>% 
  dplyr::select(theme, site, point_name = point, date, annee, saison, collection, pied = number, variable, value, units, broken)

biom19.20.cl <- biom19.20.cl %>% 
  dplyr::select(theme, site, point_name = point, date, annee, saison, collection, pied, feuille = number, variable, value, units, broken)

summary(marben_biom_after_2016)
summary(marben_biom_before_2016)

biom19.20.cl <-  biom19.20.cl %>%
  mutate(units = str_replace_all(.$units, "243", "%"), value = str_replace_all(.$value, "mm.0", "0")) %>% 
  mutate(value = as.numeric(value)) %>%
  mutate_if(is.character, as.factor)

summary(biom19.20.cl)

# Which variables measured in one quadrat ?
marben_biom_before_2016 %>%
  filter(collection == "Quadrat 1") %>% 
  distinct(collection, variable)

# Which variables measured in 1/2 quadrat ?
marben_biom_before_2016 %>%
  filter(collection == "1/2 quadrat 1") %>% 
  distinct(collection, variable)

# Which variables measured in 1/2 quadrat ?
biom19.20.cl %>%
  filter(collection == "1/2 quadrat 1") %>% 
  distinct(collection, variable)
```


## Density data from 1/2 quadrat

Aurélien :
"Normally, a sheath length and leaf width measurement is made for each foot of the half-quadrats, and the number of these measurements is counted to obtain the number of feet.
If one of the two measurements has not been made, we take the maximum number of measurements that have been made between one and the other of the variables".

### Before 2016

```{r}
# Data selection 
#---------------
demi_bf_16 <-  marben_biom_before_2016 %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  filter(variable %in% c("Largeur moyenne des feuilles de Zostera marina", "Longueur de la gaine du pied de Zostera marina")) %>% 
  filter (annee %in% (2004:2015)) %>%
  droplevels() 

# Check duplicated lines
dup_demi_bf_16 <-  demi_bf_16 %>% dplyr::select(-value, -units)
dup_demi_bf_16 <- dup_demi_bf_16[duplicated(dup_demi_bf_16),] 
dup_demi_bf_16 <- left_join(dup_demi_bf_16, demi_bf_16) 
dup_demi_bf_16 

# Marion 09/12/20 : "The max values for duplicated lines (for each site) are actually leaves lengths
# Correction :
correction <- dup_demi_bf_16
correction$variable <- as.character(correction$variable)
correction$variable[correction$value == 46] <- "Longueur de la feuille de Zostera marina"
correction$variable[correction$value == 398] <- "Longueur de la feuille de Zostera marina"
correction$variable[correction$value == 133] <- "Longueur de la feuille de Zostera marina"
correction$variable[correction$value == 300] <- "Longueur de la feuille de Zostera marina"
correction$variable <- as.factor(correction$variable)
correction

# Corrections in marben_biom_before_2016
marben_biom_before_2016 <- marben_biom_before_2016 %>% 
  anti_join(dup_demi_bf_16) %>%
  rbind(correction)

# Data are re selected
demi_bf_16 <-  marben_biom_before_2016 %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  filter(variable %in% c("Largeur moyenne des feuilles de Zostera marina", "Longueur de la gaine du pied de Zostera marina")) %>% 
  filter (annee %in% (2004:2015)) %>%
  droplevels() %>%
  spread(variable, value)

#–> Number of shoot in 1/2 quadrat = density
# Count of measure per 1/2 quadrat 
#---------------------------------
dens_bf_16 <- demi_bf_16 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(n = n()) %>%
  ungroup

# Convert in  m2 
#----------------
dens_bf_16 <- dens_bf_16 %>% 
  mutate(n = n/0.05) %>%
  dplyr::rename(density = n)

summary(dens_bf_16)

# Mean of both 1/2 quadrats
#--------------------------
dens.moy_bf_16 <- dens_bf_16 %>%
  group_by(theme, site, point_name, date, annee, saison) %>%
  summarise_at(vars(density), mean, na.rm= TRUE) %>%
  ungroup()
```

### After 2016
2016-2017 : “1/2 quadrat 1” and “1/2 quadrat 2” are substituted for “Biométrie 1” “Biométrie 2” which correspond to quadrats of 0.1m2 (there are not 1/2 quadrat)
After 2017 : normal measures => 1/2 quadrat = 0.05 m2

```{r}
# Data selection 
demi_16_17 <- marben_biom_after_2016 %>%
  filter(collection %in% c("Biométrie 1", "Biométrie 2")) %>%
  filter(variable %in% c("Largeur moyenne des feuilles de Zostera marina", "Longueur de la gaine du pied de Zostera marina")) %>% 
  droplevels()

demi_sup_17 <- marben_biom_after_2016 %>%
  # 2019-2020 are removed
  filter(annee != 2019 & annee != 2020) %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  filter(variable %in% c("Largeur moyenne des feuilles de Zostera marina", "Longueur de la gaine du pied de Zostera marina")) %>% 
  droplevels()

demi_19_20 <- biom19.20.cl %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  filter(variable %in% c("Largeur moyenne des feuilles de Zostera marina", "Longueur de la gaine du pied de Zostera marina")) %>% 
  droplevels()

# Check for duplicated lines
dup_demi_16_17 <-  demi_16_17 %>% dplyr::select(-value, -units)
dup_demi_16_17 <- dup_demi_16_17[duplicated(dup_demi_16_17),]
dup_demi_16_17 <- left_join(dup_demi_16_17, demi_16_17)
dup_demi_16_17 %>% 
  kbl(.) %>% 
  kable_paper %>% 
  scroll_box(width = "900px")

dup_demi_sup_17 <-  demi_sup_17 %>% dplyr::select(-value, -units)
dup_demi_sup_17 <- dup_demi_sup_17[duplicated(dup_demi_sup_17),]
dup_demi_sup_17 <- left_join(dup_demi_sup_17, demi_sup_17)
dup_demi_sup_17 %>% 
  kbl(.) %>% 
  kable_paper %>% 
  scroll_box(width = "900px")

dup_demi_19_20 <-  demi_19_20 %>% dplyr::select(-value, -units)
dup_demi_19_20 <- dup_demi_19_20[duplicated(dup_demi_19_20),]
dup_demi_19_20 <- left_join(dup_demi_19_20, demi_19_20)
dup_demi_19_20 %>% 
  kbl(.) %>% 
  kable_paper %>% 
  scroll_box(width = "900px")

# Correction
# Variable name changed with "Longueur de la feuille de Zostera marina"
correction <- dup_demi_16_17
correction$variable <- as.character(correction$variable)
correction$variable[correction$value == 100] <- "Longueur de la feuille de Zostera marina"
correction$variable[correction$value == 203] <- "Longueur de la feuille de Zostera marina"
correction$variable[correction$value == 363] <- "Longueur de la feuille de Zostera marina"
correction$feuille[correction$value == 75] <- 5
correction$feuille[correction$value == 4] <- 4
correction$variable <- as.factor(correction$variable)
correction

# Correction in marben_biom_after_2016
marben_biom_after_2016 <- marben_biom_after_2016 %>% 
  anti_join(dup_demi_16_17) %>% 
  rbind(correction)

# Variable name changed with "Longueur de la feuille de Zostera marina" in demi_sup_17 
correction <- dup_demi_sup_17
correction$variable <- as.character(correction$variable)
correction$variable[correction$value == 4.5] <- "Largeur moyenne des feuilles de Zostera marina"
correction$variable[correction$value == 390] <- "Longueur de la feuille de Zostera marina"
correction$variable[correction$value == 127] <- "Longueur de la feuille de Zostera marina"
correction$variable <- as.factor(correction$variable)
correction

# Correction in marben_biom_after_2016
marben_biom_after_2016 <- marben_biom_after_2016 %>% 
  anti_join(dup_demi_sup_17) %>% 
  rbind(correction)

# Data re selected
demi_16_17 <- marben_biom_after_2016 %>%
  filter(collection %in% c("Biométrie 1", "Biométrie 2")) %>%
  filter(variable %in% c("Largeur moyenne des feuilles de Zostera marina", "Longueur de la gaine du pied de Zostera marina")) %>% 
  droplevels() %>% 
  spread(variable, value)

demi_sup_17 <- marben_biom_after_2016 %>%
  # 2019-2020 are removed
  filter(annee != 2019 & annee != 2020) %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  filter(variable %in% c("Largeur moyenne des feuilles de Zostera marina", "Longueur de la gaine du pied de Zostera marina")) %>% 
  droplevels() %>%
  spread(variable, value)

demi_19_20 <- demi_19_20 %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  filter(variable %in% c("Largeur moyenne des feuilles de Zostera marina", "Longueur de la gaine du pied de Zostera marina")) %>% 
  droplevels() %>%
  spread(variable, value)

# Count of measures per 1/2 quadrat
#--------------------------------
dens_16_17 <- demi_16_17 %>% 
  distinct(theme, site, point_name, date, annee, saison, collection, pied) %>%
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(n = n()) %>%
  ungroup

dens_sup_17 <- demi_sup_17 %>%
  distinct(theme, site, point_name, date, annee, saison, collection, pied) %>%
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(n = n()) %>%
  ungroup

dens_19_20 <- demi_19_20 %>%
  distinct(theme, site, point_name, date, annee, saison, collection, pied) %>%
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(n = n()) %>%
  ungroup

# Convert in  m2 
#----------------
dens_16_17 <- dens_16_17 %>%
  mutate(n = n/0.1) %>% 
  dplyr::rename(density = n)

summary(dens_16_17)

dens_sup_17 <- dens_sup_17 %>%
  mutate(n = n/0.05) %>% 
  dplyr::rename(density = n)

summary(dens_sup_17)

dens_19_20 <- dens_19_20 %>%
  mutate(n = n/0.05) %>%
  dplyr::rename(density = n)

summary(dens_19_20)

# Mean of both 1/2 quadrats
#--------------------------
dens.moy_16_17 <- dens_16_17 %>%
  group_by(theme, site, point_name, date, annee, saison) %>%
  summarise_at(vars(density), mean, na.rm= TRUE) %>%
  ungroup()

dens.moy_sup_17 <- dens_sup_17 %>%
  group_by(theme, site, point_name, date, annee, saison) %>%
  summarise_at(vars(density), mean, na.rm= TRUE) %>%
  ungroup()

dens.moy_19_20 <- dens_19_20 %>%
  group_by(theme, site, point_name, date, annee, saison) %>%
  summarise_at(vars(density), mean, na.rm= TRUE) %>%
  ungroup()
```

### Binding

```{r}
density.tot <- rbind(dens.moy_bf_16, dens.moy_16_17, dens.moy_sup_17, dens.moy_19_20)
summary(density.tot)

# Plot of densities
ggplot(data = density.tot, aes(x= date, y = density)) + 
  geom_point(aes(color = annee)) + 
  facet_wrap(~ site)

# Check outlayers
demi_bf_16 %>% filter(site == "Arcouest" & annee == 2008) %>% distinct (point_name, annee, pied)
```

## Leaves biomass per m2 (from 1/2 quadrats)
Biomass after 2016 is found in the extraction before 2016 because the path in the database has not changed for this variable.

```{r}

# Data selection
#---------------

bf <- marben_biom_before_2016 %>% 
  # 2019 and 2020 are removed 
  filter(annee != 2019 & annee != 2020) %>%
  filter(variable == "Biomasse des feuilles de Zostera marina (Herbier REBENT)") %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  droplevels()

bf_16_17 <- marben_biom_before_2016 %>% 
  filter(variable == "Biomasse des feuilles de Zostera marina (Herbier REBENT)") %>%
  filter(collection %in% c("Biométrie 1", "Biométrie 2")) %>%
  droplevels()

bf_19_20 <- biom19.20.cl %>% 
  filter(variable == "Biomasse des feuilles de Zostera marina (Herbier REBENT)") %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  droplevels()

# Convertion in m2
#-----------------
bf.m2 <- bf %>% 
  mutate (bf.m2 = value / 0.05) %>% 
  dplyr:: select(-value, -units, -pied, -variable)

bf_16_17.m2 <- bf_16_17 %>% 
  mutate (bf.m2 = value / 0.1) %>% 
  dplyr:: select(-value, -units, -pied, -variable)

bf_19_20.m2 <- bf_19_20 %>% 
  mutate (bf.m2 = value / 0.05) %>% 
  dplyr:: select(-value, -units, -pied, -variable, -feuille)

# Check for outlayers 
#--------------------
bf %>% filter(site == "Arcouest" & annee == 2009)
bf %>% filter(site == "Glenan" & annee == 2009)
```

## Root biomass per m2 (from 1/2 quadrats)

```{r}
# Data selection
#---------------
br <- marben_biom_before_2016 %>% 
  # 2019 and 2020 are removed
  filter(annee != 2019 & annee != 2020) %>%
  filter(variable == "Biomasse des racines de Zostera marina") %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  droplevels()

br_16_17 <- marben_biom_before_2016 %>% 
  filter(variable == "Biomasse des racines de Zostera marina") %>%
  filter(collection %in% c("Biométrie 1", "Biométrie 2")) %>%
  droplevels()

br_19_20 <- biom19.20.cl %>% 
  filter(variable == "Biomasse des racines de Zostera marina") %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  droplevels()

# Conversion in m2
#-----------------
br.m2 <- br %>% 
  mutate (br.m2 = value / 0.05) %>% 
  dplyr:: select(-value, -units, -pied, -variable)

br_16_17.m2 <- br_16_17 %>% 
  mutate (br.m2 = value / 0.1) %>% 
  dplyr:: select(-value, -units, -pied, -variable)

br_19_20.m2 <- br_19_20 %>% 
  mutate (br.m2 = value / 0.05) %>% 
  dplyr:: select(-value, -units, -pied, -variable, -feuille)

# Check for outlayers
#---------------------
br %>% filter(site == "Callot" & annee == 2018)
br %>% filter(site == "Molene" & annee  %in% c(2017, 2019))
br %>% filter(site == "Roscanvel" & annee %in% c(2012, 2018))
```

## Biomasses binded

```{r}
# Binding
#----------------
bf_tot <- rbind(bf.m2, bf_16_17.m2, bf_19_20.m2)
br_tot <- rbind(br.m2, br_16_17.m2, br_19_20.m2)

biomasse <- full_join(bf_tot, br_tot)

# Mean on both demi-quadrats
#---------------------------
biomasse <- biomasse %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(bf.m2, br.m2), mean, na.rm= TRUE) %>% 
  ungroup()

# Plot of biomasses 
#------------------

ggplot(data = biomasse, aes(x = date, y = bf.m2)) + 
  geom_point(aes(color = annee)) +
  facet_wrap(~ site)

ggplot(data = biomasse, aes(x = date, y = br.m2)) + 
  geom_point(aes(color = annee)) +
  facet_wrap(~ site)
```

## Average sheathes length (from) 1/2 quadrats)

### Before 2016

```{r}
# Data selection
#---------------
long.g_bf_16 <- marben_biom_before_2016 %>% 
  filter(annee %in% c(2004:2015)) %>%
  filter (variable == "Longueur de la gaine du pied de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  # Remove broken sheathes
  filter (broken != 1)

# Average length per point
#----------------------
long.g.moy_bf16 <- long.g_bf_16 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(long.g.moy = value) %>%
  ungroup()
```

### After 2016

```{r}
# Data selection
#---------------
long.g_16_17 <- marben_biom_after_2016 %>% 
  filter (variable == "Longueur de la gaine du pied de Zostera marina") %>%
  filter (collection %in% c("Biométrie 1", "Biométrie 2")) %>%
  # Remove broken sheathes
  filter (broken != 1)

long.g_sup_17 <- marben_biom_after_2016 %>% 
  # 2019 and 2020 removed 
  filter(annee != 2019 & annee != 2020) %>%
  filter (variable == "Longueur de la gaine du pied de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  # Remove broken sheathes
  filter (broken != 1)

long.g_19_20 <- biom19.20.cl %>% 
  filter (variable == "Longueur de la gaine du pied de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  # Remove broken sheathes
  filter (broken != 1)

# Average length per point
#-------------------------
long.g.moy_16_17 <- long.g_16_17 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(long.g.moy = value) %>%
  ungroup()

long.g.moy_sup_17 <- long.g_sup_17 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(long.g.moy = value) %>%
  ungroup()

long.g.moy_19_20 <- long.g_19_20 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(long.g.moy = value) %>%
  ungroup()
```

### Binding

```{r}
long.g.moy <- rbind(long.g.moy_bf16, long.g.moy_16_17, long.g.moy_sup_17, ... = long.g.moy_19_20)
  summary(long.g.moy)

# Plot of sheath lenght
#----------------------
ggplot(data = long.g.moy, aes(x= date, y = long.g.moy)) + 
  geom_point(aes(color = annee)) + 
  facet_wrap(~ site)
```

## Average leaves length

Aurélien :
"Here we take into account the broken leaves in the average because represents real complexity of the habitat and the mechanical effect of the environment + physio"

### Before 2016

```{r}
# Data selection
#---------------
long.f_bf_16 <- marben_biom_before_2016 %>% 
  filter(annee %in% c(2004:2015)) %>%
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2"))

# Average length per point
#-------------------------
long.f.moy_bf16 <- long.f_bf_16 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(long.f.moy = value) %>%
  ungroup()
```

### After 2016

```{r}
# Data selection
#---------------
long.f_16_17 <- marben_biom_after_2016 %>% 
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("Biométrie 1", "Biométrie 2")) 

long.f_sup_17 <- marben_biom_after_2016 %>% 
  # 2019 and 2020 removed
  filter(annee != 2019 & annee != 2020) %>%
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) 

long.f_19_20 <- biom19.20.cl %>% 
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) 

# Average length per point
#-------------------------
long.f.moy_16_17 <- long.f_16_17 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(long.f.moy = value) %>%
  ungroup()

long.f.moy_sup_17 <- long.f_sup_17 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(long.f.moy = value) %>%
  ungroup()

long.f.moy_19_20 <- long.f_19_20 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(long.f.moy = value) %>%
  ungroup()
```

### Binding

```{r}
long.f.moy <- rbind(long.f.moy_bf16, long.f.moy_16_17, long.f.moy_sup_17, long.f.moy_19_20)
summary(long.f.moy)

# Plot of leaves length
#----------------------
ggplot(data = long.f.moy, aes(x= date, y = long.f.moy)) + 
  geom_point(aes(color = annee)) + 
  facet_wrap(~ site)
```

## Average leaves width

### Before 2016

```{r}
# Data selection
#---------------
larg.f_bf_16 <- marben_biom_before_2016 %>% 
  filter(annee %in% c(2004:2015)) %>%
  filter (variable == "Largeur moyenne des feuilles de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2"))

# Average width per point
#-------------------------
larg.f.moy_bf16 <- larg.f_bf_16 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(larg.f.moy = value) %>%
  ungroup()
```

### After 2016

```{r}
# Data selection
#---------------
larg.f_16_17 <- marben_biom_after_2016 %>% 
  filter (variable == "Largeur moyenne des feuilles de Zostera marina") %>%
  filter (collection %in% c("Biométrie 1", "Biométrie 2")) 

larg.f_sup_17 <- marben_biom_after_2016 %>% 
  # 2019 and 2020 are removed
  filter(annee != 2019 & annee != 2020) %>%
  filter (variable == "Largeur moyenne des feuilles de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) 

larg.f_19_20 <- biom19.20.cl %>% 
  filter (variable == "Largeur moyenne des feuilles de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) 

# Average width per point
#-------------------------
larg.f.moy_16_17 <- larg.f_16_17 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(larg.f.moy = value) %>%
  ungroup()

larg.f.moy_sup_17 <- larg.f_sup_17 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(larg.f.moy = value) %>%
  ungroup()

larg.f.moy_19_20 <- larg.f_19_20 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(value), mean, na.rm = TRUE) %>% 
  dplyr::rename(larg.f.moy = value) %>%
  ungroup()
```

### Binding

```{r}
larg.f.moy <- rbind(larg.f.moy_bf16, larg.f.moy_16_17, larg.f.moy_sup_17, larg.f.moy_19_20)
summary(larg.f.moy)

# Plot of leaves width
#---------------------
ggplot(data = larg.f.moy, aes(x= date, y = larg.f.moy)) + 
  geom_point(aes(color = annee)) + 
  facet_wrap(~ site)
```

## Number of leaves per shoot

### Before 2016

```{r}
# Data selection
#---------------
n.f_bf_16 <- marben_biom_before_2016 %>% 
  filter(annee %in% c(2004:2015)) %>%
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2"))

# Number of leaves per shoot for each 1/2 quadrat
#------------------------------------------------
n.f_bf_16 <- n.f_bf_16 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection, pied) %>%
  # Number of leaves per shoot (number of lines grouped by "pied")
  dplyr::summarise(n.f = n()) %>% 
  ungroup %>% 
  # Number of leaves per 1/2 quadrat and number of shoot per 1/2 quadrat
  group_by (theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(pied.tot = n(), n.f.tot = sum(n.f)) %>% 
  ungroup %>% 
  # Number of leaves/shoot per 1/2 quadrat
  group_by (theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(n.f.pied = n.f.tot/pied.tot) %>%
  ungroup()

# Mean of bothe 1/2 quadrats
#----------------------------
n.f.pied.moy_bf_16 <- n.f_bf_16  %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(n.f.pied), mean, na.rm = TRUE) %>% 
  ungroup()
```

### After 2016

```{r}
# Data selection
#---------------
n.f_16_17 <- marben_biom_after_2016 %>% 
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("Biométrie 1", "Biométrie 2")) 

n.f_sup_17 <- marben_biom_after_2016 %>% 
  # 2019 and 2020 are removed 
  filter(annee != 2019 & annee != 2020) %>%
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) 

n.f_19_20 <- biom19.20.cl %>% 
  mutate(pied = str_replace_all(.$pied, "I", "i")) %>%
  mutate_if(is.character, as.factor) %>%
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) 

n.f_19_20 %>%
  filter(pied == "")

# /!\ Numerotation of some shoot for Roscanvem_Autumn_2019 are missing, this sample is removed 

n.f_19_20 <- n.f_19_20 %>% 
  filter(pied != "")

# Number of leaves per shoot for each 1/2 quadrat
#------------------------------------------------
n.f_16_17 <- n.f_16_17 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection, pied) %>%
  # Number of leaves per shoot
  dplyr::summarise(n.f = n()) %>%
  ungroup %>% 
  # Number of leaves per 1/2 quadrat and number of shoot per 1/2 quadrat
  group_by (theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(pied.tot = n(), n.f.tot = sum(n.f)) %>%
  ungroup %>% 
  # Number of leaves/shoot per 1/2 quadrat
  group_by (theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(n.f.pied = n.f.tot/pied.tot) %>%
  ungroup()

n.f_sup_17 <- n.f_sup_17 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection, pied) %>%
  # Number of leaves per shoot
  dplyr::summarise(n.f = n()) %>%
  ungroup %>% 
  # Number of leaves per 1/2 quadrat and number of shoot per 1/2 quadrat
  group_by (theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(pied.tot = n(), n.f.tot = sum(n.f)) %>%
  ungroup %>% 
  # Number of leaves/shoot per 1/2 quadrat
  group_by (theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(n.f.pied = n.f.tot/pied.tot) %>%
  ungroup()

n.f_19_20 <- n.f_19_20 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection, pied) %>%
  # Number of leaves per shoot
  dplyr::summarise(n.f = n()) %>%
  ungroup %>% 
  # Number of leaves per 1/2 quadrat and number of shoot per 1/2 quadrat
  group_by (theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(pied.tot = n(), n.f.tot = sum(n.f)) %>%
  ungroup %>% 
  # Number of leaves/shoot per 1/2 quadrat
  group_by (theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(n.f.pied = n.f.tot/pied.tot) %>%
  ungroup()

# Mean of both 1/2 quadrat
#-------------------------
n.f.pied.moy_16_17 <- n.f_16_17  %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(n.f.pied), mean, na.rm = TRUE) %>% 
  ungroup()

n.f.pied.moy_sup_17 <- n.f_sup_17  %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(n.f.pied), mean, na.rm = TRUE) %>% 
  ungroup()

n.f.pied.moy_19_20 <- n.f_19_20  %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(n.f.pied), mean, na.rm = TRUE) %>% 
  ungroup()
```

### Binding

```{r}
n.f.pied.moy <- rbind(n.f.pied.moy_bf_16, n.f.pied.moy_16_17, n.f.pied.moy_sup_17, n.f.pied.moy_19_20)
summary(n.f.pied.moy)

# Plot of leaves number per shoot
ggplot(data = n.f.pied.moy, aes(x= date, y = n.f.pied)) + 
  geom_point(aes(color = annee)) + 
  facet_wrap(~ site)
```

## Percentage of broken leaves

### Before 2016

```{r}
# Data selection
#---------------
brk_bf_16 <- marben_biom_before_2016 %>% 
  filter(annee %in% c(2004:2015)) %>%
  filter (variable == "Longueur de la feuille de Zostera marina") %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2"))

# Percentage of broken leaves per 1/2 quadrat
#--------------------------------------------
brk_bf_16 <- brk_bf_16 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  # Number of leaves per 1/2 quadrat and number of broken leaves per 1/2 quadrat
  dplyr::summarise(n.f = n(), brk = sum(broken)) %>% 
  ungroup() %>%
  # Percentage of broken leaves
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(pourc.brok = (brk/n.f)*100) %>%
  ungroup()

# Mean of both 1/2 quadrats
#--------------------------
brk.moy_bf_16 <- brk_bf_16   %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(pourc.brok), mean, na.rm = TRUE) %>% 
  ungroup()
```

### After 2016

```{r}
# Data selection 
#---------------
brk_16_17 <-  marben_biom_after_2016 %>% 
  filter(variable == "Longueur de la feuille de Zostera marina")%>%   filter(collection %in% c("Biométrie 1", "Biométrie 2")) 

brk_sup_17 <- marben_biom_after_2016 %>% 
  # 2019 and 2020 are removed
  filter(annee != 2019 & annee != 2020) %>%
  filter (variable == "Longueur de la feuille de Zostera marina")%>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) 

brk_19_20 <- biom19.20.cl %>%
  mutate(site=as.character(site))%>%
  filter (variable == c("Longueur de la feuille de Zostera marina","Etat extremite"))%>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2"))

# Percentage of broken leaves per 1/2 quadrat
#--------------------------------------------
brk_16_17 <- brk_16_17 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  # Number of leaves per 1/2 quadrat and number of broken leaves per 1/2 quadrat
  dplyr::summarise(n.f = n(), brk = sum(broken)) %>% 
  ungroup() %>%
  # Percentage of broken leaves
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(pourc.brok = (brk/n.f)*100) %>%
  ungroup()

brk_sup_17 <- brk_sup_17 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  # Number of leaves per 1/2 quadrat and number of broken leaves per 1/2 quadrat
  dplyr::summarise(n.f = n(), brk = sum(broken)) %>% 
  ungroup() %>%
  # Percentage of broken leaves
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(pourc.brok = (brk/n.f)*100) %>%
  ungroup()

brk_19_20 <- brk_19_20 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  # Number of leaves per 1/2 quadrat and number of broken leaves per 1/2 quadrat
  dplyr::summarise(n.f = n(), brk = sum(broken)) %>% 
  ungroup() %>%
  # Percentage of broken leaves 
  group_by(theme, site, point_name, date, annee, saison, collection) %>%
  dplyr::summarise(pourc.brok = (brk/n.f)*100) %>%
  ungroup()

# Mean of both 1/2 quadrats
#--------------------------
brk.moy_16_17 <- brk_16_17   %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(pourc.brok), mean, na.rm = TRUE) %>% 
  ungroup()

brk.moy_sup_17 <- brk_sup_17   %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(pourc.brok), mean, na.rm = TRUE) %>% 
  ungroup()

brk.moy_19_20 <- brk_19_20 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(pourc.brok), mean, na.rm = TRUE) %>% 
  ungroup()
```

### Binding

```{r}
brk.tot <- rbind(brk.moy_bf_16, brk.moy_16_17, brk.moy_sup_17, brk.moy_19_20)
summary(brk.tot)

# Plot of number of broken leaves
#--------------------------------
ggplot(data = brk.tot, aes(x = date, y = pourc.brok)) + 
  geom_point(aes(color = annee)) + 
  facet_wrap(~ site)
```

# Macroalgae and epiphytes data 

Macroalgae and epiphyte data up to 2018 are to be extracted from "marben_biom_before_2016". 
Macroalgae data from 2016 and 2017 are to be extracted from the "Biométrie 1" and "Biométrie 2" collections: area 0.1m2. 
Macroalgae data from 2005 to 2015 and 2018 are to be extracted from the "1/2 quadrat 1" and "1/2 quadrat 1" collections: area 0.5m2. 

## Macroalgae biomass

### Before 2016 + 2018 

```{r}
# Data selection 
#---------------
biom.macro_bf_16 <- marben_biom_before_2016 %>% 
  filter(annee %in% c(2005 : 2015, 2018)) %>%
  filter (variable %in% c("biomasse algues rouges", "biomasse algues vertes", "biomasse algues brunes")) %>%
  filter(collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  dplyr :: select(-pied)
 
# Sum of the three algae types per 1/2 quadrat
#---------------------------------------------
biom.macro_bf_16 <- biom.macro_bf_16 %>% 
  group_by(theme, site, point_name, date, annee, saison, collection) %>% 
  summarise_at(vars(value), sum, na.rm = TRUE) %>% 
  dplyr::rename(biom.macro = value) %>%
  ungroup()

# Conversion in m2
#-----------------
biom.macro_bf_16 <- biom.macro_bf_16 %>%
  mutate(biom.macro = biom.macro/0.5) 

# Average macroalgae biomass per point
#-------------------------------------
biom.macro.moy_bf_16 <- biom.macro_bf_16 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(biom.macro), mean, na.rm = TRUE) %>% 
  dplyr ::rename(biom.macro.moy = biom.macro) %>%
  ungroup()

```

### After 2016 

```{r}
# Data selection
#---------------
biom.macro_16_17 <- marben_biom_before_2016 %>% 
  filter(annee %in% c(2016, 2017)) %>%
  filter (variable %in% c("biomasse algues rouges", "biomasse algues vertes", "biomasse algues brunes")) %>%
  filter (collection %in% c("Biométrie 1", "Biométrie 2")) %>%
  dplyr :: select(-pied, -units, -broken)

biom.macro_19_20 <- biom19.20.cl %>% 
  filter (variable %in% c("biomasse algues rouges", "biomasse algues vertes", "biomasse algues brunes")) %>%
  filter (collection %in% c("1/2 quadrat 1", "1/2 quadrat 2")) %>%
  dplyr :: select(-pied, -feuille)

# Sum of the three algae types per 1/2 quadrat
#---------------------------------------------
biom.macro_16_17 <- biom.macro_16_17%>%
  group_by(theme, site, point_name, date, annee, saison, collection) %>% 
  summarise_at(vars(value), sum, na.rm = TRUE) %>% 
  dplyr::rename(biom.macro = value) %>%
  ungroup()

biom.macro_19_20 <- biom.macro_19_20 %>% 
  group_by(theme, site, point_name, date, annee, saison,collection) %>% 
  summarise_at(vars(value), sum, na.rm = TRUE) %>% 
  dplyr::rename(biom.macro = value) %>%
  ungroup()

# Conversion in m2
#-----------------
biom.macro_16_17 <- biom.macro_16_17%>%
mutate(biom.macro = biom.macro/0.1)

biom.macro_19_20<- biom.macro_19_20 %>% 
mutate(biom.macro = biom.macro/0.1)

# Average macroalgae biomass per point
#-------------------------------------
biom.macro.moy_16_17 <- biom.macro_16_17 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(biom.macro), mean, na.rm = TRUE) %>% 
  dplyr::rename(biom.macro.moy = biom.macro) %>%
  ungroup()

biom.macro.moy_19_20 <- biom.macro_19_20 %>% 
  group_by(theme, site, point_name, date, annee, saison) %>% 
  summarise_at(vars(biom.macro), mean, na.rm = TRUE) %>% 
  dplyr::rename(biom.macro.moy = biom.macro) %>%
  ungroup()
```

### Binding

```{r}
biom.macro.moy <-
  rbind(biom.macro.moy_bf_16, biom.macro.moy_16_17, biom.macro.moy_19_20)

# Plot of macroalgae biomass
ggplot(data = biom.macro.moy, aes(x = date, y = biom.macro.moy)) + geom_point(aes(color = annee)) + 
facet_wrap(~ site)

# Boxplot 
ggplot(biom.macro.moy, aes(x = site, y = biom.macro.moy)) +
  geom_boxplot()
```

## Epiphyte biomasses

### Before 2019

```{r}
# Data selection
#---------------
biom.epi_bf_19<-marben_biom_before_2016 %>% 
  filter(annee %in% c(2004:2018)) %>%
  filter(variable == "Biomasse des épiphytes (Herbiers REBENT)") %>% 
  dplyr :: rename (b.epi = value) %>%
  dplyr :: select(!c(pied, variable, broken, collection, units))

biom.epi.f_bf_19 <-marben_biom_before_2016 %>% 
  filter(annee %in% c(2004:2018)) %>%
  filter(collection == "Epiphytes") %>%
  filter(variable == "Biomasse des feuilles de Zostera marina (Herbier REBENT)") %>%
  dplyr :: rename (b.epi.f = value) %>%
  dplyr :: select(!c(pied, variable, broken, collection, units))

  
# Ratio epiphytes biomass/ leaves biomass
#----------------------------------------
biom.epi.ratio_bf_19 <- biom.epi_bf_19 %>%
  left_join(biom.epi.f_bf_19) %>%
  mutate(ratio.epi = b.epi/b.epi.f)
```

### After 2019

```{r}
# Data selection
#---------------
biom.epi_19_20 <- biom19.20.cl %>% 
  filter(variable == "Biomasse des epiphytes (Herbiers REBENT)") %>%
  dplyr::rename (b.epi = value)%>%
  dplyr :: select(!c(pied, feuille, variable, broken, collection, units))

biom.epi.f_19_20 <- biom19.20.cl %>%
  filter(collection == "Epiphytes") %>%
  filter(variable == "Biomasse des feuilles de Zostera marina (Herbier REBENT)") %>%
  dplyr::rename (b.epi.f = value) %>%
  dplyr :: select(!c(pied, feuille, variable, broken, collection, units))

# Ratio epiphytes biomass/ leaves biomass
#----------------------------------------
biom.epi.ratio_19_20 <- biom.epi_19_20 %>%
  left_join(biom.epi.f_19_20) %>%
  mutate(ratio.epi = b.epi/b.epi.f) 
```

### Binding

```{r}
biom.epi.ratio <- rbind(biom.epi.ratio_bf_19, biom.epi.ratio_19_20)

# Plot
ggplot(data = biom.epi.ratio, aes(x = date, y = ratio.epi)) + 
  geom_point(aes(color = annee)) + 
  facet_wrap(~ site)

# Boxplot 
ggplot(biom.epi.ratio, aes(x = site, y = ratio.epi)) +
  geom_boxplot()
```

# Final dataset for biometry data

```{r}
biom_tab <- density.tot %>%
  left_join(., biomasse) %>%
  left_join(., long.g.moy) %>%
  left_join(., long.f.moy) %>%
  left_join(., larg.f.moy) %>%
  left_join(., n.f.pied.moy) %>%
  left_join(., brk.tot) %>%
  left_join(.,biom.macro.moy ) %>%
  left_join(.,biom.epi.ratio) %>%
  mutate_if(is.character, as.factor) %>%
  dplyr :: select(-date, -b.epi, -b.epi.f)
```

## Visualization of missing biometric values

```{r}
# Points with no missing values
#------------------------------
biom.point.disp <- na.omit(biom_tab) %>%
  distinct(theme, site, point_name, annee, saison) %>% 
  group_by(theme, site, annee, saison) %>%
  tally() %>%
  ungroup()

ggplot(data = biom.point.disp, aes(x = annee, y = site, fill = n)) +
  facet_wrap(vars(saison))+ geom_tile() + geom_text(aes(label = n)) +
  labs(x = "Année", y = "Site", title = "Nombre de points avec toutes les variables disponibles") +
  theme(
    axis.text = element_text(size = 13),
    axis.title = element_text(size = 17, face = "bold"),
    title = element_text(size = 17,face = "bold")) +
  scale_fill_viridis(option = "turbo", begin = 0.3,end = 0.7)

# NA count per site and per variables
#-------------------------------------

# Spring
nbr.NA.print <- biom_tab %>%
  filter(saison == "Printemps") %>%
  dplyr :: select(-theme, -annee, -saison, -point_name)

gg_miss_var(nbr.NA.print, facet = site) +
  labs(y = "Nombre de NA", x = "Variables", title = "Nombre de NA par variable et par site au Printemps")

# Autumne
nbr.NA.autom <- biom_tab %>%
  filter(saison == "Automne") %>%
  dplyr :: select(-theme, -annee, -saison, -point_name)

gg_miss_var(nbr.NA.autom, facet = site) +
  labs(y = "Nombre de NA", x="Variables", title="Nombre de NA par variable et par site en Automne")

# Number of variables available per site and per year
#-----------------------------------------------------
biom.var.disp <- biom_tab %>%
  gather(-theme, -site, -point_name, -annee, -saison, key = "variable", value = "value") %>%
  na.omit() %>%
  distinct(theme, site, point_name, annee, saison, variable) %>%
  group_by(theme, site, annee, saison)%>%
  tally()%>%
  ungroup()

ggplot(data = biom.var.disp, aes(x = annee, y = site, fill = n))+
  facet_wrap(vars(saison)) + geom_tile() + geom_text(aes(label = n)) +
  labs(x = "Année", y="Site", title="Nombre de variables disponibles") +
  theme(
    axis.text = element_text(size = 13),
    axis.title = element_text(size = 17, face = "bold"),
    title = element_text(size = 17, face = "bold")) +
  scale_fill_viridis(option = "turbo", begin = 0.3, end = 0.7)

# Missing variable for each point with NA
#----------------------------------------
#biom_tab_filter<-biom_tab %>%
  #filter(site=="Arradon"& annee %in% c(2005,2011,2013) & saison=="Printemps")
  #filter(site=="Arcouest"& annee %in% C(2005,2011) & saison=="Printemps")
  #filter(site=="Callot"& annee %in% c(2005,2007,2011) & saison=="Printemps")
  #filter(site=="Glenan"& annee %in% c(2005,2006,2008,2011,2014)) & saison=="Printemps")
  #filter(site=="Molene"& annee %in% c(2008,2011,2013,2014,2015,2016) &  saison=="Printemps")
  #filter(site=="Roscanvel"& annee %in% c(2005,2008,2011) &  saison=="Printemps")
  #filter(site=="Saint-Malo"& annee %in% c(2005,2011) &  saison=="Printemps")
  #filter(site=="Sainte-Marguerite"& annee == 2007 &  saison=="Printemps")
  #filter(site=="Sainte-Marguerite"& annee %in% c(2005,2006,2007,2008,2011) &  saison=="Printemps")
  #filter(site=="Sept-Iles"& annee %in% c(2005,2006,2011,2020) &  saison=="Printemps")

# Visualisation of commetn assigned to these points
#--------------------------------------------------

# Data selection
comment <- read.csv("Input/commentaires_HI-1.csv", sep = ",", dec = ".") %>%
  mutate_if(is.character, as.factor)

# Formating date to remove hours
comment$date <- as.POSIXct(comment$date , format = "%Y-%m-%d")
# Variable "date"
comment$annee <- as.factor(format.Date(comment$date,"%Y"))
# Variable "month"
comment$mois <- as.numeric(format.Date(comment$date,"%m"))
# Variable "season"
comment$saison[comment$mois<7] <- "Printemps"
comment$saison[comment$mois>=7] <- "Automne"
comment$saison<- as.factor(comment$saison)

# Comments for each interesting point
#comment_filter <- comment %>%
  #filter(site=="Arradon Herbier"& annee %in% c(2005,2011,2013) & saison=="Printemps")
  #filter(site=="L'Arcouest Herbier" & annee %in% c(2005,2011) & saison=="Printemps")
  #filter(site=="Callot Herbier"& annee %in% c(2005,2007,2011) & saison=="Printemps")
  #filter(site=="Glénan Herbier"& annee %in% c(2005,2006,2008,2011,2014) & saison=="Automne")
  #filter(site=="Molène Herbier" & annee %in% c(2008,2011,2013,2014,2015,2016) &  saison=="Printemps")
  #filter(site=="Roscanvel Herbier"& annee %in% c(2005,2008,2011) &  saison=="Printemps")
  #filter(site=="Saint-Malo Herbier"& annee %in% c(2005,2011) &  saison=="Printemps")
  #filter(site=="Sainte-Marguerite Herbier"& annee %in% c(2005,2006,2007,2008,2011) &      saison=="Printemps")
 #filter(site=="Sept-Iles Herbier"& annee %in% c(2005,2006,2011,2020) &  saison=="Printemps")
```

## Extrapolation or removal of NA

From this stage, work is only carried out for Spring.

A point for which all biometry variables are missing is deleted. For points with only one or two missing biometry variables, the values are extrapolated using the nearest-neighbor method, using only the biometry data. 

```{r}
biom_noNA <- biom_tab %>%
  filter(saison == "Printemps") %>%
  #Variables renamed
  dplyr :: rename (densite = density) %>% 
  #Extrapolation by nearest-neighbor
  kNN(variable = "long.g.moy", imp_var = F, k = 5) %>%
  kNN(variable = "biom.macro.moy", imp_var = F, k = 5) %>%
  kNN(variable = "ratio.epi", imp_var = F, k = 5) %>%
  kNN(variable = "br.m2", imp_var = F, k = 5) %>%
  #Points with all biometry variables missing are deleted 
  filter(site != "Toulindac") %>%
  filter(annee != 2011 & annee != 2004) %>%
  filter(!(site == "Arradon"& annee %in% c(2005, 2013))) %>%
  filter(!(site == "Callot" & annee == 2005))%>%
  filter(!(site == "Glenan" & annee %in% c(2005, 2006, 2008, 2014))) %>%
  filter(!(site == "Molene" & annee == 2005)) %>%
  filter(!(point_name == "Molène 2" & annee %in% c(2013, 2015, 2016))) %>%
  filter(!(point_name == "Molène 3" & annee %in% c(2008, 2013, 2014, 2015, 2016))) %>%
  filter(!(site == "Roscanvel" & annee == 2008)) %>%
  filter(!(site =="Sainte-Marguerite" & annee %in% c(2005, 2006)) ) %>%
  filter(!(point_name %in% c("Sainte-Marguerite 2", "Sainte-Marguerite 3") & annee == 2008)) %>%
  filter(!(site=="Sept-Iles" & annee %in% c(2005, 2006, 2020)))
```

## Separation of environment matrix and metadata 

```{r}
biom_print_data <- biom_noNA %>%
  as.data.frame() %>%
  set_rownames(paste(.$site, .$point_name, .$annee, sep = "_"))%>%
  dplyr :: select(-theme, -site, -point_name, -annee, -saison)

biom_print_meta <- biom_noNA %>%
  dplyr :: select(site, point_name, annee)
```

## PCA biometry

```{r}
# PCA
#----
pca_biom <-biom_print_data %>% 
  dudi.pca (scannf = FALSE, nf = 2) 

#Point coordinates
#-----------------
point_scores_biom <-scores(pca_biom, display = "sites", scaling = 1)%>%
  as.data.frame() %>%
  rownames_to_column("location") %>%
  separate(location, c("site","point_name","annee"), sep="_") %>%
  right_join(biom_print_meta,.)

# Get centroids for each point
#----------------------------
point_scores_cent_biom <- point_scores_biom %>%
  group_by(site)%>% 
  summarise_at(vars(Axis1, Axis2), mean, na.rm=TRUE) %>%
  dplyr :: rename(cent1 = Axis1)%>%
  dplyr :: rename(cent2 = Axis2)%>%
  ungroup %>%
  left_join(point_scores_biom,.)

# Get centroids for each labels
#------------------------------
label_biom <- point_scores_cent_biom %>%
  group_by(site) %>%
  summarise_at(vars(Axis1, Axis2),mean, na.rm=TRUE) %>%
  dplyr :: rename(cent1 = Axis1)%>%
  dplyr :: rename(cent2 = Axis2)%>%
  ungroup 

# Variance for each axis
#-----------------------
var_axes_biom <- round(pca_biom$eig/sum(pca_biom$eig)*100, 2)

# Plot of individuals
#--------------------
p <- ggplot(point_scores_biom, aes(x = Axis1, y = Axis2, col = site, fill = site))
p <- p + geom_point(size = 4, shape = 21, alpha = 0.7) + geom_point(data = point_scores_cent_biom, aes(x = cent1, y = cent2), size = 4, col = "black", shape=17 )
p <- p + geom_label_repel(data = label_biom, aes(x = cent1, y = cent2, label = site, fill = site), fontface = "bold", size = 6, segment.colour = "black", box.padding = 1, inherit.aes  =FALSE, show.legend = FALSE, alpha = 0.8, nudge_x = -0.4, nudge_y = 0.2)
p <- p + xlab(paste("PCA1:", var_axes_biom[1], "%")) + ylab(paste("PCA2:", var_axes_biom[2], "%"))
p <- p + theme_bw() + theme(text = element_text(size = 20), strip.text.x = element_text(face = "bold"), legend.position = "none")
p

#Plot of variables 
#--------------
fviz_pca_var(pca_biom, repel = TRUE, col.var = factor(c("Biométrie", "Biométrie", "Biométrie","Biométrie","Biométrie","Biométrie","Biométrie","Biométrie","Externe","Externe")), legend.title = list(color = "Catégories"), palette = c("#8fd175","#56B4E9")) + theme(axis.title = element_text(size = 16))
```

## Correlation matrix

```{r}
correl_biom <- biom_print_data%>%
  cor(.,method = "pearson", use = "pairwise.complete.obs")
corrplot(correl_biom,method = "number", type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
```

# Fauna data

## Load data

```{r}
load("Input/Faune03.20_19_01_2022.Rdata")
name_accepted <- read.csv("Input/name_accepted.csv", sep = ";")
```

## Formating

```{r}
# Duplicated deleted 
#-------------------
faune <- faune03.20 %>%
  # Sum of abundance for duplicated species
  group_by(theme, site, point_name, replicat, method_echant, annee, saison, date, species) %>%
  summarise_at(vars(abondance), sum, na.rm = TRUE) %>%
  ungroup() %>%
  # "Azoiques" are removed
  spread(species, abondance, fill = 0) %>%
  dplyr:: select(-"Azoïque azoïque") %>%
  gather(species, abondance, -theme, -site, -point_name, -replicat, -method_echant, -annee, -saison, -date)
```

## Infauna

### Aggregation at point scale

```{r}
# Sum of the three replicates abundances per point
#-------------------------------------------------
endofaune_print <- faune %>% filter(abondance != 0) %>%
  left_join(name_accepted) %>%
  filter(theme == "Herbiers Intertidaux" & saison == "Printemps" & method_echant == "Carotte")%>%
  #Point with issues of preservation are deleted
  filter(!(point_name == "Arradon 1" & annee == 2009)) %>%
  filter(!(point_name == "Molene 3" & annee == 2008)) %>%
  filter(!(point_name == "Sept-Iles 2" & annee == 2009)) %>%
  group_by(theme, site, method_echant, saison, date, annee, point_name, scientificname_accepted)%>%
  summarise_at(vars(abondance), sum, na.rm = TRUE) %>%
  ungroup()%>%
  spread(scientificname_accepted, abondance, fill = 0) %>%
  as.data.frame()

# Separation of infauna matrix and metadata
#-------------------------------------------

endofaune_print_data <- endofaune_print %>%
  set_rownames(paste(.$site, .$point_name, .$annee, sep = "_"))%>%
  dplyr :: select(-theme, -site, -point_name, -annee, -saison, -method_echant, -date)

endofaune_print_meta <- endofaune_print %>%
  dplyr :: select(site, point_name, annee)%>%
  mutate(annee = as.factor(annee))
```

### PCA infauna

```{r}
# PCA 
# ---
pca_endofaune <-endofaune_print_data %>% 
  decostand(method = "hellinger")%>%
  dudi.pca (scannf = FALSE, scale  = FALSE, nf = 2) 

# Points coordinates
#-------------------
point_scores_endofaune <-scores(pca_endofaune, display = "sites", scaling = 1)%>%
  as.data.frame() %>%
  rownames_to_column("location") %>%
  separate(location, c("site", "point_name", "annee"), sep = "_") %>%
  right_join(endofaune_print_meta,.)

# Get centroids for each point
#----------------------------
point_scores_cent_endofaune <- point_scores_endofaune %>%
  group_by(site)%>% 
  summarise_at(vars(Axis1, Axis2),mean,na.rm=TRUE) %>%
  dplyr :: rename(cent1=Axis1)%>%
  dplyr :: rename(cent2=Axis2)%>%
  ungroup %>%
  left_join(point_scores_endofaune,.)

# Get centroids for each labels
#------------------------------
label_endofaune <- point_scores_cent_endofaune %>%
  group_by(site) %>%
  summarise_at(vars(Axis1, Axis2),mean,na.rm=TRUE) %>%
  dplyr :: rename(cent1=Axis1)%>%
  dplyr :: rename(cent2=Axis2)%>%
  ungroup 

# Variance for each axis
#-----------------------
var_axes_endofaune <- round(pca_endofaune$eig/sum(pca_endofaune$eig)*100,2)

# Plot of individuals
#--------------------
p <-ggplot(point_scores_endofaune,aes(x = Axis1, y = Axis2, col = site, fill = site, group = site))

p <- p + geom_point(aes(color = site), size = 4, alpha = 0.7) + 
  #scale_shape_manual(values=c(0,1,4,3,5,6,15,21,23)) + 
  scale_fill_manual(values = c("#25FDE9", "black", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855")) +
  scale_color_manual(values = c("#25FDE9", "black", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855")) + 
  geom_point(data = point_scores_cent_endofaune, aes(x = cent1, y = cent2), size = 4, shape = 17, col = "black")

p <- p + geom_label_repel(data = label_endofaune, aes(x = cent1, y = cent2,label = site),col = c("#25FDE9", "black", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855"), fill = "white", fontface = "bold", size = 10, segment.colour = "black", box.padding = 1, inherit.aes = FALSE, show.legend = FALSE, alpha = 0.8, nudge_x = 0.3, nudge_y = 0.3, label.size = 1) 

p <- p + xlab("PCA 1: 25.7 %") + ylab("PCA 2: 11.81 %") + scale_y_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) +
scale_x_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.'))
pca_ind_infauna_plot <- p + theme_bw() + theme(
  text=element_text(size=30),
  strip.text.x = element_text(face="bold"),
  legend.position="none",
  axis.title=element_text(size=30,face="bold"))
pca_ind_infauna_plot

#Plot of variables 
#--------------
pca_var_infauna_plot <- fviz_pca_var(pca_endofaune, repel = TRUE,select.var = list(contrib = 5), labelsize = 6) + xlab("PCA 1: 25.7 %") + ylab("PCA 2: 11.81 %") + theme_bw() + scale_y_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) + scale_x_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.'))+ theme(text=element_text(size = 20), legend.title=element_text(size = 0), axis.title=element_text(size = 20, face = "bold")) + ggtitle("")
pca_var_infauna_plot

```

### Diversity index for infauna

```{r}
#Abundance per point
#--------------------
abondance_endo <- faune%>%
  filter(theme == "Herbiers Intertidaux" & saison == "Printemps" & method_echant == "Carotte") %>%
  filter(!(point_name == "Arradon 1" & annee==2009)) %>%
  filter(!(point_name == "Molene 3" & annee==2008)) %>%
  filter(!(point_name == "Sept-Iles 2" & annee==2009)) %>%
  droplevels() %>%
  group_by(theme, site,method_echant, saison, date, annee, point_name) %>%
  summarise_at(vars(abondance), sum, na.rm = TRUE) %>% ungroup()%>% mutate(annee = as.factor(annee))

#Hill index
#----------
hill_0 <- diversity(as.matrix(endofaune_print_data), type = "td", q = 0)%>%
  dplyr :: rename(hill_0 = hill.numbers)
hill_1 <- diversity(as.matrix(endofaune_print_data), type = "td", q=1)%>%
  dplyr :: rename(hill_1 = hill.numbers)
hill_2 <- diversity(as.matrix(endofaune_print_data), type = "td", q=2)%>%
  dplyr :: rename(hill_2 = hill.numbers)

# Additional diversity index
#---------------------------
diversity_index_endo <-diversity(as.matrix(endofaune_print_data), q = 0) %>%
  dplyr :: select(entropy, variety, simpson.I, evenness) %>%
  dplyr :: rename (rich_sp = variety) %>%
  dplyr :: rename (shannon = entropy) %>%
  rownames_to_column("location") %>%
  separate(location, c("site", "point_name", "annee"), sep="_") %>%
  inner_join(abondance_endo) %>%
  dplyr :: select(-theme, -date, -saison, -method_echant)

#Plots
#----

#Abundance 
ggplot(diversity_index_endo, aes(x = annee, y = abondance))+
  labs(x = "Année", y = "Abondance") +
  geom_point(aes(color = annee), size = 3)+
  facet_wrap(vars(site),scale = "free")+
  theme(
    axis.text.x = element_text(size = 8, angle=90),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 15, face = "bold"),
    strip.text = element_text(size = 12, face = "bold"),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12))

#Species richness
ggplot(diversity_index_endo, aes(x = annee, y = rich_sp)) +
  labs(x = "Année", y = "Richesse spécifique") +
  geom_point(aes(color = annee), size = 3) +
  facet_wrap(vars(site),scale = "free") +
  theme(
    axis.text.x = element_text(size = 8, angle = 90),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 15, face = "bold"),
    strip.text = element_text(size = 12, face = "bold"),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12))

#Shannon
ggplot(diversity_index_endo, aes(x = annee, y = shannon)) +
  labs(x = "Année", y = "Indice de Shannon") +
  geom_point(aes(color = annee), size = 3) +
  facet_wrap(vars(site), scale = "free") +
  theme(
    axis.text.x = element_text(size = 8, angle = 90),
    axis.text.y=element_text(size = 12),
    axis.title = element_text(size = 15, face = "bold"),
    strip.text = element_text(size = 12, face = "bold"),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12))

#Simpson
ggplot(diversity_index_endo, aes(x = annee, y = simpson.I))+
  labs(x = "Année", y = "Indice de Simpson") +
  geom_point(aes(color = annee), size = 3) +
  facet_wrap(vars(site),scale = "free") +
  theme(
    axis.text.x = element_text(size = 8, angle = 90),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 15,face = "bold"),
    strip.text = element_text(size = 12, face = "bold"),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12))

#Pielou
ggplot(diversity_index_endo, aes(x = annee,y = evenness)) +
  labs(x = "Année", y = "Indice de Piélou") +
  geom_point(aes(color = annee), size = 3) +
  facet_wrap(vars(site), scale = "free") +
  theme(
    axis.text.x = element_text(size = 8, angle=90),
    axis.text.y = element_text(size = 12), 
    axis.title = element_text(size = 15,face = "bold"), 
    strip.text = element_text(size = 12, face = "bold"),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12))

#Boxplots
#--------

#Species richness
ggplot(diversity_index_endo, aes(x = point_name, y = rich_sp)) +
  labs(x = "Point", y = "Richesse spécifique") +
  geom_boxplot(aes(color = site), scale = "free") +
  theme(
    axis.text.x = element_text(size = 15, angle = 90),
    axis.text.y = element_text(size = 12), 
    axis.title = element_text(size = 20,face = "bold"),
    strip.text = element_text(size = 12, face = "bold"),
    legend.title = element_text(size = 15), 
    legend.text = element_text(size = 12))

#Shannon
ggplot(diversity_index_endo, aes(x = point_name, y = shannon)) +
 labs(x = "Point", y = "Indice de Shannon") +
  geom_boxplot(aes(color = site), scale = "free")+
  theme(
    axis.text.x = element_text(size = 15, angle=90),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 20, face = "bold"),
    strip.text = element_text(size = 12, face = "bold"),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12))

#Simpson
ggplot(diversity_index_endo, aes(x = point_name, y = simpson.I)) +
  labs(x = "Point", y = "Indice de Simpson") +
  geom_boxplot(aes(color = site), scale = "free") +
  theme(
    axis.text.x = element_text(size = 15, angle = 90),
    axis.text.y = element_text(size = 12),
    axis.title=element_text(size = 20, face = "bold"),
    strip.text = element_text(size = 12, face = "bold"),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12))
```

### Correlation matrix for infauna

```{r}
correl_endofaune <- diversity_index_endo %>%
  dplyr :: select ( -site, -point_name, -annee) %>% 
  cor(.,method = "pearson", use = "pairwise.complete.obs")
corrplot(correl_endofaune, method = "number", type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
```

## Epifauna

### Aggregation at point scale

```{r}
# Sum of the three replicates abundances per point
#-------------------------------------------------
epifaune_print <- faune %>% filter(abondance != 0) %>%
  left_join(name_accepted) %>%
  filter(theme == "Herbiers Intertidaux" & saison == "Printemps" & method_echant == "Haveneau") %>%
  filter(!(site == "Arradon" & annee == 2009)) %>%
  group_by(theme, site, method_echant, saison, date, annee, point_name, species) %>%
  summarise_at(vars(abondance), sum, na.rm = TRUE) %>%
  ungroup() %>%
  spread(species, abondance, fill = 0) %>%
  as.data.frame()

# Separation of infauna matrix and metadata
#-------------------------------------------
epifaune_print_data <- epifaune_print %>%
  set_rownames(paste(.$site, .$point_name, .$annee, sep = "_")) %>%
  dplyr :: select(-theme, -site, -point_name, -annee,-saison, -method_echant, -date)

epifaune_print_meta <- epifaune_print %>%
  dplyr :: select(site, point_name, annee) %>%
  mutate(annee = as.factor(annee))
```

### PCA epifauna

```{r}
#PCA
#---
pca_epifaune <- epifaune_print_data %>%
  decostand(method = "hellinger") %>%
  dudi.pca (scannf = FALSE, scale=FALSE, nf = 2) 

# Points coordinates
#-------------------
point_scores_epifaune <-scores(pca_epifaune, display = "sites", scaling = 1)%>%
  as.data.frame() %>%
  rownames_to_column("location") %>%
  separate(location,c("site", "point_name", "annee"), sep = "_") %>%
  right_join(epifaune_print_meta, .)

# Get centroids for each point
#----------------------------
point_scores_cent_epifaune <- point_scores_epifaune %>%
  group_by(site) %>% 
  summarise_at(vars(Axis1, Axis2),mean, na.rm = TRUE) %>%
  dplyr :: rename(cent1 = Axis1) %>%
  dplyr :: rename(cent2 = Axis2) %>%
  ungroup %>%
  left_join(point_scores_epifaune, .)

# Get centroids for each labels
#------------------------------
label_epifaune <- point_scores_cent_epifaune %>%
  group_by(site) %>%
  summarise_at(vars(Axis1, Axis2), mean, na.rm=TRUE) %>%
  dplyr :: rename(cent1 = Axis1) %>%
  dplyr :: rename(cent2 = Axis2) %>%
  ungroup 

# Variance for each axis
#-----------------------
var_axes_epifaune <- round(pca_epifaune$eig/sum(pca_epifaune$eig)*100, 2)

#Plot of individuals
#------------------

#Calcul de l'ACP 
#---------------
pca_epifaune <- epifaune_print_data %>%
  decostand(method = "hellinger") %>%
  dudi.pca (scannf = FALSE, scale=FALSE, nf = 2) 

#On récupère les coordonnées des points
#--------------------------------------
point_scores_epifaune <-scores(pca_epifaune, display = "sites", scaling = 1)%>%
  as.data.frame() %>%
  rownames_to_column("location") %>%
  separate(location,c("site", "point_name", "annee"), sep = "_") %>%
  right_join(epifaune_print_meta, .)

#On récupère les centroides de chaque point
#------------------------------------------
point_scores_cent_epifaune <- point_scores_epifaune %>%
  group_by(site) %>% 
  summarise_at(vars(Axis1, Axis2),mean, na.rm = TRUE) %>%
  dplyr :: rename(cent1 = Axis1) %>%
  dplyr :: rename(cent2 = Axis2) %>%
  ungroup %>%
  left_join(point_scores_epifaune, .)

# On récupère les centroides pour les labels
#-------------------------------------------
label_epifaune <- point_scores_cent_epifaune %>%
  group_by(site) %>%
  summarise_at(vars(Axis1, Axis2), mean, na.rm=TRUE) %>%
  dplyr :: rename(cent1 = Axis1) %>%
  dplyr :: rename(cent2 = Axis2) %>%
  ungroup 

# Calcul de la variance représentée par chaque axe
#-------------------------------------------------
var_axes_epifaune <- round(pca_epifaune$eig/sum(pca_epifaune$eig)*100, 2)

#Plot individus
#--------------

options(ggrepel.max.overlaps = Inf)
p <-ggplot(point_scores_granul_hydro, aes(x = Axis1, y = Axis2, col = site, fill = site, group = site))

p <- p + geom_point(aes(color = site), size = 3, alpha = 0.7) + 
  #scale_shape_manual(values=c(0,1,4,3,5,6,15,21,23)) + 
  scale_fill_manual(values = c("#25FDE9", "black", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855")) +
  scale_color_manual(values = c("#25FDE9", "black", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855")) +
  geom_point(data = point_scores_cent_granul_hydro, aes(x = cent1, y = cent2), size = 4, shape = 17, col = "black")

p <- p <- p + geom_label_repel(data = label_granul_hydro, aes(x = cent1, y = cent2,label = site), col = c("#25FDE9", "black", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855"), fill = "white", fontface = "bold",size = 10, segment.colour = "black", box.padding = 1, inherit.aes = FALSE, show.legend = FALSE, alpha = 0.8, nudge_x = 0.3, nudge_y = 0.3, label.size = 1)

p <- p + xlab(paste("PCA 1: 31.08 %")) + ylab(paste("PCA 2: 23.63 %"))

pca_ind_env_plot <- p + scale_y_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) +
scale_x_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) + 
  theme_bw() + theme(
    text=element_text(size=30), 
    strip.text.x = element_text(face = "bold"),
    legend.position="none",
    axis.title=element_text(size=30,face="bold"))
pca_ind_env_plot

# Plot variables
#---------------

pca_var_epifauna_plot <- fviz_pca_var(pca_epifaune, repel = TRUE, select.var = list(contrib=5), labelsize=6) + xlab(paste("PCA 1: 21.59 %")) + ylab(paste("PCA 2: 13.4 %")) +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) +
  scale_x_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) +
  theme_bw() + 
  theme(
    text = element_text(size = 20),
    legend.title = element_text(size = 0),
    axis.title = element_text(size = 20, face = "bold")) + ggtitle("")
  
pca_var_env_plot <- fviz_pca_var(pca_granul_hydro, repel = TRUE, labelsize = 4) + 
  xlab("PCA 1: 31.08 %") + ylab("PCA 2: 23.63 %") +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) +
  scale_x_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) + 
  ggtitle("") + theme_bw() + 
  theme(
  text = element_text(size = 15),
  legend.title = element_text(size = 0),
  legend.text = element_text(size = 6), 
  legend.position = "bottom",
  axis.title = element_text(size = 15,face = "bold"))
pca_var_env_plot 

```

### Diversity index for epifauna

```{r}
#Abundances per point
#--------------------
abondance_epi <- faune%>%
  filter(theme == "Herbiers Intertidaux" & saison == "Printemps" & method_echant == "Haveneau") %>%
  filter(!(site == "Arradon" & annee == 2009)) %>%
  droplevels()%>%
  group_by(theme, site, method_echant, saison, date, annee, point_name) %>%
  summarise_at(vars(abondance), sum, na.rm = TRUE) %>% ungroup() %>%
  mutate(annee = as.factor(annee))

#Hill index
#----------
hill_0_epi <- diversity(as.matrix (epifaune_print_data), type = "td", q = 0)%>%
  dplyr :: rename(hill_0 = hill.numbers)
hill_1_epi <- diversity(as.matrix (epifaune_print_data), type = "td", q = 1)%>%
  dplyr :: rename(hill_1 = hill.numbers)
hill_2_epi <- diversity(as.matrix (epifaune_print_data), type = "td", q = 2)%>%
  dplyr :: rename(hill_2 = hill.numbers)

#Additional diversity index
#--------------------------
diversity_index_epi <- diversity(as.matrix (epifaune_print_data), q = 0) %>%
  dplyr :: select(entropy, variety, simpson.I, evenness) %>%
  dplyr :: rename (rich_sp = variety) %>%
  dplyr :: rename (shannon = entropy) %>%
  rownames_to_column("location") %>%
  separate(location, c("site", "point_name", "annee"), sep="_") %>%
  inner_join(abondance_epi) %>%
  mutate(annee = as.factor(annee)) %>%
  dplyr :: select(-theme, -date, -saison, -method_echant)
  
#Plots
#----

#Abundance 
ggplot(diversity_index_epi, aes(x = annee, y = abondance)) +
labs(x = "Année", y = "Abondance") +
geom_point(aes(color = annee), size = 3)+
facet_wrap(vars(site), scale = "free")+
theme(
  axis.text.x = element_text(size = 8, angle = 90),
  axis.text.y = element_text(size = 12),
  axis.title = element_text(size = 15, face = "bold"),
  strip.text = element_text(size = 12, face = "bold"),
  legend.title = element_text(size = 15),
  legend.text = element_text(size = 12))

#Species richness
ggplot(diversity_index_epi, aes(x = annee, y = rich_sp)) +
labs(x = "Année", y = "Richesse spécifique") +
geom_point(aes(color = annee), size = 3)+
facet_wrap(vars(site), scale = "free")+
theme(
  axis.text.x = element_text(size = 8, angle=90),
  axis.text.y = element_text(size = 12),
  axis.title = element_text(size = 15, face = "bold"),
  strip.text = element_text(size = 12, face = "bold"),
  legend.title = element_text(size = 15),
  legend.text = element_text(size = 12))

#Shannon
ggplot(diversity_index_epi, aes(x = annee, y = shannon)) +
labs(x = "Année", y = "Indice de Shannon") +
geom_point(aes(color = annee), size = 3)+
facet_wrap(vars(site), scale = "free")+
theme(
  axis.text.x = element_text(size = 8, angle=90),
  axis.text.y = element_text(size = 12),
  axis.title = element_text(size = 15, face = "bold"),
  strip.text = element_text(size = 12, face = "bold"),
  legend.title = element_text(size = 15),
  legend.text = element_text(size = 12))

#Simspon
ggplot(diversity_index_epi, aes(x = annee, y = simpson.I)) +
labs(x = "Année", y = "Indice de Simpson") +
geom_point(aes(color = annee), size = 3)+
facet_wrap(vars(site), scale = "free")+
theme(
  axis.text.x = element_text(size = 8, angle=90),
  axis.text.y = element_text(size = 12),
  axis.title = element_text(size = 15, face = "bold"),
  strip.text = element_text(size = 12, face = "bold"),
  legend.title = element_text(size = 15),
  legend.text = element_text(size = 12))

#Piélou
ggplot(diversity_index_epi, aes(x = annee, y = evenness)) +
labs(x = "Année", y = "Indice de Piélou") +
geom_point(aes(color = annee), size = 3)+
facet_wrap(vars(site), scale = "free")+
theme(
  axis.text.x = element_text(size = 8, angle=90),
  axis.text.y = element_text(size = 12),
  axis.title = element_text(size = 15, face = "bold"),
  strip.text = element_text(size = 12, face = "bold"),
  legend.title = element_text(size = 15),
  legend.text = element_text(size = 12))

#Boxplots
#--------

#Species richness
ggplot(diversity_index_epi, aes(x = point_name,y = rich_sp)) +
geom_boxplot(aes(color = site)) +
labs(x = "Point", y = "Richesse spécifique") +
theme(
  axis.text.x = element_text(size = 15, angle = 90), 
  axis.text.y = element_text(size = 12), 
  axis.title=element_text(size = 20,face = "bold"), 
  strip.text = element_text(size = 12, face = "bold"),
  legend.title = element_text(size = 15),
  legend.text = element_text(size = 12))

#Shannon
ggplot(diversity_index_epi, aes(x = point_name, y = shannon)) +
labs(x = "Point", y="Indice de Shannon") +
geom_boxplot(aes(color=site))+
theme(
  axis.text.x = element_text(size = 15, angle = 90),
  axis.text.y = element_text(size = 12),
  axis.title = element_text(size=20,face="bold"), 
  strip.text = element_text(size = 12, face = "bold"), 
  legend.title = element_text(size = 15),
  legend.text = element_text(size = 12))

#Simpson
ggplot(diversity_index_epi, aes(x = point_name, y = simpson.I)) +
labs(x = "Point", y = "Indice de Simpson") +
geom_boxplot(aes(color = site))+
theme(
  axis.text.x = element_text(size = 15, angle = 90),
  axis.text.y = element_text(size = 12),
  axis.title = element_text(size=20,face="bold"), 
  strip.text = element_text(size = 12, face = "bold"), 
  legend.title = element_text(size = 15),
  legend.text = element_text(size = 12))
```

### Correlation matrix for epifauna

```{r}
correl_epifaune <-diversity_index_epi %>%
  dplyr :: select (-site, -point_name, -annee) %>% 
  cor(.,method="pearson", use="pairwise.complete.obs")

corrplot(correl_epifaune, method = "number", type = "upper",order = "hclust", tl.col = "black", tl.srt = 45)
```

### Abundance grazers 

```{r}
# Abundance grazers
#------------------

abond_brouteur_epi <- faune_trophique %>%
  filter(method_echant == "Haveneau") %>%
  group_by(site, point_name, annee,method_echant, description) %>%
  summarise_at(vars(abondance), sum) %>% ungroup %>%
  filter(description == "Microbrouteurs") %>%
  dplyr :: rename(abond_brouteur = abondance) %>%
  dplyr :: select (site, point_name, annee, abond_brouteur)
```

# Granulomerty

## Load data

```{r}
load("Input/granulo03.20_21_01_2022.Rdata")
```

## Formating

```{r}
granul_point <- granulo_MO_all_point %>%
  filter(theme == "Herbiers Intertidaux" & saison == "Printemps") %>%
  mutate(annee = as.factor(annee)) %>%
  group_by(site, point_name,annee) %>%
  summarise_if(is.numeric,mean) %>% ungroup() %>% clean_names() %>%
  dplyr :: rename(distribution_moy = mean_fw_um, kurtosis = kurtosis_fw_um, distribution_var = d50_um, trask = trask_so, matiere_organique = pourcentage_mo, vase = mud, gravier = gravel, sable = sand)
```

## Visualisation of missing data for granulometry

```{r}
# Points with all variables available 
#------------------------------------
granul_point_disp<-na.omit(granul_point) %>%
  distinct(site, point_name, annee) %>% 
  group_by(site, annee) %>% tally() %>%
  ungroup()

ggplot(data = granul_point_disp, aes(x = annee, y = site, fill = n)) + 
  geom_tile() + geom_text(aes(label = n)) + 
  labs(x = "Année", y = "Site", title="Nombre de points disponible par site et année") +
  theme(
    axis.text = element_text(size = 13),
    axis.title=element_text(size = 17, face = "bold"),
    title=element_text(size = 17, face = "bold")) +
  scale_fill_viridis(option = "turbo", begin = 0.3, end = 0.7)

#Number of variables available per site and per year
#---------------------------------------------------
granul_var <- granul_point %>%
  gather(-site, -point_name, -annee, key = "variable", value = "value") %>%
  na.omit() %>%
  distinct(site, point_name, annee, variable) %>%
  group_by(site, annee) %>%
  tally() %>% ungroup()

ggplot(data = granul_var, aes(x = annee, y = site, fill = n)) + geom_tile() + geom_text(aes(label = n)) + 
  labs(x = "Année", y="Site", title="Nombre de variables disponibles") + 
  theme(
    axis.text = element_text(size = 13),
    axis.title = element_text(size = 17, face = "bold"), 
    title = element_text(size = 17, face = "bold")) +
  scale_fill_viridis(option = "turbo", begin = 0.3, end = 0.7)
```

## Extrapolation or removal of NA 

A point for which all grain size variables are missing is deleted. For points with only one or two missing granulometry variables, the values are extrapolated using the nearest neighbor method, using only the granulometry data.

```{r}
granul_noNA <- granul_point %>%
  #Extrapolation with nearest neighbor
  kNN(variable = "matiere_organique", imp_var = F, k = 5) %>%
  #Points for which all granulometry variables are missinf are deleted
  filter(annee!=2011 & annee !=  2004) %>%
  filter(!(site == "Arradon"& annee %in% c(2005, 2013))) %>%
  filter(!(site == "Callot" & annee == 2005)) %>%
  filter(!(site == "Glenan" & annee %in% c(2005, 2006, 2008, 2014))) %>%
  filter(!(site == "Molene" & annee == 2005))%>%
  filter(!(point_name == "Molène 2" & annee %in% c(2013, 2015, 2016))) %>%
  filter(!(point_name == "Molène 3" & annee %in% c(2008, 2013, 2014, 2015, 2016))) %>%
  filter(!(site == "Roscanvel" & annee == 2008)) %>%
  filter(!(site =="Sainte-Marguerite" & annee %in% c(2005, 2006))) %>%
  filter(!(point_name %in% c("Sainte-Marguerite 2","Sainte-Marguerite 3") & annee == 2008)) %>%
  filter(!(site == "Sept-Iles" & annee %in% c(2005, 2006, 2020)))
```

## Separation of environmental matrix and metadata

```{r}
granul_print_data <- granul_noNA %>%
  as.data.frame()%>%
  set_rownames(paste(.$site, .$point_name, .$annee, sep = "_")) %>%
  dplyr :: select(-site, -point_name, -annee)

granul_print_meta <- granul_noNA %>%
  dplyr :: select(site, point_name, annee)
```

## Check the consistency of extrapolated points

```{r}
#Extrapolation of organic matter percentage data
#-----------------------------------------------

#Glénan
Glenan_granul<-granul_noNA %>%
  dplyr :: select(site, point_name, annee, matiere_organique) %>%
  filter(site == "Glenan")

G1 <- granul_noNA %>%
  filter(site == "Glenan" & point_name == "Glénan 3" & annee == 2010)

g1 <- ggplot(Glenan_granul, aes(x = annee, y = matiere_organique)) +
  labs(x = "Année", y = "matiere organique")+
  geom_point(col="black", size=3)+ geom_point(data=G1, colour="red",size=3) +
  facet_wrap(vars(site))+
  theme(axis.text.x=element_text(size=12,angle=90),axis.text.y=element_text(size=12),
  axis.title=element_text(size=12,face="bold"), strip.text = element_text(size = 15, 
  face = "bold"),legend.title = element_text(size = 20),legend.text = element_text(size =
  16))

#Roscanvel
Roscanvel_granul<-granul_noNA %>%
  dplyr :: select(site,point_name, annee, matiere_organique) %>%
  filter(site=="Roscanvel")

R1<-granul_noNA%>%
  filter(site=="Roscanvel" & point_name %in% c("Roscanvel 1","Roscanvel 2","Roscanvel 3") & annee %in% c(2005,2007)) %>%
  filter(!(site=="Roscanvel" & point_name == "Roscanvel 3" & annee == 2007))

g2<-ggplot(Roscanvel_granul, aes(x=annee, y=matiere_organique))+
  labs(x = "Année", y="pourcentage de MO") +
  geom_point(col="black", size=3)+ geom_point(data=R1, colour="red",size=3) +
  facet_wrap(vars(site)) +
  theme(axis.text.x=element_text(size=12, angle=90),axis.text.y=element_text(size=12), axis.title=element_text(size=12,face="bold"), strip.text = element_text(size = 15, face = "bold"), legend.title = element_text(size = 20), legend.text = element_text(size = 16))

#Arcouest
Arcouest_granul<-granul_noNA%>%
  dplyr :: select(site,point_name, annee, matiere_organique)%>%
  filter(site=="Arcouest")

A1<-granul_noNA%>%
  filter(site=="Arcouest" & point_name %in% c("Arcouest 1","Arcouest 2","Arcouest 3") & annee == 2005)

g3<-ggplot(Arcouest_granul, aes(x=annee, y=matiere_organique)) +
  labs(x = "Année", y="pourcentage de MO") +
  geom_point(col="black", size=3)+ geom_point(data=A1, colour="red",size=3) +
  facet_wrap(vars(site)) +
  theme(axis.text.x=element_text(size=12,angle=90),
axis.text.y=element_text(size=12), axis.title=element_text(size=12,face="bold"), strip.text = element_text(size = 15, face = "bold"),legend.title = element_text(size = 20), legend.text = element_text(size = 16))

#Saint-Malo
St_Malo_granul<-granul_noNA%>%
  dplyr :: select(site,point_name, annee, matiere_organique) %>%
  filter(site=="Saint-Malo")

SMALO1<-granul_noNA%>%
  filter(site=="Saint-Malo" & point_name %in% c("Saint-Malo 2","Saint-Malo 3") & annee == 2005)

g4<-ggplot(St_Malo_granul, aes(x=annee, y= matiere_organique))+
  labs(x = "Année", y="pourcentage de MO") +
  geom_point(col="black", size=3)+ geom_point(data=SMALO1, colour="red",size=3) + facet_wrap(vars(site))+
  theme(axis.text.x=element_text(size=12, angle=90), axis.text.y=element_text(size=12),
        axis.title=element_text(size=12,face="bold"), strip.text = element_text(
        size = 15, face = "bold"),legend.title = element_text(size = 20),
  legend.text = element_text(size = 16))
g4

#Sainte-Marguerite
St_Marg_granul<-granul_noNA%>%
  dplyr :: select(site,point_name, annee, matiere_organique)%>%
  filter(site=="Sainte-Marguerite")

SMARG1<-granul_noNA%>%
  filter(site=="Sainte-Marguerite" & point_name %in% c("Sainte-Marguerite 1","Sainte-Marguerite 2","Sainte-Marguerite 3") & annee == 2007)

g5<-ggplot(St_Marg_granul, aes(x=annee, y= matiere_organique))+
  labs(x = "Année", y="pourcentage de MO")+
  geom_point(col="black", size=3)+ geom_point(data=SMARG1, colour="red",size=3) +
  facet_wrap(vars(site))+
  theme(axis.text.x=element_text(size=12, angle=90),axis.text.y=element_text(size=12),
        axis.title=element_text(size=12,face="bold"), strip.text = element_text(
        size = 15, face = "bold"),legend.title = element_text(size = 20),
  legend.text = element_text(size = 16))
g5

#Plot
grid.arrange(g1,g2,g3,g4,g5, ncol=3, nrow=2)
```

## PCA granulometry

```{r}
#PCA
#---
pca_granul <- granul_print_data %>%
  dplyr :: select(-boulder, -sable, -matiere_organique) %>%
  dudi.pca (scannf = FALSE, nf = 2) 

# Get point coordinates
#----------------------
point_scores_granul <-scores(pca_granul, display = "sites", scaling = 1) %>%
  as.data.frame() %>%
  rownames_to_column("location") %>%
  separate(location, c("site", "point_name", "annee"), sep = "_")

# Get centroides of each points
#------------------------------
point_scores_cent_granul <- point_scores_granul %>%
  group_by(site) %>% 
  summarise_at(vars(Axis1, Axis2), mean, na.rm = TRUE) %>%
  dplyr :: rename(cent1 = Axis1) %>%
  dplyr :: rename(cent2 = Axis2) %>%
  ungroup %>%
  left_join(point_scores_granul, .)

# Get centroides of each labels
#------------------------------
label_granul <- point_scores_cent_granul %>%
  group_by(site) %>%
  summarise_at(vars(Axis1, Axis2),mean, na.rm = TRUE) %>%
  dplyr :: rename(cent1 = Axis1)%>%
  dplyr :: rename(cent2 = Axis2)%>%
  ungroup 

# Variance for each axis
#-----------------------
var_axes_granul <- round(pca_granul$eig/sum(pca_granul$eig)*100,2)

#Plot of individuals
#------------------
p <- ggplot(point_scores_granul, aes(x = Axis1, y = Axis2, col = site, fill = site))
p <- p + geom_point(size = 4,shape = 21, alpha = 0.7) + geom_point(data = point_scores_cent_granul, aes(x = cent1, y = cent2), size = 4, shape = 17, col = "black")
p <- p + geom_label_repel(data = label_granul, aes(x = cent1, y = cent2, label = site, fill = site), fontface = "bold", size = 6, segment.colour = "black", box.padding = 1, inherit.aes = FALSE, show.legend = FALSE, alpha = 0.8, nudge_x = 0.2, nudge_y = -0.2)
p <- p + xlab(paste("PCA1:", var_axes_granul[1], "%")) + ylab(paste("PCA2:", var_axes_granul[2], "%"))
p <- p + theme_bw() + theme(
  text = element_text(size = 20),
  strip.text.x = element_text(face = "bold"),
  legend.position = "none")
p

#Plot of variables
#--------------
fviz_pca_var(pca_granul, repel = TRUE, col.var = factor(c("Indice", "Indice", "Indice", "Indice", "Taille des grains", "Taille des grains")), legend.title = list(color = "Catégories"), palette = c("#56B4E9", "#E69F00")) + theme(axis.title = element_text(size = 16))
```

## Correlation matrix

```{r}
correl_granulo <- granul_print_data %>%
  dplyr :: select(-boulder) %>%
  cor(.,method = "pearson", use = "pairwise.complete.obs")

corrplot(correl_granulo,method = "number", type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
```

## Temporal series

```{r}
#Trask
ggplot(granul_noNA, aes(x = annee, y = trask)) +
  labs(x = "Année", y = "Trask") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"), 
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

#Sand
ggplot(granul_noNA, aes(x = annee, y = sable)) +
  labs(x = "Année", y = "Sable") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"), 
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

#Gravel
ggplot(granul_noNA, aes(x = annee, y = gravier)) +
  labs(x = "Année", y = "Gravier") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"), 
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

#Mud
ggplot(granul_noNA, aes(x = annee, y = vase)) +
  labs(x = "Année", y = "Vase") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"), 
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))


#Organic matter
ggplot(granul_noNA, aes(x = annee, y = matiere_organique)) +
  labs(x = "Année", y = "Matiere organique") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"), 
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

#Kurtosis
ggplot(granul_noNA, aes(x = annee, y = kurtosis)) +
  labs(x = "Année", y = "Kurtosis") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"), 
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))
```

# Environmental data

## Load of data 

```{r}
# Temperature, SST (extraction from surface layer), bottomT, zos and mlotst (daily mean data)
load("Input/copern_phy_one_layer_min_buffer.Rdata")

# Temperature, Salinity data averaged on the total water column (daily mean data)
load("Input/copern_phy_temp_sal_min_buffer.Rdata")

# Daily current averaged from surface uo et vo
load ("Input/copern_phy_curr_min_buffer.Rdata")

# Air temperature, precipitation and wind (daily mean data)
load("Input/Meteo_30_11_2020.Rdata")

# Daily temperature min/max
load("Input/Temp_minmax_22_01_2021.Rdata")

#Fetch
load("Input/Fetch_tot_13-11-20.Rdata")
```

## Formating

Dataframe with variables of interest

```{r}
# Character to factor
copern_phy_curr <-  copern_phy_curr %>% mutate_if(is.character, as.factor)
copern_phy_t_s <-  copern_phy_t_s %>% mutate_if(is.character, as.factor)
copern_phy_one_layer <- copern_phy_one_layer %>% mutate_if(is.character, as.factor)

# Temporal series until 2002
copern_phy_curr <- copern_phy_curr %>% filter(date >= "2002-01-01")
copern_phy_t_s <- copern_phy_t_s %>% filter(date >= "2002-01-01")
copern_phy_one_layer <- copern_phy_one_layer %>% filter(date >= "2002-01-01")

# We retain bottomT + salinity averaged on the water colum + mix layer depth
bottomT_ml <- copern_phy_one_layer %>% 
  filter (variable %in% c("bottomT", "mlotst")) %>% 
  droplevels()
so <- copern_phy_t_s %>% 
  filter(variable =="so") %>% 
  droplevels()

# Rename variables
copern_phy_curr <- copern_phy_curr %>% 
  dplyr::rename (value = daily_curr_mean) %>% 
  mutate (variable = "current") %>%
  mutate_if(is.character, as.factor) %>%
  dplyr::select(value, site_name, variable, date)

# Binding for seawater variables
hydro_daily <- rbind(bottomT_ml, copern_phy_curr, so)

# Variables "theme", "site", "point_name" and "method_echant" are created and only theme = "Herbiers intertidaux" is filtered
hydro_daily_point <- hydro_daily %>%
  separate(site_name, into = c("theme", "site", "method_echant", "point_name"), sep="_") %>%
  dplyr:: select(theme, site, point_name, method_echant, date, variable, value) %>%
  mutate_if (is.character, as.factor) %>%
  filter(theme=="Herbiers Intertidaux") %>%
  spread(variable, value)

# Inspect hydro_daily 
head(hydro_daily_point)
summary(hydro_daily_point)
hydro_daily_point %>% filter(is.na(mlotst)) %>% distinct (theme, site, date)

# Meteo
summary(meteo_points)
summary(temp_minmax_points)
dim(meteo_points)
dim(temp_minmax_points)

# Binding of seawater variables + meteo and filter on "herbiers intertidaux" = intertidal meadows
env_point <- meteo_points %>% 
  left_join(.,temp_minmax_points)%>%
  left_join(hydro_daily_point)%>%
  filter(theme=="Herbiers Intertidaux")
summary(env_point)
dim(env_point)
```

### Season variable added

```{r}
# Year variable
env_point$annee <- as.factor(format.Date(env_point$date,"%Y"))

# Month variable
env_point$mois <- as.numeric(format.Date(env_point$date,"%m"))

# Season variable
env_point$saison[env_point$mois<7] <- "Printemps"
env_point$saison[env_point$mois>=7] <- "Automne"
env_point$saison<- as.factor(env_point$saison)

#Environmental data for spring
env_point_print<-env_point%>%
  filter(saison=="Printemps")


```

### Format fauna data for spring

```{r}
# Extraction of spring samples ID 
#--------------------------------
 
id_echant_print <- faune %>% 
  filter(saison == "Printemps") %>%
  filter(theme == "Herbiers Intertidaux") %>%
  droplevels() %>%
  distinct(theme, site, point_name, method_echant, date, annee) %>%
  mutate_if(is.character, as.factor)

# New column : id_echant_print 
id_echant_print$site_name <- paste(id_echant_print$theme, id_echant_print$site, id_echant_print$point_name, id_echant_print$method_echant, sep = "_")
# Same new column in environmental dataframe
env_point_print$site_name <- paste(env_point_print$theme, env_point_print$site, env_point_print$point_name, env_point_print$method_echant, sep = "_")
```

## Integration 

### Back to 01-11 of the previous year : for spring data

Function for integrating data back to the previous year (integration back to 01-11 of the previous year for spring). And for future work --> data back to 01-01 for autumn data

```{r}
integr <- function(env_point_print, id_echant_print, varnames=c("bottomT", "current", "so", "mlotst"), verbose = T, winter_start = "11-01", summer = F, summer_date = c("07-01","09-30")){
 
  # Package
  require(lubridate)
  require(chron)
  require(dplyr)
  require(magrittr)
  
  # Prevision out
  prev_out <- NULL
  
  # Description
    # id_echant : with sample ID and date
  
  # Check date format
  if (!is.Date(id_echant_print$date)){
     stop("Date not formatted as a date object")
  } else {
    
    # Loop on each sampling site
    for( i in 1:nrow(id_echant_print))
    {
      # Sampling date extraction
      sampling_date <- id_echant_print$date[i]
      
      # Sampling site extraction
      sampling_point <- id_echant_print$site_name[i]
      
      if(verbose){
        print(paste(sampling_point,sampling_date,sep="-"))
      }
      
      # Extraction of the year preceding sampling
      prev_year <- as.numeric(as.character(year(sampling_date)))-1
      
      # Define the date from which we extract the data
      starting_date <- as.Date(paste(prev_year, winter_start, sep="-"),format="%Y-%m-%d")
      
      # Extraction of winter data (between starting_date et sampling_date)
      prev_data_winter <- env_point_print %>% 
        mutate_if(is.factor,as.character) %>% filter(site_name==sampling_point & date >= starting_date & date <= sampling_date)
      
      # Calculation of variables for winter and formating
      res_winter <- prev_data_winter %>% 
        # Selection of variables to integrate
        dplyr::select(varnames) %>% 
        # Integration by sum, quantile, sd ....
        summarise_all(funs(min,Q25=quantile(.,0.25), mean,Q75=quantile(.,0.75), max, sd)) %>%
        # Add a mention of "winter" in the variable name
        set_names(str_c(names(.),"w",sep="_")) %>%
        # Add sampling ID
        cbind(id_echant_print[i,],.)
      
      if(summer){
        
      # Define summer start and summer end
        summer_start <- as.Date(paste(prev_year,summer_date[1],sep="-"),format="%Y-%m-%d")
        summer_end <- as.Date(paste(prev_year,summer_date[2],sep="-"),format="%Y-%m-%d")
        
      # Extraction of summer data (between summer_start and summer_end)
        prev_data_summer <- hydro_daily_point %>% 
          mutate_if(is.factor,as.character) %>% filter(site_name==sampling_point & date >= summer_start & date <= summer_end)
        
      # Calculation of variables for summer and formating
        res_summer <- prev_data_summer %>% 
          # Selection of variables to integrate
          dplyr::select(varnames) %>% 
          # Integration by sum, quantile, sd ....
          summarise_all(funs(min,Q25=quantile(.,0.25), mean,Q75=quantile(.,0.75), max, sd)) %>%
          # Add a mention of "summer" in the variable name
          set_names(str_c(names(.),"s",sep="_"))
        
        res_sample <- cbind(res_winter,res_summer)  
      }else{
        res_sample <- res_winter
      }
      
      # Save results
      prev_out <- rbind(prev_out,res_sample)

    }
    
     return(prev_out)
    
  }
}

# Function apply on environmental data
integr_print <- integr(env_point_print,id_echant_print,varnames=c("bottomT", "current", "so", "mlotst","FF_Q", "PRETOTQ_Q", "T_Q", "TINF_H_Q", "TSUP_H_Q", "drange_T"),verbose=T)
integr_print$annee<-as.factor(integr_print$annee)
colnames(integr_print) <- str_remove(colnames(integr_print), "[_]w")
```

## Fetch data added

This variable is only spatial so not integrated

```{r}
fetch_site <- fetch_point %>%
  filter(theme == "Herbiers Intertidaux") %>%
  dplyr::select(-North,-South, -East, -West) %>%
  distinct(theme, site, point_name, Average) %>%
  dplyr::rename(fetch = Average)
```

## Environmental dataframe 

```{r}
env_print <- integr_print %>%
  left_join(fetch_site) %>%
  dplyr :: select(!ends_with(c("min", "Q25", "Q75", "max"))) %>% dplyr :: select(-site_name) %>%
  filter(!annee == 2020) %>%
  dplyr :: rename(courant_moy = current_mean, courant_sd = current_sd, vent_moy = FF_Q_mean, vent_sd = FF_Q_sd, temperature_eau_moy = bottomT_mean, temperature_eau_sd = bottomT_sd,temperature_air_moy = T_Q_mean, temperature_air_sd = T_Q_sd, precipitation_moy = PRETOTQ_Q_mean, precipitation_sd = PRETOTQ_Q_sd, salinite_moy=so_mean) %>%
  group_by(site,point_name,annee) %>% summarise_if(is.numeric,  mean) %>% ungroup()
```

## Separation of environmental matrix and metadata

```{r}
env_print_data <- env_print %>%
  as.data.frame() %>%
  set_rownames(paste(.$site, .$point_name, .$annee, sep="_")) %>%
  dplyr :: select( -site, -point_name,-annee)

env_print_meta <- env_print %>%
  dplyr :: select(site, point_name, annee)
```

## Temporal series 

```{r}
# Mean Air temperature (after integration)
ggplot(env_print, aes(x = annee, y = temperature_air_moy)) +
  labs(x = "Année", y = "Température de l'air moyenne (°C)") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"),
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

# Mean seawater temperature (after integration)
ggplot(env_print, aes(x = annee, y = temperature_eau_moy)) +
  labs(x = "Année", y = "Température de l'eau moyenne (°C)") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"),
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

# Mean Wind velocity (after integration)
ggplot(env_print, aes(x = annee, y = vent_moy)) +
  labs(x = "Année", y = "Vitesse du vent moyenne (m/s)") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"),
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

# Mean total precipation
ggplot(env_print, aes(x = annee, y = precipitation_moy)) +
  labs(x = "Année", y = "Précipitations totales moyennes (mm)") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"),
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

# Mean current velocitie (after integration)
ggplot(env_print, aes(x = annee, y = courant_moy)) +
  labs(x = "Année", y = "Vitesse du courant (m/s)") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"),
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))

# Mean salinity (after integration)
ggplot(env_print, aes(x = annee, y = salinite_moy)) +
  labs(x = "Année", y = "Salinité moyenne") +
  geom_point(aes(color = annee), size = 5) +
  facet_wrap(vars(site)) +
  theme(
    axis.text.x = element_text(size = 16, angle = 90),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 25, face = "bold"),
    strip.text = element_text(size = 20, face = "bold"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 16))
```

# Final environmental dataframe : granulo + hydrology + meteo

## Separation of environmental matrix and metadata

```{r}
granul_hydro_print <- env_print %>%
  inner_join(granul_noNA)

granul_hydro_print_data <- granul_hydro_print %>%
  as.data.frame() %>%
  set_rownames(paste(.$site, .$point_name, .$annee, sep = "_")) %>%
  dplyr :: select(-site, -point_name, -annee)

granul_hydro_print_meta <- granul_hydro_print %>%
  dplyr :: select(site, point_name, annee)
```

## PCA Hydrology_granulometry

```{r}
# PCA 
# ---

pca_granul_hydro <- granul_hydro_print_data %>%
  dplyr :: select (courant_moy, fetch, vent_moy, temperature_eau_moy, temperature_air_moy, temperature_eau_sd, temperature_air_sd, gravier, vase, trask, matiere_organique) %>%
  dplyr :: rename (current_mean = courant_moy, wind_mean = vent_moy, water_temp_mean = temperature_eau_moy, air_temp_mean = temperature_air_moy, water_temp_sd = temperature_eau_sd, air_temp_sd = temperature_air_sd, gravel = gravier, mud = vase, organic_matter = matiere_organique ) %>% 
  dudi.pca (scannf = FALSE, nf = 2)

pca_granul_hydro$li$Axis2 <- pca_granul_hydro$li$Axis2 * (-1)
pca_granul_hydro$co$Comp2 <- pca_granul_hydro$co$Comp2 * (-1)

pca_granul_hydro$li$Axis1 <- pca_granul_hydro$li$Axis1 * (-1)
pca_granul_hydro$co$Comp1 <- pca_granul_hydro$co$Comp1 * (-1)

# Get points coordinates
# ----------------------
point_scores_granul_hydro <- scores(pca_granul_hydro, display = "sites", scaling = 1) %>%
  as.data.frame() %>%
  rownames_to_column("location") %>%
  separate(location, c("site", "point_name", "annee"), sep = "_") %>%
  right_join(granul_hydro_print_meta, .)

# Get centroids of each point
# ---------------------------
point_scores_cent_granul_hydro <- point_scores_granul_hydro %>%
  group_by(site) %>% 
  summarise_at(vars(Axis1, Axis2), mean, na.rm = TRUE) %>%
  dplyr :: rename(cent1 = Axis1) %>%
  dplyr :: rename(cent2 = Axis2) %>%
  ungroup %>%
  left_join(point_scores_granul_hydro, .)

# Get centroids of each labels
#-----------------------------
label_granul_hydro <- point_scores_cent_granul_hydro %>%
  group_by(site) %>%
  summarise_at(vars(Axis1, Axis2),mean,na.rm = TRUE) %>%
  dplyr :: rename(cent1 = Axis1) %>%
  dplyr :: rename(cent2 = Axis2) %>%
  ungroup 

# Variance on each axis
#----------------------
var_axes_granul_hydro <- round(pca_granul_hydro$eig/sum(pca_granul_hydro$eig)*100,2)

# Plot of individuals
#--------------------
options(ggrepel.max.overlaps = Inf)
p <-ggplot(point_scores_granul_hydro, aes(x = Axis1, y = Axis2, col = site, fill = site, group = site))

p <- p + geom_point(aes(shape = site, color = site), size = 3, alpha = 0.7) + 
  scale_shape_manual(values=c(0,1,4,3,5,6,15,21,23)) + 
  scale_fill_manual(values = c("#25FDE9", "grey", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855")) +
  scale_color_manual(values = c("#25FDE9", "grey", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855")) +
  geom_point(data = point_scores_cent_granul_hydro, aes(x = cent1, y = cent2), size = 4, shape = 17, col = "black")

p <- p <- p + geom_label_repel(data = label_granul_hydro, aes(x = cent1, y = cent2,label = site), col = c("#25FDE9", "grey", "mediumpurple1", "#FF00FF", "#FF0000", "#1FA055", "#2C75FF", "#ED7F10", "#9FE855"), fill = "white", fontface = "bold",size = 6, segment.colour = "black", box.padding = 1, inherit.aes = FALSE, show.legend = FALSE, alpha = 0.8, nudge_x = 0.3, nudge_y = 0.3, label.size = 1)

p <- p + xlab(paste("PCA 1: 31.08 %")) + ylab(paste("PCA 2: 23.63 %"))

pca_ind_env_plot <- p + scale_y_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) +
scale_x_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) + 
  theme_bw() + theme(
    text=element_text(size=20), 
    strip.text.x = element_text(face = "bold"),
    legend.position="none",
    axis.title=element_text(size=23,face="bold"))
pca_ind_env_plot

# Plot of variables
#------------------

pca_var_env_plot <- fviz_pca_var(pca_granul_hydro, repel = TRUE, labelsize = 4) + 
  xlab("PCA 1: 31.08 %") + ylab("PCA 2: 23.63 %") +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) +
  scale_x_continuous(labels = function(x) format(x, scientific = FALSE, decimal.mark = '.')) + 
  ggtitle("") + theme_bw() + 
  theme(
  text = element_text(size = 15),
  legend.title = element_text(size = 0),
  legend.text = element_text(size = 6), 
  legend.position = "bottom",
  axis.title = element_text(size = 15,face = "bold"))
pca_var_env_plot 

```


## Correlation matrix 

```{r, results='asis'}

# Kendall

correl_granul_hydro <- granul_hydro_print_data %>%
  dplyr :: select(-boulder) %>%
  cor(.,method = "kendall", use = "pairwise.complete.obs")

corrplot(correl_granul_hydro, method = "number", type = "upper", order = "hclust", tl.col = "black", number.cex = 1.5, tl.cex = 2, cl.cex = 2, title = "Kendall's coefficient")

par(mfrow = c(1, 1))

# Pearson

correl_granul_hydro_pearson <- granul_hydro_print_data %>%
  dplyr :: select(-boulder) %>%
  cor(.,method = "pearson", use = "pairwise.complete.obs")

corrplot(correl_granul_hydro_pearson, method = "number", type = "upper", order = "hclust", tl.col = "black", number.cex = 1.5, tl.cex = 2, cl.cex = 2, title = "Pearson's coefficient")

par(mfrow = c(1, 1))
```

# Final dataframes to fit SEM

Dataframes are generated with different set of variables to fit SEM : for Infauna and Epifauna 
Some variables are log-transformed or standardised. 

## Endofaune

```{r}
# Species richness and Shannon
diversity_index_endo_fin <- diversity_index_endo %>%
  dplyr :: select(site, point_name, annee, rich_sp, shannon, abondance) %>%
  mutate(abondance = log1p(abondance)) %>%
  mutate(rich_sp = log1p(rich_sp)) %>%
  mutate(shannon = log1p(shannon)) %>%
  mutate_if(is.numeric, scale) %>%
  mutate_if(is.numeric,as.vector)

# Axis 1 : PCA infauna
point_scores_endofaune_1_fin <- point_scores_endofaune %>%
  dplyr :: select(site, point_name, annee, Axis1) %>%
  dplyr :: rename(axis1_endo = Axis1)

# Axis 2 : PCA infauna
point_scores_endofaune_2_fin <- point_scores_endofaune %>%
  dplyr :: select(site, point_name, annee, Axis2) %>%
  dplyr :: rename(axis2_endo = Axis2)

# Environment
granul_hydro_fin <- granul_hydro_print %>%
  mutate(matiere_organique = log(matiere_organique)) %>%
  mutate_if(is.numeric, scale) %>%
  mutate_if(is.numeric, as.vector)

# Biometry
biom_print_fin <- biom_noNA %>%
  dplyr :: select(-theme, -saison) %>%
  mutate(densite = log1p(densite)) %>%
  mutate(br.m2 = log1p(br.m2)) %>%
  mutate(bf.m2 = log1p(bf.m2)) %>%
  mutate(ratio.epi = sqrt(ratio.epi)) %>%
  mutate(biom.macro.moy = sqrt(biom.macro.moy)) %>%
  mutate_if(is.numeric, scale) %>%
  mutate_if(is.numeric,as.vector)

# Final dataframe
tab_final_endo <- diversity_index_endo_fin %>%
  inner_join(point_scores_endofaune_1_fin) %>%
  inner_join(point_scores_endofaune_2_fin) %>%
  inner_join(biom_print_fin) %>%
  inner_join(granul_hydro_fin) %>%
  dplyr:: select(-n.f.pied,-pourc.brok, -long.g.moy, -boulder)

# Correlation matrix : Pearson
#-----------------------------
correl_endo_fin <- tab_final_endo %>%
  dplyr :: select ( -site, -point_name, -annee) %>% 
  cor(.,method = "pearson", use = "pairwise.complete.obs")

corrplot(correl_endo_fin, method = "number", type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

# Correlation matrix : Kendall
#-----------------------------
correl_endo_fin <- tab_final_endo %>%
  dplyr :: select ( -site, -point_name, -annee) %>% 
  cor(.,method = "kendall", use = "pairwise.complete.obs")
corrplot(correl_endo_fin,method = "number", type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
```

## Epifauna

```{r}
# Species richness, Shannon et abundance 
diversity_index_epi_fin <- diversity_index_epi %>%
  dplyr :: select(site, point_name, annee, shannon, rich_sp, abondance) %>%
  mutate(abondance = log1p(abondance)) %>%
  mutate(rich_sp = log1p(rich_sp)) %>%
  mutate(shannon = log1p(shannon)) %>%
  mutate_if(is.numeric, scale) %>%
  mutate_if(is.vector, as.vector)

# Abundance grazers
abond_brouteur <- abond_brouteur_epi %>%
  dplyr :: select (site, point_name, annee, abond_brouteur) %>%
  mutate(annee = as.factor(annee)) %>%
  mutate(abond_brouteur = scale(log1p(abond_brouteur))) %>%
  mutate_if(is.numeric, as.vector)
         
# Axis 1 : PCA epifauna
point_scores_epifaune_1_fin <- point_scores_epifaune %>%
  dplyr :: select(site, point_name, annee, Axis1) %>%
  dplyr::rename(axis1_epi = Axis1)

# Axis 2 : PCA epifaune
point_scores_epifaune_2_fin <- point_scores_epifaune %>%
  dplyr :: select(site, point_name, annee, Axis2) %>%
  dplyr :: rename(axis2_epi = Axis2)

# Environment
granul_hydro_fin <- granul_hydro_print %>%
  mutate(matiere_organique = log(matiere_organique)) %>%
  mutate_if(is.numeric, scale) %>%
  mutate_if(is.numeric, as.vector)

# Biometry
biom_print_fin <- biom_noNA %>%
  dplyr :: select(-theme,-saison) %>%
  mutate(densite = log1p(densite)) %>%
  mutate(br.m2 = log1p(br.m2)) %>%
  mutate(bf.m2 = log1p(bf.m2)) %>%
  mutate(ratio.epi = sqrt(ratio.epi)) %>%
  mutate(biom.macro.moy = sqrt(biom.macro.moy)) %>%
  mutate_if(is.numeric, scale) %>%
  mutate_if(is.numeric,as.vector)

# Final dataframe
tab_final_epi <- diversity_index_epi_fin %>%
  inner_join(abond_brouteur) %>%
  inner_join(point_scores_epifaune_1_fin) %>%
  inner_join(point_scores_epifaune_2_fin) %>%
  inner_join(biom_print_fin) %>%
  inner_join(granul_hydro_fin) %>%
  dplyr :: select(-n.f.pied,-pourc.brok, -long.g.moy, -boulder)

# Correlation matrix : Pearson
#-----------------------------
correl_epi_fin <- tab_final_epi %>%
  dplyr :: select (-site, -point_name, -annee) %>% 
  cor(.,method = "pearson", use = "pairwise.complete.obs")

corrplot(correl_epi_fin,method = "number", type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

# Correlation matrix : Kendall
#-----------------------------
correl_epi_fin <- tab_final_epi %>%
  dplyr :: select ( -site, -point_name, -annee) %>% 
  cor(.,method = "kendall", use = "pairwise.complete.obs")

corrplot(correl_epi_fin, method = "number", type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
```

