The script of the function spat_temp_index can be found in the following repository and supplementary material

Link to the functions spat_temp_index: https://github.com/Cunillera-Montcusi/Quantifyinig-SpaTem-connectivity Link to supplementary material:

# Load functions to run the all the simulations
setwd("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/1. FEHM coses al DIA/4. Mecodispers Spa-Tem/MECODISPER SpaTem")
source("SpaTemp_HOBOS_function.R")

1. Charging & depurating HOBOS Dataset

The initial step is to charge the datasets containing the names of each monitored site containing the stream name, the ID code and the latitude and longitude. This dataset will be used to localize each one of the data loggers along the river (the “Longlat_Rius.csv” can be found in the repository - LINK)

The following step is to upload the HOBOS_sites information, which corresponds to the Water permanence database (Figure 1). We upload individually the tables for each stream (one sheet for each stream). These “.csv” documents can be found in the online repository - LINK and contain the information about each daily wet/dry status of each stream.

IMPORTANT INFORMATION: For each one of these “.csv” documents the data loggers must be ordered in Upstream to Downstream direction and must have the following structure:

data.frame(
MonitoredDays= c("Day1","Day2","Day3","Day4"), 
UpstreamSite1=c(0,1,1,1), 
StreamSite2=c(0,1,1,1),StreamSite3=c(1,1,1,1),StreamSite4=c(1,1,1,1),
DownstreamSite5=c(1,1,1,1) 
)

Finally, we upload the “Stream_order” database. This database contains the Upstream - Downstream “order” of each datalogger (1 being Upstream). The file “Stream Order.csv” can be found in the online repository - LINK. Later, we check and correct some coordinates to make sure that each datalogger has a unique coordinate and plot them individually to assess if they are correctly placed in space.

With this initial process we will depurate and charge the two most important informations required for the function: 1. Inermitence_dataset 2. Sites_coordinates

Sites <- read.csv("Raw_HOBOS_Database/Longlat_Rius.csv", header = T, sep = ";", dec = ",")
colnames(Sites) <- c("Riera", "Codi_HOBO","Latitud","Longitud")

#Correction for matching HOBOS -- collecting the  coordinates of the HOBOS that we have  
length(which(as.matrix(dist(Sites[,3:4]))==0))
length(diag(as.matrix(dist(Sites[,3:4]))))
plot(Sites$Longitud, Sites$Latitud)

# We upload the HOBO dataset (1 and 0 defining wet/dry moments)
# We upload the order of the rivers from Upstream to Downstream to arrange the HOBOS into the desired order
## IMPORTANT: This order will also define the direction of the river! 
## The first HOBO in the row of the "Sites_list" must be the first HOBO in the Column of "HOBOS_sites"

HOBOS_sites <- list(
  read.csv("Raw_HOBOS_Database/CA_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/M_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/R_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/SA_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/SC_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/T_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/VH_HOBOS_data.csv", header = T, sep = ";"))

#Stream order to check
Stream_order <- read.csv("Raw_HOBOS_Database/Stream Order.csv", header = T, sep = ";")

library(dplyr)
# PLotting HOBOS altogether (make the plot window bigger)
par(mfrow=c(3,3))
local_hobos_list <- list()
Sites_list <- list()
for (sit in 1:length(HOBOS_sites)) {
  colnames(HOBOS_sites[[sit]])[1] <- c("Day")
  names_hobos <- colnames(HOBOS_sites[[sit]])[2:length(colnames(HOBOS_sites[[sit]]))]
  local_hobos <- c()
  Sites_list[[sit]] <- Sites%>%filter(Codi_HOBO%in%names_hobos)%>%
                               left_join(Stream_order, by="Codi_HOBO")%>%
                               arrange(UtoD)%>%
                               dplyr::select(Riera,Codi_HOBO,Longitud,Latitud)
#plot(Sites_list[[sit]]$Longitud, Sites_list[[sit]]$Latitud)
}
par(mfrow=c(1,1))

2. Dates selection

To select different “time windows” in order to calculate ST indicies for the the whole monitored period (513 days) or until the Autumn sampling (18-21/11/2019; 480 days) we upload the biological information files and retrieve and modify the dates in order to easily select them and calculate the difference between the beginning (“date_HOBOS” = 2018-07-26) and the other two dates, the end (“date_HOBOS_end”=2019-12-20) and the autumn sampling (“date_correct_Autumn”).

As each row of the matrix “HOBOS_sites” corresponds to a unique day, the difference (in days) between two dates (beginning vs end or beginning vs Autumn sampling) will deliver the number of rows that we have to select from the “HOBOS_sites” in order to account for the desired “time-window”. Therefore, after this second step we will have the correct HOBOS_sites matrix corresponding to the time-window that we want to focus.

### Autumn dataset - November 2019
BDD <- read.csv2("BiolData/Matriz_otoño.csv", sep=";")
date_correct_Autumn <- min(unique(BDD$Date))
# Edit the dates according the desired format to be entered in the following steps
date_correct_Autumn <- strsplit(date_correct_Autumn,"/") # Split the numbers by the "/"
date_correct_Autumn <- paste(date_correct_Autumn[[1]][3],"-", # Paste the numbers in the correct order and separated by "-" 
                             date_correct_Autumn[[1]][2],"-",
                             date_correct_Autumn[[1]][1],sep = "")

date_correct_Autumn
## All hobos begun at the same date
date_HOBOS <- "2018-07-26"

## All hobos finish at the same date
date_HOBOS_end <- "2019-12-20"

#Beginning
bd <- as.Date(date_HOBOS)

#End
# We change the end for the total dataset (513 days) or for comparing only with Autumn
ed <- as.Date(date_correct_Autumn)

# Difference in number of days will correspond to the number of rows to be selected. 
time_window_beg <- as.numeric(difftime(bd+1, date_HOBOS, units = "days")) 
time_window_beg <- round(time_window_beg)
time_window_end <- as.numeric(difftime(ed, date_HOBOS, units = "days"))
time_window_end <- round(time_window_end)

for (site in 1:length(HOBOS_sites)) {
  HOBOS_sites[[site]] <- HOBOS_sites[[site]][time_window_beg:time_window_end,]
}

3. Distance matrix creation

If we wish to calculate weighted indices we have to create a weighted pairwise matrix with the desired distance. In our case we caluclate the euclidean distances based on the longitude and latitude. We use the coordinates information contained in “Sites_list” object.

By doing this we obtain another required information for running the function, the dist_matrices, here named “eucl_dist_matrices” as we used the direct euclidean distances between monitored sites.

# Must be a list of distances matrix for each river. 
# Must be a symmetric matrix (upper and lower triangles must be equal). The function already selects which distances 
#are selected in each case.  
# Distances must be small or the function my increase a lot computational demands! In this case they are in km
eucl_dist_matrices <- list()
for (river in 1:length(Sites_list)) {
  eucl_dist_matrices[[river]] <-as.matrix(dist(Sites_list[[river]][,4:3],diag = T,upper = T))/1000 
}

4. Network structure creation

Finally, we need to set the network structure. This is a representation of how our sampled sites are organized in space. We do this by creating an adjacency matrix that corresponds to a head to downstream direction or directed structure and anothe where all sites can be reached from everywhere or undirected structure. By doing this we finalize all the required information to run the function.

# Must be a list with the network strucutre for each river.
# For directed scenarios
Dir_netw_struct_matrices <- list()
for (river in 1:length(HOBOS_sites)) {
  Dir_SpatConnect <- seq(1:c(ncol(HOBOS_sites[[river]])-1))
  netw_struct_matrix <- matrix(nrow =length(Dir_SpatConnect),ncol = length(Dir_SpatConnect),0)
  for (site_step in 1:length(Dir_SpatConnect)-1) {
    netw_struct_matrix[Dir_SpatConnect[site_step],Dir_SpatConnect[site_step]+1] <- 1}
Dir_netw_struct_matrices[[river]] <- netw_struct_matrix
}

# For undirected scenarios 
UNdir_netw_struct_matrices <- list()
for (river in 1:length(HOBOS_sites)) {
  Dir_SpatConnect <- seq(1:c(ncol(HOBOS_sites[[river]])-1))
  netw_struct_matrix <- matrix(nrow =length(Dir_SpatConnect),ncol = length(Dir_SpatConnect),0)
  for (site_step in 1:length(Dir_SpatConnect)) {
    netw_struct_matrix[Dir_SpatConnect[site_step],
                       c(Dir_SpatConnect[1]:length(Dir_SpatConnect))[-site_step]]<- 1}
UNdir_netw_struct_matrices[[river]] <- netw_struct_matrix
}

Function spat_temp_index elements to consider before running it:

LINK and NO-LINK values: The meaning of a “LINK” (Input matrix value= 1) or a “NO-LINK” (Input matrix value= 0) can be modified inside the function according to the different interpretations that one once to give to the values (e.g. do we want to quantify connectivity? dispersal resistance?).

Direction:direction can be either “directed” or “undirected”. This feature will modify the way the graph is being controlled.

Sense: this is referring to the direction that will be considered for the graph connectivity considered. Specially for centrality metrics. It can be either “in”, “out”, “all”. See ?igraph or ?igraph::closeness for a better understanding. In summary, when considering “out” sense we are only taking into account the “connection leaving” the site (source effect), when considering “in” sense we consider the connectivity entering each site (sink effect), and with “all” sense we consider everyting. The direction and sense can be combined in order to precisely define how connectivity is being understood.

Weighting: weighting either FALSE or TRUE. The value of the weight is defined by the matrix added (must be in a distance matrix format) and can consider any type of distance between pairs of monitored sites (euclidean, environmental, topographic, …). The dist_matrices (alreday created in this example, see section 3) is mandatory when WEIGHTING=TRUE.

Network_structure: Already created (see section 4). It is a matrix that corresponds to the “basic” connections that can be possible. An adjacency matrix where all sites are connected among them. Should be equivalent to the network that you would expect from a “fully” connected network.

weighting_links=FALSE & link_weights: If instead of using the values used as “LINKS/NO-LINKS” we want to use specific data wuantifying for every time unit the connections “strength” between sites (e.g. flow data, wind strength, etcetera) we must activate these two elements. This information can be incorporated by selecting weighting_links=TRUE and then providing a list of data.frames with exactly the same information as the Intermitence_dataset, having the same amount of rows and columns (e.g. daily flow data in each site).

LINKS / NO-LINKS: The value_LINK is the value attributed for each “effective link”, which is a connection between two nodes (a line). The value_S_LINK for spatial links and the value_T_LINK for temporal links. The value_NO_link is the value for each “effective disconnection”, which is a “connection void” between two sites (a white space). The value_NO_S_LINK for spatial links and the value_NO_T_LINK for temporal links.

Legacy effects & Legacy lenght:Legacy effects are a way to quantify spatiotemporal connectivity during several time units all together. This means that we can quantify temporal links for more than 1 time unit and modify the weight or relevance that each time unit has. For example: - Legacy lenght defines the “number” of time units considered. By default it is set to 1 which means that the function will consider only 1 time unit (the connectivity of time T1 will only be considered for T2). If value is set at legacy lenght=3 the connectivity of T1 will be considered for T2, T3, and T4 (the connections present at T1 will be also considered for the 3 following T). - Legacy effect defines the weight given to the LINK / NO-LINK values for each time unit. It is a vector with values ranging from 0 to 1 that will modulate the relevance of a connection through time.

Ploting and other elements of the function: - Virid_option defines which palette the user wants to use for the plots (“A”,“B”,“C”,…) from the viridis package (see ?viridis for more information). If not specified the default is set to the “viridis” colour palette
- Network_variables=F activates the calculation of centrality metrics for each network for each analysed day. These metrics include closennes, outclosennes and betweenness. They will be included in the output as a data.frame with the values for each site (colums) for each day (rows). - print.plots=F defines if we want to directly print the plots or not in a directory (FOR BIG NETWORKS with hundreds of nodes it can be computationally demanding and it can even block the session). - print.directory=“Figure/” sets the directory where plots will be saved.

Finally, we define all these items for each of the four scenarios:

1. Scenario DirBin - Binnary directed stream: Scenario considering a directed network (upstream-downstream direction; direction=“directed”) and ponderating the “links” as “1”( value_S_LINK=1 & value_T_LINK=1 ) and the “no links” as “0”( value_NO_S_link=0 & value_NO_T_link=0 ). Therefore, evaluating the connectivity ( weighting=FALSE ) of each monitored site according to water absence/presence following stream flow direction.

2. Scenario DirWei - Weighted directed stream: Scenario considering a directed network (upstream-downstream direction; direction=“directed”), ponderating the “links” as “0.1”( value_S_LINK=0.1 & value_T_LINK=0.1 ) and the “no links” as “1”( value_NO_S_link=1 & value_NO_T_link=1 ) and weighting them according to euclidean distance between each pair of sites ( weighting=TRUE & dist_matrices = eucl_dist_matrices ). Therefore, evaluating the dispersal resistance of each monitored site according to water absence/presence and distance between sites (when dry, sites are considered further) following stream flow direction.

3. Scenario UndBin - Binnary undirected stream: Scenario considering an undirected network (all nodes can be reached; direction=“undirected”) and ponderating the “links” as “1”( value_S_LINK=1 & value_T_LINK=1 ) and the “no links” as “0”( value_NO_S_link=0 & value_NO_T_link=0 ). Therefore, evaluating the connectivity ( weighting=FALSE ) of each monitored site according to water absence/presence without considering any direction.

4. Scenario UndWei - Weighted undirected stream: Scenario considering an undirected network (all nodes can be reached; direction=“undirected”), ponderating the “links” as “0.1”( value_S_LINK=0.1 & value_T_LINK=0.1 ) and the “no links” as “1”( value_NO_S_link=1 & value_NO_T_link=1 ) and weighting them according to euclidean distance between each pair of sites ( weighting=TRUE & dist_matrices = eucl_dist_matrices ). Therefore, evaluating the dispersal resistance of each monitored site according to water absence/presence and distance between sites (when dry, sites are considered further) without considering any direction.

Scenario DirBin - Binnary directed stream

# Building a directed non weighted network
# Value of link=1
# Value of NO link=0
Dir_NonW_Net <- spat_temp_index(Inermitence_dataset = HOBOS_sites,
                                Sites_coordinates=Sites_list,
                                direction="directed", sense = "out",
                                Network_stru =Dir_netw_struct_matrices,
                                weighting=FALSE,dist_matrices = NULL,
                                weighting_links =FALSE,link_weights = NULL,
                                legacy_effect =1, legacy_lenght = 1,
                                value_S_LINK=1,
                                value_T_LINK=1,
                                value_NO_S_link=0,
                                value_NO_T_link=0,
                                Network_variables=F,print.plots=F,print.directory="Figure/")

Scenario DirWei - Weighted directed stream

# Building a directed Weighted network
# Value of link= 0.1
# Value of NO link=1
# - I am evaluationg "resistance" to move from one to another. 
# SO, higher numbers mean High resistance
Dir_WEIG_Net <- spat_temp_index(Inermitence_dataset =HOBOS_sites,
                                Sites_coordinates=Sites_list,
                                direction="directed", sense = "out",
                                weighting=TRUE, dist_matrices = eucl_dist_matrices,
                                Network_stru =Dir_netw_struct_matrices,
                                weighting_links =FALSE,link_weights = NULL,
                                legacy_effect =1, legacy_lenght = 1,
                                value_S_LINK=0.1,
                                value_T_LINK=0.1,
                                value_NO_S_link=1,
                                value_NO_T_link=1,Virid_option = "B",
                                Network_variables=F,print.plots=F,print.directory="Figure/")

Scenario UndBin - Binnary undirected stream

# Building a Undirected non weighted network
# Value of link=1
# Value of NO link=0
UnD_NonW_Net <- spat_temp_index(Inermitence_dataset =HOBOS_sites,
                                Sites_coordinates=Sites_list,
                                direction="undirected", sense = "all",
                                weighting=FALSE,dist_matrices = NULL,
                                Network_stru =UNdir_netw_struct_matrices, 
                                weighting_links =FALSE,link_weights = NULL,
                                legacy_effect =1, legacy_lenght = 1,
                                value_S_LINK=1,
                                value_T_LINK=1,
                                value_NO_S_link=0,
                                value_NO_T_link=0,
                                Network_variables=F,print.plots=F,print.directory="Figure/")

Scenario UndWei - Weighted undirected stream

# Building a Undirected Weighted network
# Value of link= 0.1
# Value of NO link=1
# - I am evaluationg "resistance" to move from one to another. 
# SO, higher numbers mean High resistance
UnD_WEIG_Net <- spat_temp_index(Inermitence_dataset = HOBOS_sites,
                                Sites_coordinates=Sites_list,
                                direction="undirected",  sense = "all",
                                weighting=TRUE, dist_matrices = eucl_dist_matrices,
                                Network_stru =UNdir_netw_struct_matrices,
                                weighting_links =FALSE,link_weights = NULL,
                                legacy_effect =1, legacy_lenght = 1,
                                value_S_LINK=0.1,
                                value_T_LINK=0.1,
                                value_NO_S_link=1,
                                value_NO_T_link=1,Virid_option = "B",
                                Network_variables=F,print.plots=F,print.directory="Figure/")

DirBin Scenario - Binnary directed stream

# STcon 
NonW_ST_connectivity_value <- Dir_NonW_Net$STcon
NonW_ST_directed_output <- data.frame(NonW_Dir_con=unlist(NonW_ST_connectivity_value))
# STconmat 
NonW_ST_matrix_out_out <- Dir_NonW_Net$STconmat

DirWei Scenario - Weighted directed stream

# STcon
WEIG_ST_connectivity_value <- Dir_WEIG_Net$STcon
WEIG_ST_directed_output <- data.frame(WEIG_Dir_con=unlist(WEIG_ST_connectivity_value))
# STconmat 
WEIG_ST_matrix_out_out <- Dir_WEIG_Net$STconmat

UndBin Scenario- Binnary undirected stream

# STcon 
Un_NonW_ST_connectivity_value <- UnD_NonW_Net$STcon
Un_NonW_ST_output <- data.frame(Un_NonW_con=unlist(Un_NonW_ST_connectivity_value))
# STconmat
Un_NonW_ST_matrix_out_out <- UnD_NonW_Net$STconmat

UndWei Scenario- Weighted undirected stream

# STcon
Un_WEIG_ST_connectivity_value <- UnD_WEIG_Net$STcon
Un_WEIG_ST_output <- data.frame(Un_WEIG_con=unlist(Un_WEIG_ST_connectivity_value))
# STconmat 
Un_WEIG_ST_matrix_out_out <- UnD_WEIG_Net$STconmat
                                Un_WEIG_Bet=unlist(Un_WEIG_ST_Betclo_mean))

# We add the stream ID and the order (UtoD) to each output to later analyse them. 
# We retrieve this information from the "Stream_order" dataset where each HOBO code, Stram and its position on the 
#upstream downstream direction were registered
# This two following vectors are created for being used in other scripts
HOB_riv_ID<- Stream_order$ï..Riera
ups_dos <- Stream_order$UtoD

NonW_ST_directed_out <- data.frame(ID=Stream_order$ï..Riera,DtoU=Stream_order$UtoD,
                                   NonW_ST_directed_output)%>%
                        mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))

WEIG_ST_directed_out <- data.frame(ID=Stream_order$ï..Riera,DtoU=Stream_order$UtoD,
                                   WEIG_ST_directed_output)%>%
                        mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))

Un_NonW_ST_out <- data.frame(ID=Stream_order$ï..Riera,DtoU=Stream_order$UtoD,
                             Un_NonW_ST_output)%>%
                        mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))

Un_WEIG_ST_out <- data.frame(ID=Stream_order$ï..Riera,DtoU=Stream_order$UtoD,
                             Un_WEIG_ST_output)%>%
                        mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))


# Final values HOBOS dataframe
NonW_ST_matrix_out_out
WEIG_ST_matrix_out_out
Un_NonW_ST_matrix_out_out
Un_WEIG_ST_matrix_out_out

NonW_ST_directed_out
WEIG_ST_directed_out
Un_NonW_ST_out
Un_WEIG_ST_out

Initial and required data to proceed (previously crated)

# Cunillera_palette - Available at https://github.com/Cunillera-Montcusi/Cunillera_Pallete
source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/CUNILLERA_palette.R")
setwd("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/1. FEHM coses al DIA/4. Mecodispers Spa-Tem/MECODISPER SpaTem")

## IMPORTANT MESSAGE: YOU NEED TO RUN the treatment script
# Needed values from the script "SpaTemp_HOBOS_treatment.R"
## RUN BEFORE PROCEEDING WITH THIS ONE

## You need to run the previous script in order to built the HOBOS dataset (a list with the HOBOS info for each river)
HOBOS_sites
# This two other objects are also built in the "SpaTemp_HOBOS_treatment.R" script. They contain the names 
## of each river and the order upstream to downstream within it. 
HOB_riv_ID
ups_dos

# Matrix ST outputs obtained from "SpaTemp_HOBOS_treatment.R"
NonW_ST_matrix_out_out
WEIG_ST_matrix_out_out
Un_NonW_ST_matrix_out_out
Un_WEIG_ST_matrix_out_out

# INDIVIDUAL ST outputs obtained from "SpaTemp_HOBOS_treatment.R"
NonW_ST_directed_out
WEIG_ST_directed_out
Un_NonW_ST_out
Un_WEIG_ST_out

Individual HOBO values and comparisons with old metrics - STcon tratment

We calculate a different set of metrics that where previously used for this type of dataset to evaluate drought patterns. For each site we calculated the total duration of drying events (i.e. the total number of dry days; hereafter TotDur), the frequency of drying events (i.e. the number of changes from wet to dry days; hereafter TotNum), and the average length of each drying event (the average number of dry days for each drying event; hereafter TotLeng).

The script to calculate these metrics can be found in the “Old_HOBOS_calculation.R” file found in the following repository - LINK

source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/1. FEHM coses al DIA/4. Mecodispers Spa-Tem/MECODISPER SpaTem/FlowIntermittence_indices.R")
Old_HOBOS_comb

Once obtained, we compare them with the STcon indices obtained for the four scenarios studied (PCA analysis-Figure 3) and analyse the correlation among them (Figure S3 ). More information in the article DOI:

Figure 3: PCA plot for STcon indices for each one of the four scenarios: Directed binary (DirBin), Directed Weighted (DirWei), Undirected Binary (UndBin), Undirected Weighted (UndWei) and the previously existing indices for water permanence data: the total duration of drying events (TotDur), the frequency of drying events (TotNum), and the average length of each drying event (TotLeng). All the data is calculated based on the total monitored time window monitored (513 days).

# Figure 3
data_PCA <- data.frame(data.frame("DirBin"=NonW_ST_directed_out$NonW_Dir_con,
                                  "DirWei"=WEIG_ST_directed_out$WEIG_Dir_con,
                                  "UndBin"=Un_NonW_ST_out$Un_NonW_con, 
                                  "UndWei"=Un_WEIG_ST_out$Un_WEIG_con),
            Old_HOBOS_comb[,3:5])


png(filename =paste("Figure/Data_treat/","Figure3",".png"), 
    width = 700*2.5, height = 550*2.5, 
    units = "px",res = 300)

PCA_info <- prcomp(data_PCA,center = T, scale. = T)
a <- data.frame(
     adonis2(PCA_info$x[,2]~NonW_ST_directed_out$ID*NonW_ST_directed_out$DtoU,
        method = "euclidean")
      )
rownames(a) <- c("Stream ID","Upstr-Downstr", "Interaction","Residual","Total")
caption_text <- paste(rownames(a)[1],"pvalue=",a[1,5]," ",
                       rownames(a)[2],"pvalue=",a[2,5]," ",
                       rownames(a)[3],"pvalue=",a[3,5], sep = " ")
autoplot(prcomp(data_PCA,
  center = T, scale. = T),
  loadings=T, loadings.label = TRUE, loadings.label.size = 6, shape = F, label=F,
  loadings.colour = 'grey15', loadings.label.colour="black")+
  geom_point(alpha=0.5,shape=21,
             color="grey30", 
             aes(size=NonW_ST_directed_out$DtoU ,fill=NonW_ST_directed_out$ID))+
  stat_ellipse(aes(colour=NonW_ST_directed_out$ID),type = "t", linetype=2,size=1, alpha=0.5)+
  scale_fill_CUNILLERA(palette = "LGTBI", name="Stream ID")+
  scale_color_CUNILLERA(palette = "LGTBI", name="Stream ID")+
  scale_size(name="Upstr-Downstr")+
  scale_x_continuous(limits = c(-0.35,0.4))+
  scale_y_continuous(limits = c(-0.32,0.32))+
  labs(caption = caption_text)+
  theme_bw()
dev.off()

Figure S5: Correlation plot among all the indices obtained for the four scenarios and the previously used drought-related metrics.

# Figure S5
corrmorant::corrmorant(cbind("DirBin"=NonW_ST_directed_out$NonW_Dir_con,
                             "log(DirWei)"=log(WEIG_ST_directed_out$WEIG_Dir_con+1),
                             "STcon UndBin"=Un_NonW_ST_out$Un_NonW_con,
                             "log(UndWei)"=log(Un_WEIG_ST_out$Un_WEIG_con+1),
                             "log(TotDur)"= log(Old_HOBOS_comb[,3]+1),
                             "log(TotNum)"= log(Old_HOBOS_comb[,4]+1),
                             "log(TotLeng)"= log(Old_HOBOS_comb[,5]+1)),
                            rescale = "by_sd",corr_method = "spearman")

Matrix HOBOS values NMDS building - STconmat and stream characterization using both indices

For each scenario (DirBin, DirWei, UndBin and UndWei) we calculate an NMDS using euclidean distances and ploting their values in order to characterize each one of the analysed streams in a bidimensional space.

The obtained indices can be used to characterise the river spatiotemporally by ploting their STcon vlaues along the stream direction for example, wich allows to detect more disconnected sites or the impact of determined sites on the rest of the network (e.g. directed networks).

Finally, we merge the obtained results in a unique figure where the two developed indices are shown. Figure 4 of the article (DOI: LINK) provides a complete image of streams spatiotemporal patterns based on both indices.

png(filename =paste("Figure/Data_treat/Figure4.png"), 
      width = 700*6, height = 650*6, 
      units = "px",res = 350)
  rows_filter <- unlist(c(four_approxim[[1]]%>%
                            mutate(files=1:nrow(.))%>%
                            group_by(ID)%>%
                            mutate(maxim=max(DtoU), value=DtoU/maxim)%>%
                            filter(value==1)%>%
                            ungroup()%>%
                            dplyr::select(files)))
  
  alpha_stream <- as.numeric(unlist(c(four_approxim[[1]]%>%
                  mutate(alpha_stream=ifelse(ID=="VH", 0.8, 0.8))%>%
                  ungroup()%>%
                  dplyr::select(alpha_stream))))
  alpha_HOBOS3 <- c(rep(0,58),1,rep(0,7))
  
  grid.arrange(
  #A1 line plot  
  arrangeGrob(
    ggplot(data = four_approxim[[1]][-rows_filter,])+
      #geom_point(aes(x=DtoU,y=four_approxim[[1]][-rows_filter,3]),color="red",alpha=alpha_HOBOS3[-rows_filter], shape=16, size=5,)+
      geom_point(aes(x=DtoU,y=four_approxim[[1]][-rows_filter,3],fill=ID),color="grey30", alpha=alpha_stream[-rows_filter], shape=21, size=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[1]][-rows_filter,3],color=ID),alpha=alpha_stream[-rows_filter], linetype=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[1]][-rows_filter,3],color=ID, alpha=ID), 
                size=1.3,stat="smooth",method = "lm")+
      scale_alpha_manual(values = c(rep(0.3,6),0.9))+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "DirBin - STcon", title = "A)")+
      xlab("HOBO relative position Upstream-Downstream")+
      ylab("ST Connecticity")+
      theme_classic()+theme(legend.position="none"),
    top = ""),
  
  #A1 NMDS plot  
  arrangeGrob(    
    ggplot(NonW_ST_directed_MatrixOut)+
      #geom_point(aes(x=x,y=y, fill=ID,size=DtoU),color="red",alpha=alpha_HOBOS3, shape=16, size=5,)+
      geom_point(aes(x=x, y=y, fill=ID,size=DtoU), shape=21,colour="grey10", alpha=alpha_stream)+
      stat_ellipse(aes(x=x, y=y, colour=ID),type = "t", linetype=2,size=1, alpha=0.5)+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "DirBin - STconmat", title = "B)")+
      scale_y_continuous(label=NULL)+
      scale_x_continuous(label=NULL)+
      theme_classic()+theme(legend.position="none",axis.ticks = element_blank())+
      facet_grid(. ~ ID),
    top = ""),
    
  
  #A2 line plot  
  arrangeGrob( 
    ggplot(data = four_approxim[[2]][,])+
      #geom_point(aes(x=DtoU,y=four_approxim[[2]][,3]),color="red",alpha=alpha_HOBOS3, shape=16, size=5,)+
      geom_point(aes(x=DtoU,y=four_approxim[[2]][,3],fill=ID),color="grey30", alpha=alpha_stream, shape=21, size=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[2]][,3],color=ID),alpha=alpha_stream, linetype=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[2]][,3],color=ID, alpha=ID), 
                size=1.3,stat="smooth",method = "lm")+
      scale_alpha_manual(values = c(rep(0.3,6),0.9))+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "DirWei - STcon")+ 
      xlab("HOBO relative position Upstream-Downstream")+
      ylab("ST Connecticity")+
      scale_y_continuous(limits = c(0,15000))+
      theme_classic()+theme(legend.position="none"),
    top=""),
    
  #A2 NMDS plot
  arrangeGrob(
    ggplot(WEIG_ST_MatrixOut)+
      geom_point(aes(x=x, y=y, fill=ID,size=DtoU), shape=21,colour="grey10", alpha=alpha_stream)+
      stat_ellipse(aes(x=x, y=y, colour=ID),type = "t", linetype=2,size=1, alpha=0.5)+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      scale_y_continuous(label=NULL)+
      scale_x_continuous(label=NULL)+
      labs(subtitle = "DirWei - STconmat")+
      theme_classic()+theme(legend.position="none",axis.ticks = element_blank())+
    facet_grid(. ~ ID), top = ""),
  
  #B1 line plot
  arrangeGrob(
    ggplot(data = four_approxim[[3]][,])+
      #geom_point(aes(x=DtoU,y=four_approxim[[3]][,3]),color="red",alpha=alpha_HOBOS3, shape=16, size=5,)+
      geom_point(aes(x=DtoU,y=four_approxim[[3]][,3],fill=ID),color="grey30", alpha=alpha_stream, shape=21, size=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[3]][,3],color=ID),alpha=alpha_stream, linetype=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[3]][,3],color=ID, alpha=ID), 
                size=1.3,stat="smooth",method = "lm")+
      scale_alpha_manual(values = c(rep(0.3,6),0.9))+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "UndBin - STcon")+
      xlab("HOBO relative position Upstream-Downstream")+
      ylab("ST Connecticity")+
      scale_y_continuous(limits = c(0,11.5))+
      theme_classic()+theme(legend.position="none"),
    top=""),
  
  #B1 NMDS plot
  arrangeGrob(
    ggplot(Un_NonW_ST_MatrixOut)+
      geom_point(aes(x=x, y=y, fill=ID,size=DtoU), shape=21,colour="grey10", alpha=alpha_stream)+
      stat_ellipse(aes(x=x, y=y, colour=ID),type = "t", linetype=2,size=1, alpha=0.5)+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "UndBin - STconmat")+
      scale_y_continuous(label=NULL)+
      scale_x_continuous(label=NULL)+
      theme_classic()+theme(legend.position="none",axis.ticks = element_blank())+
      facet_grid(. ~ ID),top = ""),
  
  #B2 line plot  
  arrangeGrob(
    ggplot(data = four_approxim[[4]][,])+
      #geom_point(aes(x=DtoU,y=four_approxim[[4]][,3]),color="red",alpha=alpha_HOBOS3, shape=16, size=5,)+
      geom_point(aes(x=DtoU,y=four_approxim[[4]][,3],fill=ID),color="grey30", alpha=alpha_stream, shape=21, size=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[4]][,3],color=ID),alpha=alpha_stream, linetype=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[4]][,3],color=ID, alpha=ID), 
                size=1.3,stat="smooth",method = "lm")+
      scale_alpha_manual(values = c(rep(0.3,6),0.9))+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "UndWei - STcon")+
      xlab("HOBO relative position Upstream-Downstream")+
      ylab("ST Connecticity")+
      scale_y_continuous(limits = c(0,60000))+
      theme_classic()+theme(legend.position="none"),
    top=""),
  
  #B1 NMDS plot
  arrangeGrob(
    ggplot(Un_WEIG_ST_MatrixOut)+
      geom_point(aes(x=x, y=y, fill=ID,size=DtoU), shape=21,colour="grey10", alpha=alpha_stream)+
      stat_ellipse(aes(x=x, y=y, colour=ID),type = "t", linetype=2,size=1, alpha=0.5)+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      scale_y_continuous(label=NULL)+
      scale_x_continuous(label=NULL)+
      labs(subtitle = "UndWei - STconmat")+
      theme_classic()+theme(legend.position="none",axis.ticks = element_blank())+
      facet_grid(. ~ ID),top = ""),
  

    get_legend(ggplot(NonW_ST_directed_out)+
                 geom_point(aes(x=DtoU,y=four_approxim[[1]][,3],fill=ID),shape=21, size=5)+
                 scale_fill_CUNILLERA(palette = "LGTBI", name="Stream ID")+
                 theme_classic()+theme(legend.direction = "horizontal",legend.box="vertical")),
    nrow=5 ,ncol=2,widths=c(0.6,1.1))
  dev.off()

Charging and depurating biological data

### IMPORTANT: TO PROPERLY RUN THIS SCRIPT YOU NEED TO RUN the following scripts: 

## SpaTemp_HOBOS_treatment.R - The most important one as is the one calculating the ST indices and 
# selecting and creating the HOBOS dataset of interest and then preparing all the data in order to be analysed. 

## Old_HOBOS_calculation.R - Not that important but interesting. Here the Old HOBOS indices are calculated. This 
# indices are calcualted based on the amounts of "1" or "0" of the HOBOS table. They are interesting as they offer
# an extra perspective on the whole ST dynamics. 
## IMPORTANT: If you do not want to use them you will need to remove the parts where they are "called". Generally,
# they are named "Old_HOBOS" or similar names (TotDur, TotNum and TotLeng).

## SpaTemp_comparison.R - It is not essential to run this script. However, it is interesting to check how the 
# HOBOS obtained indices behave among them. Specially to understand what they mean and what they do not mean. So, 
# my recomendation is to run it. 

# SUGGESTED ORDER TO RUN the scripts: 
### 1 SpaTemp_HOBOS_treatment.R
### 2 FlowIntermittence_indices.R
### 3 SpaTemp_comparison.R
### 4 SpaTemp_Biolcomparison.R

# Nice colors? CUNILLERA_palette is what you need
source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/CUNILLERA_palette.R")
setwd("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/1. FEHM coses al DIA/4. Mecodispers Spa-Tem/MECODISPER SpaTem")

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# 1. BIOLOGICAL Data uploading ####
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

# Autumn dataset - November 2019
BDD <- read.csv2("BiolData/Matriz_otoño.csv", sep=";")
# We correct an error as there are two SA4 for October sampling, we sum them. 
BDD[which(BDD$Code=="SA4")[1],3:ncol(BDD)] <- apply(BDD[which(BDD$Code=="SA4"),3:ncol(BDD)],2,sum)
BDD <- BDD[-which(BDD$Code=="SA4")[2],]
BDD <- mutate(BDD, 
              Riera = case_when(
                str_detect(Code,"CA") ~ "C", 
                str_detect(Code,"H") ~ "VH",
                str_detect(Code,"MU") ~ "M",
                str_detect(Code,"R") ~ "R",
                str_detect(Code,"SA") ~ "SA",
                str_detect(Code,"SC") ~ "SC",
                str_detect(Code,"T") ~ "T",
                TRUE ~ "ERROR" )) 


Traits_val <- read.csv2("BiolData/traits.csv", row.names = 1)
bloc<-read.table("BiolData/trait.blocks.txt",h=T)

# 1.1 TRAITS implementation based on rel. abundances of each trait ####
# Filter species that are present or not present existing 
out <- c()
for (i in 1:length(rownames(Traits_val))) {
  temp_out <- which(rownames(Traits_val)==colnames(BDD%>%dplyr::select(-Date,-Riera,-Code,))[i])
  out[i] <- ifelse(length(temp_out)>0,temp_out,0)
}
Trait_matr_prep<-prep.fuzzy.var(Traits_val[out,],bloc$bloc)
apply(Trait_matr_prep,1,sum)

BDD_croped <- BDD%>%dplyr::select(matches(rownames(Trait_matr_prep)))
rownames(BDD_croped) <-BDD$Code

for(i in 1:nrow(BDD_croped)){BDD_croped[i,]<- BDD_croped[i,]/sum(BDD_croped[i,])}
apply(BDD_croped,1,sum)

Trait_abundances<-as.matrix(BDD_croped)%*%as.matrix(Trait_matr_prep)
  
Trait_abund_f4 <- data.frame("Code"=BDD$Code,"Date"=BDD$Date,"Riera"=BDD$Riera, "f4"=as.data.frame(Trait_abundances)$f4)
Trait_abund_f1 <- data.frame("Code"=BDD$Code,"Date"=BDD$Date,"Riera"=BDD$Riera, "f1"=as.data.frame(Trait_abundances)$f1)

# 1.2 TRAITS by filtering ####
# Traits addition by  filtering only species corresponding to the traits
Traits_val_edited <- data.frame("taxon"=row.names(Traits_val),Traits_val)

# Active dispersers -- f4==3 and 2
Trait_disp <- Traits_val_edited%>%dplyr::select(taxon,f4)%>%filter(f4==3)%>%
  rows_insert(Traits_val_edited%>%dplyr::select(taxon,f4)%>%filter(f4==2))
BDD_activ <- BDD%>%dplyr::select(Code, Date, Riera,matches(Trait_disp$taxon))

# Passive dispersers -- f1==3 and 2
Trait_disp <- Traits_val_edited%>%dplyr::select(taxon,f1)%>%filter(f1==3)%>%
  rows_insert(Traits_val_edited%>%dplyr::select(taxon,f1)%>%filter(f1==2))
BDD_pass <- BDD%>%dplyr::select(Code, Date, Riera,matches(Trait_disp$taxon))

Active dispersers community metrics calculated

Data_season <- BDD_activ %>%
                dplyr::select(-Date,-Riera)%>%
                pivot_longer(-Code)

colnames(Data_season) <- c("Code","Species","Abund")

#INDIVIDUAL diversity values active ________________________________________________________________________####
by_group_data_season <- group_by(Data_season, 
                                  Code, Species)

data_G_gamma <- by_group_data_season%>%
                group_by(Species)%>%
                summarise(n=sum(Abund))%>%
                mutate(s=1)%>%
                summarise(rich=sum(s))

data_A_alpha <- by_group_data_season%>%
                group_by(Code)%>%
                mutate(s=ifelse(Abund>0,1,0))%>%
                summarise(rich=sum(s))

library(vegan)
table_shannon <- by_group_data_season%>%
                 spread(key =Species,value = Abund, fill=0)%>%
                 bind_cols("Riera"=BDD_activ$Riera)

shan <- vegan::diversity(table_shannon%>%ungroup()%>%dplyr::select(-Code,-Riera),index = "shannon")

#SIMMILARITIES active ________________________________________________________________________####

table_shannon_PA <- by_group_data_season%>%
  mutate(PA=ifelse(Abund>0,1,0))%>%
  dplyr::select(-Abund)%>%
  spread(key =Species,value = PA, fill=0)%>%
  bind_cols("Riera"=BDD_activ$Riera)

Bray.distance <- list()
Jac.distance <- list()
for (riera in 1:length(unique(table_shannon$Riera))) {
  names_riera <- unique(table_shannon$Riera[order(table_shannon$Riera)])
  
  # Bray Curtis active
  Bray.distance[[riera]] <- vegdist(
                          decostand(table_shannon%>%
                                      filter(Riera==names_riera[riera])%>%ungroup()%>%dplyr::select(-Code,-Riera),
                          method = "hellinger"),
                          method = "bray")
  
  # Jaccard active
  Jac.distance[[riera]]<-vegdist(
                        decostand(table_shannon_PA%>%
                                    filter(Riera==names_riera[riera])%>%ungroup()%>%dplyr::select(-Code,-Riera),
                        method = "pa"),
                        method = "jaccard")
  }

names(Bray.distance) <- unique(table_shannon$Riera[order(table_shannon$Riera)])
names(Jac.distance) <- unique(table_shannon$Riera[order(table_shannon$Riera)])

BID_output_IND_act <- data.frame(Code=unique(by_group_data_season$Code),rich_f4=data_A_alpha$rich,shan_f4=shan,f4_abun=Trait_abund_f4$f4)
BID_output_DIST_act <- list(Bray.distance,Jac.distance)

Passive dispersers community metrics calculated

Data_season <- BDD_pass %>%
  dplyr::select(-Date,-Riera)%>%
  pivot_longer(-Code)

colnames(Data_season) <- c("Code","Species","Abund")

#INDIVIDUAL diversity values passive ________________________________________________________________________####
by_group_data_season <- group_by(Data_season, 
                                 Code, Species)

data_G_gamma <- by_group_data_season%>%
  group_by(Species)%>%
  summarise(n=sum(Abund))%>%
  mutate(s=1)%>%
  summarise(rich=sum(s))

data_A_alpha <- by_group_data_season%>%
  group_by(Code)%>%
  mutate(s=ifelse(Abund>0,1,0))%>%
  summarise(rich=sum(s))

library(vegan)
table_shannon <- by_group_data_season%>%
  spread(key =Species,value = Abund, fill=0)%>%
  bind_cols("Riera"=BDD_pass$Riera)

shan <- vegan::diversity(table_shannon%>%ungroup()%>%dplyr::select(-Code,-Riera),index = "shannon")

#SIMMILARITIES passive ________________________________________________________________________####

table_shannon_PA <- by_group_data_season%>%
  mutate(PA=ifelse(Abund>0,1,0))%>%
  dplyr::select(-Abund)%>%
  spread(key =Species,value = PA, fill=0)%>%
  bind_cols("Riera"=BDD_pass$Riera)

Bray.distance <- list()
Jac.distance <- list()
for (riera in 1:length(unique(table_shannon$Riera))) {
  names_riera <- unique(table_shannon$Riera[order(table_shannon$Riera)])
  
  # Bray Curtis
  Bray.distance[[riera]] <- vegdist(
    decostand(table_shannon%>%
                filter(Riera==names_riera[riera])%>%ungroup()%>%dplyr::select(-Code,-Riera),
              method = "hellinger"),
    method = "bray")
  
  # Jaccard
  Jac.distance[[riera]]<-vegdist(
    decostand(table_shannon_PA%>%
                filter(Riera==names_riera[riera])%>%ungroup()%>%dplyr::select(-Code,-Riera),
              method = "pa"),
    method = "jaccard")
}

names(Bray.distance) <- unique(table_shannon$Riera[order(table_shannon$Riera)])
names(Jac.distance) <- unique(table_shannon$Riera[order(table_shannon$Riera)])

BID_output_IND_pas <- data.frame(Code=unique(by_group_data_season$Code),rich_f1=data_A_alpha$rich,shan_f1=shan,f1_abun=Trait_abund_f1$f1)
BID_output_DIST_pas <- list(Bray.distance,Jac.distance)

Merging all community metrics calculated

# BID values total
BID_output_IND <- left_join(BID_output_IND_act, BID_output_IND_pas, by="Code")

BID_output_DIST <- list("Active"=BID_output_DIST_act,"Passive"=BID_output_DIST_pas)

Charging geographic and matching biological samples and data logger localization

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# CHARGING GEOGRAPHIC COORDINATES & MATCHING SAMPLES AND HOBOS ####
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# 3.1 We charge the field Samples and filter them according to the ones actually taken (from BDD)  ####
SampSites <- read.table("BiolData/Sites_1.txt",sep = ",",header = T)
out <- c()
for (match in 1:length(BDD$Code)) {
 a <- which(SampSites$Site==BDD$Code[match])
 out[match] <- ifelse(length(a)==0,0,a)
}
SampSites <- SampSites[out,]

# 3.2 We call the ST indices to MATCH them  ####
## You need to run the previous script in order to built the HOBOS dataset (a list with the HOBOS info for each river)
Sites_list # They are charged from SpaTemp_HOBOS_treatment.R
HOBOS_sites# They are charged from SpaTemp_HOBOS_treatment.R

# This two other objects are also built in the "SpaTemp_HOBOS_treatment.R" script. They contain the names 
## of each river and the order upstream to downstream within it. 
# They are used to Identify and locate river position between HOBOS and real Samples
HOB_riv_ID<- Stream_order$ï..Riera
ups_dos <- Stream_order$UtoD

# Calculation of OLD HOBOS values to include them in the comparisons. 
source("FlowIntermittence_indices.R")
Old_HOBOS_comb

# 3.3 Geographical matching ####
# Converting Sampling SITES to apropriate geog. data 
str(SampSites)
a <- c()
b <- c()
for (posit in 1:nrow(SampSites)) {
a[posit] <-as.numeric(sub("(.{2})(.*)", "\\1.\\2", SampSites$lat[posit]))
b[posit] <-as.numeric(sub("(.{1})(.*)", "\\1.\\2", SampSites$long[posit]))
}
SampSites$lat<- a
SampSites$long<- b

# Converting HOBOS SITES to appropriate geog. data 
### *IMPORTANT STEP ####
#Here is where we include the calculated individual ST indices by "left_join"
Sites_list_comb <- rbind(Sites_list[[1]],Sites_list[[2]],Sites_list[[3]],Sites_list[[4]],Sites_list[[5]],
                         Sites_list[[6]],Sites_list[[7]] )
Sites_list_comb <- cbind(Sites_list_comb, "ID"=HOB_riv_ID,"DtoU"=ups_dos)%>%
                   mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))
Sites_list_comb <- Sites_list_comb%>%
                   left_join(NonW_ST_directed_out)%>%
                   left_join(WEIG_ST_directed_out)%>%
                   left_join(Un_NonW_ST_out)%>%
                   left_join(Un_WEIG_ST_out)%>%
                   left_join(Old_HOBOS_comb)
a <- c()
b <- c()
for (posit in 1:nrow(Sites_list_comb)) {
  a[posit] <-as.numeric(sub("(.{2})(.*)", "\\1.\\2", as.character(round(Sites_list_comb$Latitud[posit]))))
  b[posit] <-as.numeric(sub("(.{1})(.*)", "\\1.\\2", as.character(round(Sites_list_comb$Longitud[posit]))))
}
Sites_list_comb$Latitud<- a
Sites_list_comb$Longitud<- b

# Plot to check that the values match. 
par(mfrow=c(1,2))
plot(Sites_list_comb$Longitud, Sites_list_comb$Latitud)
plot(SampSites$long,SampSites$lat)
par(mfrow=c(1,1))

# Match the values
library(RANN) 
Matching_HOBOS <- nn2(data=Sites_list_comb[,3:4], query = SampSites[,2:3], k=1, radius = 1)[[1]] 

### *MANUAL CHECKING ####
# Here you must make sure that the assignation is correct between Samp_ID, ID, Riera, and DtoU
HOB_BDD_match <- as.data.frame(cbind(Sites_list_comb[Matching_HOBOS,], "Samp_ID"=SampSites$Site)) # Check assignation
data.frame(
HOB_BDD_match$Riera,
HOB_BDD_match$Codi_HOBO,
HOB_BDD_match$DtoU,
HOB_BDD_match$Samp_ID)

# Order the sites by upstream downstream directionality
HOB_BDD_match <- HOB_BDD_match%>%arrange(Riera,DtoU)


#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# MERGING ALL THE VALUES TOGETHER ####
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

#Individual
HOB_BDD_match <- HOB_BDD_match%>%left_join(BID_output_IND, by=c("Samp_ID"="Code"))

#Matrices
### *IMPORTANT STEP ####
#Here is where we include the calculated matrix ST indices
NonW_Dir <- list()
WEIG_Dir<- list()
Un_NonW<- list()
Un_WEIG<- list()
for (riera in 1:length(unique(NonW_ST_matrix_out_out))) {
  names_riera <- unique(Sites_list_comb$ID)
  
  values_places <- unlist(HOB_BDD_match%>%filter(ID==names_riera[riera])%>%dplyr::select(DtoU))
  print(values_places)
  
  NonW_Dir[[riera]] <- as.dist(t(NonW_ST_matrix_out_out[[riera]][values_places,values_places]))
  WEIG_Dir[[riera]] <- as.dist(t(WEIG_ST_matrix_out_out[[riera]][values_places,values_places]))
  Un_NonW[[riera]] <- as.dist(t(Un_NonW_ST_matrix_out_out[[riera]][values_places,values_places]))
  Un_WEIG[[riera]] <- as.dist(t(Un_WEIG_ST_matrix_out_out[[riera]][values_places,values_places]))
}

### *MANUAL CHECKING ####
# We eliminate R and G matrices due to the lack of sites for this two sampling places
## For diversities (Bray and Jaccard) eliminate positions 2 and 4
## For ST matrices eliminate positions 3
HOB_BDD_match_matrix <- list("BrayCurtis_Act"=BID_output_DIST$Active[[1]][-c(2,4)],
                             "Jaccard_Act"=BID_output_DIST$Active[[2]][-c(2,4)],
                             "BrayCurtis_Pas"=BID_output_DIST$Passive[[1]][-c(2,4)],
                             "Jaccard_Pas"=BID_output_DIST$Passive[[2]][-c(2,4)],
                             "NonW_Dir"=NonW_Dir[-c(3)], 
                             "WEIG_Dir"=WEIG_Dir[-c(3)], 
                             "Un_NonW"=Un_NonW[-c(3)], 
                             "Un_WEIG"=Un_WEIG[-c(3)])

# We eliminate R and G in Individual due to the lack of samples from this period
HOB_BDD_match <- HOB_BDD_match%>%filter(Samp_ID !=c("R4","G1"))

# FINAL MERGED VALUES - INDIVIDUAL AND DISTANCE (MATRIX) BASED 
HOB_BDD_match
HOB_BDD_match_matrix

### *SMALL EDITIONS TO ELIMINATE Ocl/Acl/Bet ####
HOB_BDD_match <- HOB_BDD_match%>%dplyr::select(-c(NonW_Dir_Ocl,NonW_Dir_Acl,NonW_Dir_Bet,
                                           WEIG_Dir_Ocl,WEIG_Dir_Acl,WEIG_Dir_Bet,
                                           Un_NonW_Ocl,Un_NonW_Acl,Un_NonW_Bet,
                                           Un_WEIG_Ocl,Un_WEIG_Acl,Un_WEIG_Bet))

Richness

plots_HOB_BDD_total <- list()
plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
for (col_var in 1:ncol(HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                              -Latitud,-Longitud,
                                              -ID, -DtoU,-ID_UpDo,
                                              -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                              -rich_f4, -shan_f4,-f4_abun,
                                              -rich_f1,-shan_f1,-f1_abun))){

  
  variable_x_temp <- HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                       -Latitud,-Longitud,
                                       -ID, -DtoU,-ID_UpDo,
                                       -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                       -rich_f4, -shan_f4,-f4_abun,
                                       -rich_f1,-shan_f1,-f1_abun)

  variable_x_name <- colnames(variable_x_temp)[col_var]
  variable_x <- variable_x_temp[,col_var]
  
  variable_y_temp <- HOB_BDD_match%>%dplyr::select(rich_f4,rich_f1)
  plot_counter <- c(col_var, col_var+ncol(variable_x_temp))
  
for (variable in 1:ncol(variable_y_temp)) {
    
  variable_y <- variable_y_temp[,variable]
  variable_y_model <- unlist(variable_y)
  variable_y_name <- colnames(variable_y_temp)
  
  factors <- HOB_BDD_match%>%dplyr::select(ID, DtoU,ID_UpDo)
  Id_random <- unlist(factors$ID)
  
  model <- lme(log(variable_y_model+1)~log(variable_x+1), random = ~1|Id_random)
  results <- round(as.numeric(c(summary(model)[[4]]$fixed, summary(model)[[20]][2,5])),2)
  #Results table
  model_table <- round(summary(model)[[20]],3)
  rownames(model_table) <- c("Intercept", "STconmat")
  model_table<- rownames_to_column(as.data.frame(model_table))
  model_table<- cbind(Scenario=names_scenarios[col_var],Dispersion=Dispersal_group[variable],model_table)
  model_HOB_BDD_results[[plot_counter[variable]]] <- model_table 
  
  dataset <- cbind(factors, "X_var"=variable_x, variable_y)
  dataset$pred_values <- predict(model)   
  
plots_HOB_BDD[[plot_counter[variable]]] <- ggplot(dataset,aes(x=log(X_var+1), y=log(variable_y+1),fill=ID,colour=ID))+
    geom_point(color="grey30", alpha=0.5, shape=21, size=2)+
    geom_smooth(method = "lm", colour="grey60",alpha=0.1, se = TRUE, linetype=2)+
    geom_abline(slope =results[2],intercept = results[1], colour="black", size=2)+
    #geom_smooth(method = "lm",alpha=0.2, se = T, fill="grey50", colour="black", size=2)+
    scale_fill_manual(values =  CunilleraPAL_corrected)+
    scale_color_manual(values =  CunilleraPAL_corrected)+
    xlab(paste(axisXtitles[col_var]))+
    ylab(paste("log(Richness)"))+
    labs(caption = paste("Intercept=",results[1],"Slope=",results[2],"p-value=",results[3]))+
    labs(subtitle=paste(names_scenarios[col_var],"vs",variable_y_name[variable]))+
    theme_classic()+
    theme(legend.position="none",
          axis.line.x.bottom=element_line(color=axistextcolours[col_var], size=3),
          axis.text.x = element_text(color =axistextcolours[col_var] ),
          axis.title.x = element_text(color =axistextcolours[col_var] , size=8 ))
    #facet_grid(.~ID)  
  }
}

model_result <- rbind(model_HOB_BDD_results[[1]],model_HOB_BDD_results[[2]],
                      model_HOB_BDD_results[[3]],model_HOB_BDD_results[[4]],
                      model_HOB_BDD_results[[5]],model_HOB_BDD_results[[6]],
                      model_HOB_BDD_results[[7]],model_HOB_BDD_results[[8]])
colnames(model_result)[[3]] <- "Fixed effects"
write.table(model_result, "Table_Results/Biol_Rich_STcon.txt", sep = ",")

legend_plots<- get_legend(ggplot(dataset)+
                            geom_point(aes(x=X_var,y=variable_y ,fill=ID),shape=21, size=6)+
                            scale_fill_manual(values = CunilleraPAL_corrected, name="Stream ID")+
                            theme_classic()+theme(legend.direction = "vertical",legend.box="vertical"))


png(filename ="Figure/Biol_treat/Richness.png", 
    width = 1250*4, height = 1000*2, 
    units = "px",res = 300) 
grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  legend_plots,
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=5, top="Richness")
dev.off()
plots_HOB_BDD_total[[1]] <- plots_HOB_BDD

Shannon-Wiener

plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
for (col_var in 1:ncol(HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                              -Latitud,-Longitud,
                                              -ID, -DtoU,-ID_UpDo,
                                              -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                              -rich_f4, -shan_f4,-f4_abun,
                                              -rich_f1,-shan_f1,-f1_abun))){
  
  
  variable_x_temp <- HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                            -Latitud,-Longitud,
                                            -ID, -DtoU,-ID_UpDo,
                                            -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                            -rich_f4, -shan_f4,-f4_abun,
                                            -rich_f1,-shan_f1,-f1_abun)
  
  variable_x_name <- colnames(variable_x_temp)[col_var]
  variable_x <- variable_x_temp[,col_var]
  
  variable_y_temp <- HOB_BDD_match%>%dplyr::select(shan_f4,shan_f1)
  plot_counter <- c(col_var, col_var+ncol(variable_x_temp))
  
  for (variable in 1:ncol(variable_y_temp)) {
    
    variable_y <- variable_y_temp[,variable]
    variable_y_model <- unlist(variable_y)
    variable_y_name <- colnames(variable_y_temp)
    
    factors <- HOB_BDD_match%>%dplyr::select(ID, DtoU,ID_UpDo)
    Id_random <- unlist(factors$ID)
    
    model <- lme(log(variable_y_model+1)~log(variable_x+1), random = ~1|Id_random)
    results <- round(as.numeric(c(summary(model)[[4]]$fixed, summary(model)[[20]][2,5])),2)
    #Results table
    model_table <- round(summary(model)[[20]],3)
    rownames(model_table) <- c("Intercept", "STconmat")
    model_table<- rownames_to_column(as.data.frame(model_table))
    model_table<- cbind(Scenario=names_scenarios[col_var],Dispersion=Dispersal_group[variable],model_table)
    model_HOB_BDD_results[[plot_counter[variable]]] <- model_table 
    
    dataset <- cbind(factors, "X_var"=variable_x, variable_y)
    dataset$pred_values <- predict(model)   
    
    plots_HOB_BDD[[plot_counter[variable]]] <- ggplot(dataset,aes(x=log(X_var+1), y=log(variable_y+1),fill=ID,colour=ID))+
      geom_point(color="grey30", alpha=0.5, shape=21, size=2)+
      geom_smooth(method = "lm", colour="grey60",alpha=0.1, se = TRUE, linetype=2)+
      geom_abline(slope =results[2],intercept = results[1], colour="black", size=2)+
      #geom_smooth(method = "lm",alpha=0.2, se = T, fill="grey50", colour="black", size=2)+
      scale_fill_manual(values =  CunilleraPAL_corrected)+
      scale_color_manual(values =  CunilleraPAL_corrected)+
      xlab(paste(axisXtitles[col_var]))+
      ylab(paste("log(Shannon)"))+
      labs(caption = paste("Intercept=",results[1],"Slope=",results[2],"p-value=",results[3]))+
      labs(subtitle=paste(names_scenarios[col_var],"vs",variable_y_name[variable]))+
      theme_classic()+
      theme(legend.position="none",
            axis.line.x.bottom=element_line(color=axistextcolours[col_var], size=3),
            axis.text.x = element_text(color =axistextcolours[col_var] ),
            axis.title.x = element_text(color =axistextcolours[col_var], size=8  ))
    #facet_grid(.~ID)  
  }
}

model_result <- rbind(model_HOB_BDD_results[[1]],model_HOB_BDD_results[[2]],
                      model_HOB_BDD_results[[3]],model_HOB_BDD_results[[4]],
                      model_HOB_BDD_results[[5]],model_HOB_BDD_results[[6]],
                      model_HOB_BDD_results[[7]],model_HOB_BDD_results[[8]])
colnames(model_result)[[3]] <- "Fixed effects"
write.table(model_result, "Table_Results/Biol_ShaWei_STcon.txt", sep = ",")

legend_plots<- get_legend(ggplot(dataset)+
                            geom_point(aes(x=X_var,y=variable_y ,fill=ID),shape=21, size=6)+
                            scale_fill_manual(values = CunilleraPAL_corrected, name="Stream ID")+
                            theme_classic()+theme(legend.direction = "vertical",legend.box="vertical"))


png(filename ="Figure/Biol_treat/Shannon.png", 
    width = 1250*4, height = 1000*2, 
    units = "px",res = 300) 
grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  legend_plots,
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=5, top="Shannon")
dev.off()

plots_HOB_BDD_total[[2]] <- plots_HOB_BDD

Trait abundance

plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
for (col_var in 1:ncol(HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                              -Latitud,-Longitud,
                                              -ID, -DtoU,-ID_UpDo,
                                              -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                              -rich_f4, -shan_f4,-f4_abun,
                                              -rich_f1,-shan_f1,-f1_abun))){
  
  
  variable_x_temp <- HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                            -Latitud,-Longitud,
                                            -ID, -DtoU,-ID_UpDo,
                                            -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                            -rich_f4, -shan_f4,-f4_abun,
                                            -rich_f1,-shan_f1,-f1_abun)
  
  variable_x_name <- colnames(variable_x_temp)[col_var]
  variable_x <- variable_x_temp[,col_var]
  
  variable_y_temp <- HOB_BDD_match%>%dplyr::select(f4_abun,f1_abun)
  plot_counter <- c(col_var, col_var+ncol(variable_x_temp))
  
  for (variable in 1:ncol(variable_y_temp)) {
    
    variable_y <- variable_y_temp[,variable]
    variable_y_model <- unlist(variable_y)
    variable_y_name <- colnames(variable_y_temp)
    
    factors <- HOB_BDD_match%>%dplyr::select(ID, DtoU,ID_UpDo)
    Id_random <- unlist(factors$ID)
    
    model <- lme(log(variable_y_model+1)~log(variable_x+1), random = ~1|Id_random)
    results <- round(as.numeric(c(summary(model)[[4]]$fixed, summary(model)[[20]][2,5])),2)
    #Results table
    model_table <- round(summary(model)[[20]],3)
    rownames(model_table) <- c("Intercept", "STconmat")
    model_table<- rownames_to_column(as.data.frame(model_table))
    model_table<- cbind(Scenario=names_scenarios[col_var],Dispersion=Dispersal_group[variable],model_table)
    model_HOB_BDD_results[[plot_counter[variable]]] <- model_table 
    
    dataset <- cbind(factors, "X_var"=variable_x, variable_y)
    dataset$pred_values <- predict(model)   
    
    plots_HOB_BDD[[plot_counter[variable]]] <- ggplot(dataset,aes(x=log(X_var+1), y=log(variable_y+1),fill=ID,colour=ID))+
      geom_point(color="grey30", alpha=0.5, shape=21, size=2)+
      geom_smooth(method = "lm", colour="grey60",alpha=0.1, se = TRUE, linetype=2)+
      geom_abline(slope =results[2],intercept = results[1], colour="black", size=2)+
      #geom_smooth(method = "lm",alpha=0.2, se = T, fill="grey50", colour="black", size=2)+
      scale_fill_manual(values =  CunilleraPAL_corrected)+
      scale_color_manual(values =  CunilleraPAL_corrected)+
      xlab(paste(axisXtitles[col_var]))+
      ylab(paste("log(Trait abundance)"))+
      labs(caption = paste("Intercept=",results[1],"Slope=",results[2],"p-value=",results[3]))+
      labs(subtitle=paste(names_scenarios[col_var],"vs",variable_y_name[variable]))+
      theme_classic()+
      theme(legend.position="none",
            axis.line.x.bottom=element_line(color=axistextcolours[col_var], size=3),
            axis.text.x = element_text(color =axistextcolours[col_var] ),
            axis.title.x = element_text(color =axistextcolours[col_var], size=8  ))
    #facet_grid(.~ID)  
  }
}

model_result <- rbind(model_HOB_BDD_results[[1]],model_HOB_BDD_results[[2]],
                      model_HOB_BDD_results[[3]],model_HOB_BDD_results[[4]],
                      model_HOB_BDD_results[[5]],model_HOB_BDD_results[[6]],
                      model_HOB_BDD_results[[7]],model_HOB_BDD_results[[8]])
colnames(model_result)[[3]] <- "Fixed effects"
write.table(model_result, "Table_Results/Biol_TraitAb_STcon.txt", sep = ",")

legend_plots<- get_legend(ggplot(dataset)+
                            geom_point(aes(x=X_var,y=variable_y ,fill=ID),shape=21, size=6)+
                            scale_fill_manual(values = CunilleraPAL_corrected, name="Stream ID")+
                            theme_classic()+theme(legend.direction = "vertical",legend.box="vertical"))


png(filename ="Figure/Biol_treat/TraitAbun.png", 
    width = 1250*4, height = 1000*2, 
    units = "px",res = 300) 
grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  legend_plots,
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=5, top="Trait abundance")
dev.off()

plots_HOB_BDD_total[[3]] <- plots_HOB_BDD

Matching Bray & Jaccard with STconmat values

output <- matrix(nrow=length(unlist(HOB_BDD_match_matrix$BrayCurtis_Act)))
for (lis_elem in 1:length(HOB_BDD_match_matrix)) {
  out <- c()
  for (river in 1:length(HOB_BDD_match_matrix$BrayCurtis_Act)) {
    col_names_river <- names(HOB_BDD_match_matrix$BrayCurtis_Act)
    
    out_temp <- c()
    out_temp <- cbind(as.matrix(HOB_BDD_match_matrix[[lis_elem]][[river]])[
      upper.tri(as.matrix(HOB_BDD_match_matrix[[lis_elem]][[river]]))])
    out_temp <- as.data.frame(cbind("ID"=rep(col_names_river[river],nrow(out_temp)),out_temp))
    out<- rbind(out, out_temp)
  }
  colnames(out)[2] <- names(HOB_BDD_match_matrix)[lis_elem]
  output[,1] <- out[,1]
  output <- as.data.frame(cbind(output, out))
  output <- output%>%dplyr::select(-ID)
  output[,lis_elem+1] <- as.numeric(output[,lis_elem+1])
}
colnames(output)[1] <- "ID"

STmatrix_BiolDissim <- output

axisXtitles <- c(c("Isolated   <-- STconmat -->  Connected"), 
                 c("Low disp. resist <-- STconmat --> High disp. resist"),
                 c("Isolated   <-- STconmat -->  Connected"), 
                 c("Low disp. resist <-- STconmat --> High disp. resist"))

leng_dist_matrix <- ncol(as.matrix(HOB_BDD_match_matrix$BrayCurtis_Act[[1]]))+
ncol(as.matrix(HOB_BDD_match_matrix$BrayCurtis_Act[[2]]))+
ncol(as.matrix(HOB_BDD_match_matrix$BrayCurtis_Act[[3]]))+
ncol(as.matrix(HOB_BDD_match_matrix$BrayCurtis_Act[[4]]))+
ncol(as.matrix(HOB_BDD_match_matrix$BrayCurtis_Act[[5]]))+
ncol(as.matrix(HOB_BDD_match_matrix$BrayCurtis_Act[[6]]))

matrices_total_analyse <- list()
for (matrices_tot in 1:length(HOB_BDD_match_matrix)) {
Matrix_Beta <- HOB_BDD_match_matrix[[matrices_tot]]
full_matrix <- matrix(nrow = leng_dist_matrix,ncol = leng_dist_matrix,data = 0)
for (matrix_number in 1:length(Matrix_Beta)) {
ncol_leng <- ncol(as.matrix(Matrix_Beta[[matrix_number]]))
mac_ncol_leng <- which(apply(full_matrix,2,sum)==0)[ncol_leng]
min_ncol_leng <- which(apply(full_matrix,2,sum)==0)[1]
full_matrix[min_ncol_leng:mac_ncol_leng,
            min_ncol_leng:mac_ncol_leng] <- as.matrix(Matrix_Beta[[matrix_number]])
}
zeros <- which(full_matrix==0)
full_matrix[zeros] <- mean(full_matrix[-zeros])
matrices_total_analyse[[matrices_tot]] <- full_matrix
}
names(matrices_total_analyse) <- names(HOB_BDD_match_matrix)

Bray-Curtis simmilarity

plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
variable_y_name <- c("BrayCurtis_Act","BrayCurtis_Pas")
Dispersal_group <- c("Active", "Passive")
# BrayCurtis_Act
out_results <- data.frame()
for (stconmat_val in 1:4) {
a <- vegan::mantel(ydis = as.dist(matrices_total_analyse$BrayCurtis_Act),
                   xdis = as.dist(matrices_total_analyse[[4+stconmat_val]]),permutations = 1000,method = "spearman")
results <- data.frame("Group"=Dispersal_group[1],
                      "Scenario"=names_scenarios[stconmat_val],
                      "Mantel coefficient"=round(a$statistic,3), "p.value"=round(a$signif,3))

plots_HOB_BDD[[stconmat_val]] <- data.frame(
"X"=matrices_total_analyse[[4+stconmat_val]][lower.tri(matrices_total_analyse[[4+stconmat_val]])],
"Y"=matrices_total_analyse$BrayCurtis_Act[lower.tri(matrices_total_analyse$BrayCurtis_Act)]) %>% 
  ggplot(aes(x=X,y=Y))+
  geom_point(color="grey20",fill="grey40",alpha=0.8, shape=21, size=2)+
  geom_smooth(method = "lm", colour="black",alpha=0.1, se = TRUE, linetype=1, linewidth=2)+
  xlab(paste(axisXtitles[stconmat_val]))+
  ylab(paste("Bray curtis"))+
  labs(caption = paste("Mantel coefficient=",results[3],"p.value=",results[4]))+
  labs(subtitle=paste(names_scenarios[stconmat_val],"vs",variable_y_name[1]))+
  theme_classic()+
  theme(legend.position="none",
        axis.line.x.bottom=element_line(color=axistextcolours[stconmat_val], linewidth=3),
        axis.text.x = element_text(color =axistextcolours[stconmat_val] ),
        axis.title.x = element_text(color =axistextcolours[stconmat_val] , size=8 ))
out_results <- bind_rows(out_results, results)
}

# BrayCurtis_Pas
for (stconmat_val in 1:4) {
a <- vegan::mantel(ydis = as.dist(matrices_total_analyse$BrayCurtis_Pas),
                   xdis = as.dist((matrices_total_analyse[[4+stconmat_val]])),permutations = 1000,method = "spearman")
results <- data.frame("Group"=Dispersal_group[2],
                      "Scenario"=names_scenarios[stconmat_val],
                      "Mantel coefficient"=round(a$statistic,3), "p.value"=round(a$signif,3))

plots_HOB_BDD[[stconmat_val+4]] <- data.frame(
  "X"=(matrices_total_analyse[[4+stconmat_val]][lower.tri(matrices_total_analyse[[4+stconmat_val]])]),
  "Y"=matrices_total_analyse$BrayCurtis_Pas[lower.tri(matrices_total_analyse$BrayCurtis_Pas)]) %>% 
  ggplot(aes(x=X,y=Y))+
  geom_point(color="grey20",fill="grey40",alpha=0.8, shape=21, size=2)+
  geom_smooth(method = "lm", colour="black",alpha=0.1, se = TRUE, linetype=1, linewidth=2)+
  xlab(paste(axisXtitles[stconmat_val]))+
  ylab(paste("Bray curtis"))+
  labs(caption = paste("Mantel coefficient=",results[3],"p.value=",results[4]))+
  labs(subtitle=paste(names_scenarios[stconmat_val],"vs",variable_y_name[2]))+
  theme_classic()+
  theme(legend.position="none",
        axis.line.x.bottom=element_line(color=axistextcolours[stconmat_val], linewidth=3),
        axis.text.x = element_text(color =axistextcolours[stconmat_val] ),
        axis.title.x = element_text(color =axistextcolours[stconmat_val] , size=8 ))
out_results <- bind_rows(out_results, results)
}

model_result <- out_results
write.table(model_result, "Table_Results/Biol_BrayC_STconmat.txt", sep = ",")

png(filename ="Figure/Biol_treat/BrayCurtis.png", 
    width = 1250*4, height = 1000*2, 
    units = "px",res = 300) 
grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=4, top="Bray Crutis")
dev.off()
plots_HOB_BDD_total[[4]] <- plots_HOB_BDD

Jaccard

plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
variable_y_name <- c("Jaccard_Act","Jaccard_Pas")
Dispersal_group <- c("Active", "Passive")
# Jaccard_Act
out_results <- data.frame()
for (stconmat_val in 1:4) {
  a <- vegan::mantel(ydis = as.dist(matrices_total_analyse$Jaccard_Act),
                     xdis = as.dist(matrices_total_analyse[[4+stconmat_val]]),permutations = 1000,method = "spearman")
  results <- data.frame("Group"=Dispersal_group[1],
                        "Scenario"=names_scenarios[stconmat_val],
                        "Mantel coefficient"=round(a$statistic,3), "p.value"=round(a$signif,3))
  
  plots_HOB_BDD[[stconmat_val]] <- data.frame(
    "X"=matrices_total_analyse[[4+stconmat_val]][lower.tri(matrices_total_analyse[[4+stconmat_val]])],
    "Y"=matrices_total_analyse$Jaccard_Act[lower.tri(matrices_total_analyse$Jaccard_Act)]) %>% 
    ggplot(aes(x=X,y=Y))+
    geom_point(color="grey20",fill="grey40",alpha=0.8, shape=21, size=2)+
    geom_smooth(method = "lm", colour="black",alpha=0.1, se = TRUE, linetype=1, linewidth=2)+
    xlab(paste(axisXtitles[stconmat_val]))+
    ylab(paste("Jaccard"))+
    labs(caption = paste("Mantel coefficient=",results[3],"p.value=",results[4]))+
    labs(subtitle=paste(names_scenarios[stconmat_val],"vs",variable_y_name[1]))+
    theme_classic()+
    theme(legend.position="none",
          axis.line.x.bottom=element_line(color=axistextcolours[stconmat_val], linewidth=3),
          axis.text.x = element_text(color =axistextcolours[stconmat_val] ),
          axis.title.x = element_text(color =axistextcolours[stconmat_val] , size=8 ))
  out_results <- bind_rows(out_results, results)
}

# Jaccard_Pas
for (stconmat_val in 1:4) {
  a <- vegan::mantel(ydis = as.dist(matrices_total_analyse$Jaccard_Pas),
                     xdis = as.dist((matrices_total_analyse[[4+stconmat_val]])),permutations = 1000,method = "spearman")
  results <- data.frame("Group"=Dispersal_group[2],
                        "Scenario"=names_scenarios[stconmat_val],
                        "Mantel coefficient"=round(a$statistic,3), "p.value"=round(a$signif,3))
  
  plots_HOB_BDD[[stconmat_val+4]] <- data.frame(
    "X"=(matrices_total_analyse[[4+stconmat_val]][lower.tri(matrices_total_analyse[[4+stconmat_val]])]),
    "Y"=matrices_total_analyse$Jaccard_Pas[lower.tri(matrices_total_analyse$Jaccard_Pas)]) %>% 
    ggplot(aes(x=X,y=Y))+
    geom_point(color="grey20",fill="grey40",alpha=0.8, shape=21, size=2)+
    geom_smooth(method = "lm", colour="black",alpha=0.1, se = TRUE, linetype=1, linewidth=2)+
    xlab(paste(axisXtitles[stconmat_val]))+
    ylab(paste("Jaccard"))+
    labs(caption = paste("Mantel coefficient=",results[3],"p.value=",results[4]))+
    labs(subtitle=paste(names_scenarios[stconmat_val],"vs",variable_y_name[2]))+
    theme_classic()+
    theme(legend.position="none",
          axis.line.x.bottom=element_line(color=axistextcolours[stconmat_val], linewidth=3),
          axis.text.x = element_text(color =axistextcolours[stconmat_val] ),
          axis.title.x = element_text(color =axistextcolours[stconmat_val] , size=8 ))
  out_results <- bind_rows(out_results, results)
}

model_result <- out_results
write.table(model_result, "Table_Results/Biol_Jacc_STconmat.txt", sep = ",")

png(filename ="Figure/Biol_treat/Jaccard.png", 
    width = 1250*4, height = 1000*2, 
    units = "px",res = 300) 
grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=4, top="Jaccard")
dev.off()
plots_HOB_BDD_total[[5]] <- plots_HOB_BDD