diff --git a/outputs/CommunitySizeSpectra/CSSOLSfishcommunityBenthicPelagic.rda b/outputs/CommunitySizeSpectra/CSSOLSfishcommunity.rda
similarity index 100%
rename from outputs/CommunitySizeSpectra/CSSOLSfishcommunityBenthicPelagic.rda
rename to outputs/CommunitySizeSpectra/CSSOLSfishcommunity.rda
diff --git a/outputs/CommunitySizeSpectra/CSSOLSfrequencytableBenthicPelagic.rda b/outputs/CommunitySizeSpectra/CSSOLSfrequencytable.rda
similarity index 100%
rename from outputs/CommunitySizeSpectra/CSSOLSfrequencytableBenthicPelagic.rda
rename to outputs/CommunitySizeSpectra/CSSOLSfrequencytable.rda
diff --git a/outputs/SEM/sem_details_indirect_effect_all_lakes.rda b/outputs/SEM/sem_details_indirect_effect_all_lakes.rda
new file mode 100644
index 0000000..b4672ac
Binary files /dev/null and b/outputs/SEM/sem_details_indirect_effect_all_lakes.rda differ
diff --git a/outputs/SEM/sem_effect_cold_lakes.rda b/outputs/SEM/sem_effect_cold_lakes.rda
deleted file mode 100644
index 385ff0d..0000000
Binary files a/outputs/SEM/sem_effect_cold_lakes.rda and /dev/null differ
diff --git a/outputs/SEM/sem_effect_moderate_lakes.rda b/outputs/SEM/sem_effect_moderate_lakes.rda
deleted file mode 100644
index b976a4b..0000000
Binary files a/outputs/SEM/sem_effect_moderate_lakes.rda and /dev/null differ
diff --git a/outputs/SEM/sem_effect_warm_lakes.rda b/outputs/SEM/sem_effect_warm_lakes.rda
deleted file mode 100644
index affe1e4..0000000
Binary files a/outputs/SEM/sem_effect_warm_lakes.rda and /dev/null differ
diff --git a/outputs/Temperature/datasetTempCurrentBioclim_258L.rda b/outputs/Temperature/datasetTempCurrentBioclim.rda
similarity index 100%
rename from outputs/Temperature/datasetTempCurrentBioclim_258L.rda
rename to outputs/Temperature/datasetTempCurrentBioclim.rda
diff --git a/outputs/Temperature/datasetTempTrendBioclim40_258L.rda b/outputs/Temperature/datasetTempTrendBioclim40.rda
similarity index 100%
rename from outputs/Temperature/datasetTempTrendBioclim40_258L.rda
rename to outputs/Temperature/datasetTempTrendBioclim40.rda
diff --git a/rscripts/09_Computation_Mann-Kendall_temperature_trend_bioclim.R b/rscripts/09_Computation_Mann-Kendall_temperature_trend_bioclim.R
index c96b2c4..746be77 100644
--- a/rscripts/09_Computation_Mann-Kendall_temperature_trend_bioclim.R
+++ b/rscripts/09_Computation_Mann-Kendall_temperature_trend_bioclim.R
@@ -39,7 +39,7 @@ daily_temp_okp = as.data.frame(daily_temp_okp)
 
 
 ##get read of useless symbols
-daily_temp_okp$Lake_ID <- str_sub(daily_temp_okp$Lake_ID, start = 48, end = 52)
+daily_temp_okp$Lake_ID <- str_sub(daily_temp_okp$Lake_ID, start = 25, end = 29)
 daily_temp_okp$Lake_ID[daily_temp_okp$Lake_ID == "GRO21"] <- "GRO21b"
 daily_temp_okp$Lake_ID[daily_temp_okp$Lake_ID == "LER27"] <- "LER27a"
 head(daily_temp_okp)
@@ -54,7 +54,7 @@ daily_temp_okp <- daily_temp_okp %>% dplyr::select(Lake_ID, date, year, month, d
 rm(list_of_files)
 
 
-date_fish_sampling <- read.csv("./B_Thermal_trajectories/outputs/date_fish_sampling.txt", sep="")
+date_fish_sampling <- read.csv("./outputs/date_fish_sampling.txt", sep="")
 length(unique(date_fish_sampling$code_lac))  
 
 
@@ -133,7 +133,6 @@ for(i in 1:nrow(date_fish_sampling)){
     sub_monthly_temp_okp <- monthly_temp_okp %>% dplyr::filter(Lake_ID == date_fish_sampling[i, 1])
     sub_monthly_temp_okp <- sub_monthly_temp_okp %>% dplyr::filter(date >= format(as.Date(ymd(date_fish_sampling$date_pose[i]) - years(period[j])), "%Y-%m"))
     sub_monthly_temp_okp <- sub_monthly_temp_okp %>% dplyr::filter(date < format(as.Date(date_fish_sampling$date_pose[i]), "%Y-%m"))
-
     sub_monthly_temp_okp$year2 <- rep(as.numeric(min(sub_monthly_temp_okp$year)):(max(as.numeric(max(sub_monthly_temp_okp$year)))-1), each = 12)
     sub_monthly_temp_okp <- sub_monthly_temp_okp %>% dplyr::select(Lake_ID, date, year2, tmean.epi, tmin.epi, tmax.epi)
     colnames(sub_monthly_temp_okp) <- c("Lake_ID", "date", "year", "tmean.epi", "tmin.epi", "tmax.epi")
@@ -340,13 +339,13 @@ colnames(temporal_trend_bio11) <- c("lake.code", "year.samp", "id.samp", "bio11.
 
 
 
-datasetTempTrendBioclim40_258L <- left_join(temporal_trend_bio1, temporal_trend_bio2 %>% dplyr::select(id.samp, bio2.slope.50y, bio2.pvalue.50y), by = "id.samp")
-datasetTempTrendBioclim40_258L <- left_join(datasetTempTrendBioclim40_258L, temporal_trend_bio3 %>% dplyr::select(id.samp, bio3.slope.50y, bio3.pvalue.50y), by = "id.samp")
-datasetTempTrendBioclim40_258L <- left_join(datasetTempTrendBioclim40_258L, temporal_trend_bio4 %>% dplyr::select(id.samp, bio4.slope.50y, bio4.pvalue.50y), by = "id.samp")
-datasetTempTrendBioclim40_258L <- left_join(datasetTempTrendBioclim40_258L, temporal_trend_bio5 %>% dplyr::select(id.samp, bio5.slope.50y, bio5.pvalue.50y), by = "id.samp")
-datasetTempTrendBioclim40_258L <- left_join(datasetTempTrendBioclim40_258L, temporal_trend_bio6 %>% dplyr::select(id.samp, bio6.slope.50y, bio6.pvalue.50y), by = "id.samp")
-datasetTempTrendBioclim40_258L <- left_join(datasetTempTrendBioclim40_258L, temporal_trend_bio7 %>% dplyr::select(id.samp, bio7.slope.50y, bio7.pvalue.50y), by = "id.samp")
-datasetTempTrendBioclim40_258L <- left_join(datasetTempTrendBioclim40_258L, temporal_trend_bio10 %>% dplyr::select(id.samp, bio10.slope.50y, bio10.pvalue.50y), by = "id.samp")
-datasetTempTrendBioclim40_258L <- left_join(datasetTempTrendBioclim40_258L, temporal_trend_bio11 %>% dplyr::select(id.samp, bio11.slope.50y, bio11.pvalue.50y), by = "id.samp")
+datasetTempTrendBioclim40 <- left_join(temporal_trend_bio1, temporal_trend_bio2 %>% dplyr::select(id.samp, bio2.slope.50y, bio2.pvalue.50y), by = "id.samp")
+datasetTempTrendBioclim40 <- left_join(datasetTempTrendBioclim40, temporal_trend_bio3 %>% dplyr::select(id.samp, bio3.slope.50y, bio3.pvalue.50y), by = "id.samp")
+datasetTempTrendBioclim40 <- left_join(datasetTempTrendBioclim40, temporal_trend_bio4 %>% dplyr::select(id.samp, bio4.slope.50y, bio4.pvalue.50y), by = "id.samp")
+datasetTempTrendBioclim40 <- left_join(datasetTempTrendBioclim40, temporal_trend_bio5 %>% dplyr::select(id.samp, bio5.slope.50y, bio5.pvalue.50y), by = "id.samp")
+datasetTempTrendBioclim40 <- left_join(datasetTempTrendBioclim40, temporal_trend_bio6 %>% dplyr::select(id.samp, bio6.slope.50y, bio6.pvalue.50y), by = "id.samp")
+datasetTempTrendBioclim40 <- left_join(datasetTempTrendBioclim40, temporal_trend_bio7 %>% dplyr::select(id.samp, bio7.slope.50y, bio7.pvalue.50y), by = "id.samp")
+datasetTempTrendBioclim40 <- left_join(datasetTempTrendBioclim40, temporal_trend_bio10 %>% dplyr::select(id.samp, bio10.slope.50y, bio10.pvalue.50y), by = "id.samp")
+datasetTempTrendBioclim40 <- left_join(datasetTempTrendBioclim40, temporal_trend_bio11 %>% dplyr::select(id.samp, bio11.slope.50y, bio11.pvalue.50y), by = "id.samp")
 
-mysave(datasetTempTrendBioclim40_258L, dir = "./outputs/Temperature", overwrite = TRUE)
+mysave(datasetTempTrendBioclim40, dir = "./outputs/Temperature", overwrite = TRUE)
diff --git a/rscripts/11_Computation_current_mean_temperature.R b/rscripts/09_Computation_current_mean_temperature.R
similarity index 88%
rename from rscripts/11_Computation_current_mean_temperature.R
rename to rscripts/09_Computation_current_mean_temperature.R
index 65385f4..472a440 100644
--- a/rscripts/11_Computation_current_mean_temperature.R
+++ b/rscripts/09_Computation_current_mean_temperature.R
@@ -37,7 +37,7 @@ daily_temp_okp = as.data.frame(daily_temp_okp)
 
 
 ##get read of useless symbols
-daily_temp_okp$Lake_ID <- str_sub(daily_temp_okp$Lake_ID, start = 48, end = 52)
+daily_temp_okp$Lake_ID <- str_sub(daily_temp_okp$Lake_ID, start = 25, end = 29)
 daily_temp_okp$Lake_ID[daily_temp_okp$Lake_ID == "GRO21"] <- "GRO21b"
 daily_temp_okp$Lake_ID[daily_temp_okp$Lake_ID == "LER27"] <- "LER27a"
 head(daily_temp_okp)
@@ -218,13 +218,13 @@ rm(i, j, period)
 
 
 
-datasetTempCurrentBioclim_258L <- left_join(bio1, bio2 %>% dplyr::select(id.samp, bio2.current), by = "id.samp")
-datasetTempCurrentBioclim_258L <- left_join(datasetTempCurrentBioclim_258L, bio3 %>% dplyr::select(id.samp, bio3.current), by = "id.samp")
-datasetTempCurrentBioclim_258L <- left_join(datasetTempCurrentBioclim_258L, bio4 %>% dplyr::select(id.samp, bio4.current), by = "id.samp")
-datasetTempCurrentBioclim_258L <- left_join(datasetTempCurrentBioclim_258L, bio5 %>% dplyr::select(id.samp, bio5.current), by = "id.samp")
-datasetTempCurrentBioclim_258L <- left_join(datasetTempCurrentBioclim_258L, bio6 %>% dplyr::select(id.samp, bio6.current), by = "id.samp")
-datasetTempCurrentBioclim_258L <- left_join(datasetTempCurrentBioclim_258L, bio7 %>% dplyr::select(id.samp, bio7.current), by = "id.samp")
-datasetTempCurrentBioclim_258L <- left_join(datasetTempCurrentBioclim_258L, bio10 %>% dplyr::select(id.samp, bio10.current), by = "id.samp")
-datasetTempCurrentBioclim_258L <- left_join(datasetTempCurrentBioclim_258L, bio11 %>% dplyr::select(id.samp, bio11.current), by = "id.samp")
+datasetTempCurrentBioclim <- left_join(bio1, bio2 %>% dplyr::select(id.samp, bio2.current), by = "id.samp")
+datasetTempCurrentBioclim <- left_join(datasetTempCurrentBioclim, bio3 %>% dplyr::select(id.samp, bio3.current), by = "id.samp")
+datasetTempCurrentBioclim <- left_join(datasetTempCurrentBioclim, bio4 %>% dplyr::select(id.samp, bio4.current), by = "id.samp")
+datasetTempCurrentBioclim <- left_join(datasetTempCurrentBioclim, bio5 %>% dplyr::select(id.samp, bio5.current), by = "id.samp")
+datasetTempCurrentBioclim <- left_join(datasetTempCurrentBioclim, bio6 %>% dplyr::select(id.samp, bio6.current), by = "id.samp")
+datasetTempCurrentBioclim <- left_join(datasetTempCurrentBioclim, bio7 %>% dplyr::select(id.samp, bio7.current), by = "id.samp")
+datasetTempCurrentBioclim <- left_join(datasetTempCurrentBioclim, bio10 %>% dplyr::select(id.samp, bio10.current), by = "id.samp")
+datasetTempCurrentBioclim <- left_join(datasetTempCurrentBioclim, bio11 %>% dplyr::select(id.samp, bio11.current), by = "id.samp")
 
-mysave(datasetTempCurrentBioclim_258L, dir = "./outputs/Temperature", overwrite = TRUE)
+mysave(datasetTempCurrentBioclim, dir = "./outputs/Temperature", overwrite = TRUE)
diff --git a/rscripts/10_Compute_native_vs_nonindigenous_diversity.R b/rscripts/10_Compute_native_vs_nonindigenous_diversity.R
index 3469f27..72990d7 100644
--- a/rscripts/10_Compute_native_vs_nonindigenous_diversity.R
+++ b/rscripts/10_Compute_native_vs_nonindigenous_diversity.R
@@ -27,7 +27,7 @@ colnames(df.speciestype) ; colnames(df.speciestype)[1] <- "species"
 
 #--
 
-myload(ind_size, dir = "data/IndividualSize")
+myload(ind_size, dir = "outputs/IndividualSize")
 
 head(ind_size)
 df.fish <- ind_size %>% dplyr::select(id_campagne, species, fish)
diff --git a/rscripts/11_Computation_community_size_spectra_binsmethod.R b/rscripts/11_Computation_community_size_spectra_binsmethod.R
index 7670de6..7e0650a 100644
--- a/rscripts/11_Computation_community_size_spectra_binsmethod.R
+++ b/rscripts/11_Computation_community_size_spectra_binsmethod.R
@@ -88,6 +88,6 @@ for(i in 1:length(id)){
   
   rm(sub.df, size, SS.parameters, frequencytable)
 }
-CSSOLSfishcommunityBenthicPelagic <- CSSOLSfishcommunity
-CSSOLSfrequencytableBenthicPelagic <- CSSOLSfrequencytable
-mysave(CSSOLSfishcommunityBenthicPelagic, CSSOLSfrequencytableBenthicPelagic, dir = "./outputs/CommunitySizeSpectra", overwrite = TRUE)
+CSSOLSfishcommunity <- CSSOLSfishcommunity
+CSSOLSfrequencytable <- CSSOLSfrequencytable
+mysave(CSSOLSfishcommunity, CSSOLSfrequencytable, dir = "./outputs/CommunitySizeSpectra", overwrite = TRUE)
diff --git a/rscripts/12_SEManalysis.R b/rscripts/12_SEManalysis.R
index 1059c81..7491f77 100644
--- a/rscripts/12_SEManalysis.R
+++ b/rscripts/12_SEManalysis.R
@@ -285,6 +285,6 @@ data.effect$significant <- NA
 data.effect$significant[data.effect$SymbSign == "*"] <- TRUE
 data.effect$significant[data.effect$SymbSign == ""] <- FALSE
 data.effect <- data.effect %>% mutate(position = if_else(Effect > 0, UpperCI + 0.02, LowerCI - 0.02))
-#sem_effect_all_lakes <- data.effect
-#mysave(sem_effect_all_lakes, dir = "outputs/SEM", overwrite = TRUE)
+sem_effect_all_lakes <- data.effect
+mysave(sem_effect_all_lakes, dir = "outputs/SEM", overwrite = TRUE)
 
diff --git a/rscripts/13_Figure1.R b/rscripts/13_Figure1.R
index e3e304e..5fc1529 100644
--- a/rscripts/13_Figure1.R
+++ b/rscripts/13_Figure1.R
@@ -72,44 +72,13 @@ coord <- dataset.thermal.trajectories %>%
                    bio1.slope.40y = mean(bio1.slope.40y),
                    bio1.current = mean(bio1.current))
 str(coord)
-coord.sf <- st_as_sf(coord, coords = c("long", "lat"), crs = 4326)
-data <- bi_class(coord.sf, x = bio1.slope.40y, y = bio1.current, style = "quantile", dim = 4)
-
-
-map <- ggplot() +
-  geom_sf(data = world_joined, fill = "white", color = "black", size = 0.05) +
-  geom_sf(data = francemap, fill = gray(0.9), color = "black", size = 0.25) +
-  geom_sf(data = francerivers, col = '#6baed6', size = 0.25) +  
-  geom_sf(data = francelakes, col = '#6baed6', fill = '#6baed6', size = 0.05) + 
-  geom_sf(data = data, mapping = aes(fill = bi_class),  shape = 21, size = 2, color = "black", show.legend = FALSE) +
-  bi_scale_fill(pal = "Brown2", dim = 4) +
-  annotation_scale(location = "bl", width_hint = 0.1) +
-  annotation_north_arrow(which_north = "true", location = "tr", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(fill = c("black", "black"), text_size = 6)) +           
-  coord_sf(xlim = c(-5, 9.75), ylim = c(41.3, 51.5), expand = FALSE) +
-  map_theme +
-  theme(strip.background = element_rect(color = "black", size = 1, linetype = "solid"),
-        strip.text.x = element_text(size = 12, color = "black", face = "bold"),
-        panel.border = element_rect(colour = "black", fill = NA, size = 1)) + theme(legend.position = "none")
-map
-#4.5 x 4.5
-
-legend <- bi_legend(pal = "Brown2",
-                    dim = 4,
-                    xlab = "bio1.slope.40y",
-                    ylab = "bio1.current",
-                    #rotate_pal = TRUE,
-                    #flip_axes = TRUE,
-                    size = 8)
-ggdraw() + draw_plot(legend, 0.2, .65, 0.2, 0.2)
-#"Bluegill", "BlueGold", "BlueOr", "BlueYl", "Brown"/"Brown2", "DkBlue"/"DkBlue2", "DkCyan"/"DkCyan2", "DkViolet"/"DkViolet2", "GrPink"/"GrPink2", "PinkGrn", "PurpleGrn", or "PurpleOr".
-
-
-data$bio1.slope.40y <- data$bio1.slope.40y * 10 #in°C/dec
-data$class.bio1.slope.40y <- NA
-data$class.bio1.slope.40y[data$bio1.slope.40y < 0] <- "A"
-data$class.bio1.slope.40y[data$bio1.slope.40y >= 0 & data$bio1.slope.40y < 0.25] <- "B"
-data$class.bio1.slope.40y[data$bio1.slope.40y >= 0.25 & data$bio1.slope.40y < 0.50] <- "C"
-data$class.bio1.slope.40y[data$bio1.slope.40y >= 0.50] <- "D"
+coord <- st_as_sf(coord, coords = c("long", "lat"), crs = 4326)
+coord$bio1.slope.40y <- coord$bio1.slope.40y * 10 #in°C/dec
+coord$class.bio1.slope.40y <- NA
+coord$class.bio1.slope.40y[coord$bio1.slope.40y < 0] <- "A"
+coord$class.bio1.slope.40y[coord$bio1.slope.40y >= 0 & coord$bio1.slope.40y < 0.25] <- "B"
+coord$class.bio1.slope.40y[coord$bio1.slope.40y >= 0.25 & coord$bio1.slope.40y < 0.50] <- "C"
+coord$class.bio1.slope.40y[coord$bio1.slope.40y >= 0.50] <- "D"
 
 
 
@@ -118,7 +87,7 @@ map <- ggplot() +
   geom_sf(data = francemap, fill = gray(0.9), color = "black", size = 0.25) +
   geom_sf(data = francerivers, col = '#6baed6', size = 0.25) +  
   geom_sf(data = francelakes, col = '#6baed6', fill = '#6baed6', size = 0.05) + 
-  geom_sf(data = data, mapping = aes(fill = bio1.current, shape = class.bio1.slope.40y), size = 2, color = "black", show.legend = FALSE) +
+  geom_sf(data = coord, mapping = aes(fill = bio1.current, shape = class.bio1.slope.40y), size = 2, color = "black", show.legend = FALSE) +
   scale_shape_manual(values = c(25, 21, 22, 23)) +
   scale_fill_continuous_sequential(palette = "Heat") +
   annotation_scale(location = "bl", width_hint = 0.1) +
@@ -381,7 +350,7 @@ summary(mod.size)
 
 
 
-plot_grid(pTemp, psizespecies, 
+plot_grid(map, psizespecies, 
           pslope, pmidpoint,
           pconnect, pmaxtl,
           labels = c("A", "B", "C", "D", "E", "F"),
diff --git a/rscripts/13_Figure3.R b/rscripts/13_Figure3.R
index 30184db..a8ff5bf 100644
--- a/rscripts/13_Figure3.R
+++ b/rscripts/13_Figure3.R
@@ -27,11 +27,10 @@ myload(sem_effect_all_lakes,
        dir = "outputs/SEM")
 
 
-##--------
-## FIGURES
-##--------
+##---------------------------
+## FIGURE ABOUT TOTAL EFFECTS
+##---------------------------
 
-##FIGURE ABOUT TOTAL EFFECTS
 data.total.effect <- sem_effect_all_lakes %>% dplyr::filter(Type == "Total")
 
 
@@ -145,279 +144,3 @@ plot_grid(pCAMT, pTAMT,
           labels = c("A", "B", "C", "D", "E", "F"),
           ncol = 2, nrow = 3, align = "hv")
 #landscape 10x15
-
-
-
-
-#FIGURE ABOUT DIRECT AND INDIRECT EFFECTS
-#https://stackoverflow.com/questions/62393159/how-can-i-add-hatches-stripes-or-another-pattern-or-texture-to-a-barplot-in-ggp
-data.direct.indirect.effect <- sem_effect_all_lakes %>% dplyr::filter(!Type == "Total")
-
-
-pCAMT <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Climatic conditions"),
-                aes(x = Effect, y = Response, pattern = factor(Type))) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  #geom_bar(aes(alpha = factor(Type)), stat = "identity", position = "dodge", width = 0.75, colour = "black", fill = "#8d96a3") +
-  #geom_bar_pattern(stat='identity', position = position_dodge(preserve = "single"),
-  #                 color = "black", 
-  #                 pattern_fill = "black",
-  #                 fill = "#8d96a3",
-  #                 pattern_angle = 45,
-  #                 pattern_density = 0.1,
-  #                 pattern_spacing = 0.025,
-  #                 pattern_key_scale_factor = 0.6) + 
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.35, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effects", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23)) +
-  scale_fill_manual(values = c("#b8bec6", "#393f47")) +
-  #scale_pattern_manual(values = c(Indirect = "stripe", Direct = "none")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000")) 
-pCAMT
-
-
-pTAMT <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Climate warming"),
-                aes(x = Effect, y = Response, pattern = factor(Type))) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.35, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effects", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23)) +
-  scale_fill_manual(values = c("#b8bec6", "#393f47")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pTAMT
-
-
-pNISR <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Exotic sp. richness"),
-                aes(x = Effect, y = Response, pattern = factor(Type))) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.25, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effects", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23)) +
-  scale_fill_manual(values = c("#b0cebe", "#314f3f")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pNISR
-
-
-pTSR <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Total sp. richness"),
-               aes(x = Effect, y = Response, pattern = factor(Type))) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.25, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effects", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23)) +
-  scale_fill_manual(values = c("#b0cebe", "#314f3f")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pTSR
-
-
-pSlope <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Slope"),
-                 aes(x = Effect, y = Response, pattern = factor(Type))) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.15, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effects", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23)) +
-  scale_fill_manual(values = c("#f4cb8b", "#744c0b")) + 
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pSlope
-
-
-pMP <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Elevation"),
-              aes(x = Effect, y = Response, pattern = factor(Type))) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.15, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effects", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23)) +
-  scale_fill_manual(values = c("#f4cb8b", "#744c0b")) + 
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pMP
-
-
-
-plot_grid(pCAMT, pTAMT, 
-          pNISR, pTSR,
-          pSlope, pMP,
-          labels = c("A", "B", "C", "D", "E", "F"),
-          ncol = 2, nrow = 3, align = "hv")
-#landscape 10x15
-
-
-
-
-#FIGURE ABOUT ALL EFFECTS (TOTAL, DIRECT AND INDIRECT)
-pCAMT <- ggplot(sem_effect_all_lakes %>% dplyr::filter(Predictor == "Climatic conditions"),
-                aes(x = Effect, y = Response)) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  #geom_bar(aes(alpha = factor(Type)), stat = "identity", position = "dodge", width = 0.75, colour = "black", fill = "#8d96a3") +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.35, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effect size", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23, 21)) +
-  scale_fill_manual(values = c("#b8bec6", "#393f47", "#8d96a3")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000")) 
-pCAMT
-
-
-pTAMT <- ggplot(sem_effect_all_lakes %>% dplyr::filter(Predictor == "Climate warming"),
-                aes(x = Effect, y = Response)) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.35, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effect size", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23, 21)) +
-  scale_fill_manual(values = c("#b8bec6", "#393f47", "#8d96a3")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pTAMT
-
-
-pNISR <- ggplot(sem_effect_all_lakes %>% dplyr::filter(Predictor == "Exotic sp. richness"),
-                aes(x = Effect, y = Response)) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.25, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effect size", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23, 21)) +
-  scale_fill_manual(values = c("#b0cebe", "#314f3f", "#66a182")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pNISR
-
-
-pTSR <- ggplot(sem_effect_all_lakes %>% dplyr::filter(Predictor == "Total sp. richness"),
-               aes(x = Effect, y = Response)) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.25, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effect size", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23, 21)) +
-  scale_fill_manual(values = c("#b0cebe", "#314f3f", "#66a182")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pTSR
-
-
-pSlope <- ggplot(sem_effect_all_lakes %>% dplyr::filter(Predictor == "Slope"),
-                 aes(x = Effect, y = Response)) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.15, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effect size", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23, 21)) +
-  scale_fill_manual(values = c("#f4cb8b", "#744c0b", "#edae49")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pSlope
-
-
-pMP <- ggplot(sem_effect_all_lakes %>% dplyr::filter(Predictor == "Elevation"),
-              aes(x = Effect, y = Response)) +
-  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
-  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.15, fill = "black") +
-  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
-  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
-  labs(x = "Standardized effect size", y = "Reponse variable") +
-  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
-  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
-  scale_shape_manual(values = c(22, 23, 21)) +
-  scale_fill_manual(values = c("#f4cb8b", "#744c0b", "#edae49")) +
-  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
-  theme_classic() + 
-  theme(legend.position = "none",
-        axis.text = element_text(size = 14, colour = "#000000"),
-        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-pMP
-
-
-
-ggarrange(pCAMT, pTAMT, 
-          pNISR, pTSR,
-          pSlope, pMP,
-          labels = c("A", "B", "C", "D", "E", "F"),
-          ncol = 2, nrow = 3,
-          common.legend = TRUE, legend="bottom")
-
-plot_grid(pCAMT, pTAMT, 
-          pNISR, pTSR,
-          pSlope, pMP,
-          labels = c("A", "B", "C", "D", "E", "F"),
-          ncol = 2, nrow = 3, align = "hv")
-
-#landscape 10x15
-
diff --git a/rscripts/13_FigureSI1.R b/rscripts/13_FigureSI1.R
index 65fc96b..8f3fe41 100644
--- a/rscripts/13_FigureSI1.R
+++ b/rscripts/13_FigureSI1.R
@@ -58,9 +58,9 @@ ind_size <- as.data.frame(ind_size)
 
 ##MEAN SIZE OF SPECIES
 species_mean_size <- ind_size %>%
-  select(species, fish) %>%
-  group_by(species) %>%
-  summarise(mean_size = mean(fish)/10)
+  dplyr::select(species, fish) %>%
+  dplyr::group_by(species) %>%
+  dplyr::summarise(mean_size = mean(fish)/10)
 quantile(species_mean_size$mean_size, na.rm = T, probs = c(0.33333, 0.6666))
 species_mean_size$size_category <- NA
 species_mean_size$size_category[species_mean_size$mean_size <= 12.5] <- "< 12.5 cm"
@@ -75,47 +75,47 @@ species_mean_size$size_category[is.na(species_mean_size$size_category)] <- "12.5
 species_occurrence <- ind_size %>% select(code_lac, species) %>% unique(.)
 species_occurrence$presence <- 1
 species_occurrence <- species_occurrence %>%
-  select(species, presence) %>%
-  group_by(species) %>%
-  summarise(occurrence = sum(presence))
+  dplyr::select(species, presence) %>%
+  dplyr::group_by(species) %>%
+  dplyr::summarise(occurrence = sum(presence))
 species_occurrence$percentage <- (species_occurrence$occurrence / nrow(ind_size %>% select(code_lac) %>% unique(.)))*100
 
 
 ##NUMBER OF TIMES SPECIES REACH THE MAXIMUM TROPHIC LEVEL 
-length(unique(network_lake_metrics$id_campagne)) ; length(unique(dataset.thermal.trajectories$id.samp))
-network_lake_metrics <- network_lake_metrics %>% dplyr::filter(id_campagne %in% dataset.thermal.trajectories$id.samp)
+#length(unique(network_lake_metrics$id_campagne)) ; length(unique(dataset.thermal.trajectories$id.samp))
+#network_lake_metrics <- network_lake_metrics %>% dplyr::filter(id_campagne %in% dataset.thermal.trajectories$id.samp)
 
-species_maxTL <- data.frame(matrix(nrow = 0, ncol = 5))
-colnames(species_maxTL) <- c("id_campagne", "trophic_species", "code_species", "species", "TL")
+#species_maxTL <- data.frame(matrix(nrow = 0, ncol = 5))
+#colnames(species_maxTL) <- c("id_campagne", "trophic_species", "code_species", "species", "TL")
 
-for(i in 361:nrow(network_lake_metrics)){
+#for(i in 1:nrow(network_lake_metrics)){
 #issue i = 18 ; 269 ; 360
 
-sub <- network_lake_metrics[[6]][[i]]
-sub <- sub %>% top_n(1, obs_troph_level)
+#sub <- network_lake_metrics[[6]][[i]]
+#sub <- sub %>% top_n(1, obs_troph_level)
 
-sub_species_maxTL <- data.frame(matrix(nrow = 1, ncol = 5))
-colnames(sub_species_maxTL) <- c("id_campagne", "trophic_species", "code_species", "species", "TL")
-sub_species_maxTL[1, 1] <- network_lake_metrics[[1]][[i]]
-sub_species_maxTL[1, 2] <- sub$species_name[2]
-sub_species_maxTL[1, 3] <- substr(sub$species_name[2], 0, 3)
-sub_species_maxTL[1, 4] <- dfspeciestype %>% dplyr::filter(code == substr(sub$species_name[2], 0, 3)) %>% dplyr::select(species)
-sub_species_maxTL[1, 5] <- sub$obs_troph_level[2]
+#sub_species_maxTL <- data.frame(matrix(nrow = 1, ncol = 5))
+#colnames(sub_species_maxTL) <- c("id_campagne", "trophic_species", "code_species", "species", "TL")
+#sub_species_maxTL[1, 1] <- network_lake_metrics[[1]][[i]]
+#sub_species_maxTL[1, 2] <- sub$species_name[2]
+#sub_species_maxTL[1, 3] <- substr(sub$species_name[2], 0, 3)
+#sub_species_maxTL[1, 4] <- dfspeciestype %>% dplyr::filter(code == substr(sub$species_name[2], 0, 3)) %>% dplyr::select(species)
+#sub_species_maxTL[1, 5] <- sub$obs_troph_level[2]
 
-species_maxTL <- rbind(species_maxTL, sub_species_maxTL)
+#species_maxTL <- rbind(species_maxTL, sub_species_maxTL)
 
-rm(sub, sub_species_maxTL)
+#rm(sub, sub_species_maxTL)
 
-}
-rm(i)
+#}
+#rm(i)
 
-species_maxTL$count <- 1
+#species_maxTL$count <- 1
 
-maxTL <- species_maxTL %>%
-  select(species, count) %>%
-  group_by(species) %>%
-  summarise(maxTL = sum(count))
-maxTL$maxTL.percentage <- (maxTL$maxTL/length(unique(dataset.thermal.trajectories$id.samp)))*100
+#maxTL <- species_maxTL %>%
+#  select(species, count) %>%
+#  group_by(species) %>%
+#  summarise(maxTL = sum(count))
+#maxTL$maxTL.percentage <- (maxTL$maxTL/length(unique(dataset.thermal.trajectories$id.samp)))*100
 
 
 ##COMBINE ALL INFORMATIONS
diff --git a/rscripts/13_FigureSI2.R b/rscripts/13_FigureSI2.R
index 9e3e03a..ca77eb3 100644
--- a/rscripts/13_FigureSI2.R
+++ b/rscripts/13_FigureSI2.R
@@ -56,9 +56,9 @@ native <- ind_size %>% dplyr::filter(type == "native") %>% dplyr::select(fish)
 mean(nis[,1])/mean(native[,1])
 
 summary_ind_size <- ind_size %>%
-  select(id_campagne, type, fish) %>%
-  group_by(id_campagne, type) %>%
-  summarise(mean_size = mean(fish),
+  dplyr::select(id_campagne, type, fish) %>%
+  dplyr::group_by(id_campagne, type) %>%
+  dplyr::summarise(mean_size = mean(fish),
             median_size = median(fish),
             max_size = max(fish))
 summary_ind_size_nis <- summary_ind_size %>% dplyr::filter(type == "non-indigenous") %>% as.data.frame(.)
diff --git a/rscripts/13_FigureSI3.R b/rscripts/13_FigureSI3.R
new file mode 100644
index 0000000..2d5fa0f
--- /dev/null
+++ b/rscripts/13_FigureSI3.R
@@ -0,0 +1,168 @@
+rm(list = ls()) #Removes all objects from the current workspace (R memory)
+
+
+##------------------------------
+##LOADING PACKAGES AND FUNCTIONS
+##------------------------------
+##PACKAGES##
+library(cowplot)
+library(dplyr)
+library(ggh4x)
+library(ggpattern)
+library(ggplot2)
+library(ggpubr)
+library(lemon)
+library(stringr)
+library(tidyr)
+
+
+##FUNCTIONS##
+source("./rfunctions/misc.R")
+
+
+##--------------
+## LOAD DATASETS
+##--------------
+myload(sem_effect_all_lakes,
+       dir = "outputs/SEM")
+
+
+##-----------------------------------------
+## FIGURE ABOUT DIRECT AND INDIRECT EFFECTS
+##-----------------------------------------
+#https://stackoverflow.com/questions/62393159/how-can-i-add-hatches-stripes-or-another-pattern-or-texture-to-a-barplot-in-ggp
+
+data.direct.indirect.effect <- sem_effect_all_lakes %>% dplyr::filter(!Type == "Total")
+
+
+pCAMT <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Climatic conditions"),
+                aes(x = Effect, y = Response, pattern = factor(Type))) +
+  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
+  #geom_bar(aes(alpha = factor(Type)), stat = "identity", position = "dodge", width = 0.75, colour = "black", fill = "#8d96a3") +
+  #geom_bar_pattern(stat='identity', position = position_dodge(preserve = "single"),
+  #                 color = "black", 
+  #                 pattern_fill = "black",
+  #                 fill = "#8d96a3",
+  #                 pattern_angle = 45,
+  #                 pattern_density = 0.1,
+  #                 pattern_spacing = 0.025,
+  #                 pattern_key_scale_factor = 0.6) + 
+  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.35, fill = "black") +
+  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
+  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
+  labs(x = "Standardized effects", y = "Reponse variable") +
+  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
+  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
+  scale_shape_manual(values = c(22, 23)) +
+  scale_fill_manual(values = c("#b8bec6", "#393f47")) +
+  #scale_pattern_manual(values = c(Indirect = "stripe", Direct = "none")) +
+  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
+  theme_classic() + 
+  theme(legend.position = "none",
+        axis.text = element_text(size = 14, colour = "#000000"),
+        axis.title = element_text(size = 16, face = "bold", colour = "#000000")) 
+pCAMT
+
+
+pTAMT <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Climate warming"),
+                aes(x = Effect, y = Response, pattern = factor(Type))) +
+  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
+  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.35, fill = "black") +
+  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
+  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
+  labs(x = "Standardized effects", y = "Reponse variable") +
+  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
+  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
+  scale_shape_manual(values = c(22, 23)) +
+  scale_fill_manual(values = c("#b8bec6", "#393f47")) +
+  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
+  theme_classic() + 
+  theme(legend.position = "none",
+        axis.text = element_text(size = 14, colour = "#000000"),
+        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
+pTAMT
+
+
+pNISR <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Exotic sp. richness"),
+                aes(x = Effect, y = Response, pattern = factor(Type))) +
+  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
+  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.25, fill = "black") +
+  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
+  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
+  labs(x = "Standardized effects", y = "Reponse variable") +
+  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
+  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
+  scale_shape_manual(values = c(22, 23)) +
+  scale_fill_manual(values = c("#b0cebe", "#314f3f")) +
+  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
+  theme_classic() + 
+  theme(legend.position = "none",
+        axis.text = element_text(size = 14, colour = "#000000"),
+        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
+pNISR
+
+
+pTSR <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Total sp. richness"),
+               aes(x = Effect, y = Response, pattern = factor(Type))) +
+  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
+  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.25, fill = "black") +
+  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
+  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
+  labs(x = "Standardized effects", y = "Reponse variable") +
+  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
+  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
+  scale_shape_manual(values = c(22, 23)) +
+  scale_fill_manual(values = c("#b0cebe", "#314f3f")) +
+  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
+  theme_classic() + 
+  theme(legend.position = "none",
+        axis.text = element_text(size = 14, colour = "#000000"),
+        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
+pTSR
+
+
+pSlope <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Slope"),
+                 aes(x = Effect, y = Response, pattern = factor(Type))) +
+  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
+  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.15, fill = "black") +
+  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
+  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
+  labs(x = "Standardized effects", y = "Reponse variable") +
+  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
+  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
+  scale_shape_manual(values = c(22, 23)) +
+  scale_fill_manual(values = c("#f4cb8b", "#744c0b")) + 
+  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
+  theme_classic() + 
+  theme(legend.position = "none",
+        axis.text = element_text(size = 14, colour = "#000000"),
+        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
+pSlope
+
+
+pMP <- ggplot(data.direct.indirect.effect %>% dplyr::filter(Predictor == "Elevation"),
+              aes(x = Effect, y = Response, pattern = factor(Type))) +
+  geom_vline(xintercept = 0, linetype = 'dashed', color = 'black', size = 1) +
+  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI, fill = Type), position = position_dodge(0.75), width = 0.15, fill = "black") +
+  geom_point(aes(shape = Type, fill = Type), position = position_dodge(0.75), width = 0.25, colour = "black", size = 4) +
+  geom_text(aes(label = ifelse(significant, "*", ""), x = position, fill = Type), position = position_dodge(0.75), size = 20 / .pt) +
+  labs(x = "Standardized effects", y = "Reponse variable") +
+  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6), limits = c(-0.4, 0.627), position = "top") +
+  coord_flex_cart(left=brackets_vertical(), top = capped_horizontal('both')) +
+  scale_shape_manual(values = c(22, 23)) +
+  scale_fill_manual(values = c("#f4cb8b", "#744c0b")) + 
+  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line()) +
+  theme_classic() + 
+  theme(legend.position = "none",
+        axis.text = element_text(size = 14, colour = "#000000"),
+        axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
+pMP
+
+
+
+plot_grid(pCAMT, pTAMT, 
+          pNISR, pTSR,
+          pSlope, pMP,
+          labels = c("A", "B", "C", "D", "E", "F"),
+          ncol = 2, nrow = 3, align = "hv")
+#landscape 10x15
diff --git a/rscripts/13_FigureSI5.R b/rscripts/13_FigureSI5.R
index ac1152c..9b24331 100644
--- a/rscripts/13_FigureSI5.R
+++ b/rscripts/13_FigureSI5.R
@@ -42,8 +42,8 @@ myload(dataset_9BSCBenthicPelagicGillnetSelectivity,
 ##----------------
 summary(dataset_9BSCBenthicPelagicGillnetSelectivity)
 dataset.thermal.trajectories <- dataset_9BSCBenthicPelagicGillnetSelectivity
-count.per.lake <- dataset.thermal.trajectories %>% group_by(lake.code) %>% 
-  summarise(total_count = n(),.groups = 'drop') %>%
+count.per.lake <- dataset.thermal.trajectories %>% dplyr::group_by(lake.code) %>% 
+  dplyr::summarise(total_count = n(),.groups = 'drop') %>%
   as.data.frame()
 
 
diff --git a/rscripts/14_SEManalysisSIDetailsIndirectEffects.R b/rscripts/14_SEManalysisSIDetailsIndirectEffects.R
index 05f4569..d90a5a5 100644
--- a/rscripts/14_SEManalysisSIDetailsIndirectEffects.R
+++ b/rscripts/14_SEManalysisSIDetailsIndirectEffects.R
@@ -13,6 +13,9 @@ library(semEff)
 library(stringr)
 library(tidyr)
 library(dplyr)
+library(tidyverse)
+library(ggbeeswarm)
+library(gghalves)
 
 ##FUNCTIONS##
 source("./rfunctions/misc.R")
@@ -176,13 +179,6 @@ tot.sp$response <- "Total species richness"
 DIE <- rbind(connectance, maxtl, css.slope, css.elevation, tot.sp)
 colnames(DIE) <- c('path', 'effect', 'Predictor', 'response')
 str(DIE)
-
-
-library(tidyverse)
-library(ggbeeswarm)
-library(gghalves)
-
-# default
 ggplot(DIE, aes(x = response, y = effect)) +
   geom_abline(slope = 0, intercept = 0, linetype = 2) +
   geom_half_boxplot(nudge = 0.05, outlier.color = NA, color = "black") +
@@ -193,4 +189,7 @@ ggplot(DIE, aes(x = response, y = effect)) +
   labs(y = "Standardized effet", x = "Response variables") +
   theme(axis.text = element_text(size = 14, colour = "#000000"),
         axis.title = element_text(size = 16, face = "bold", colour = "#000000"))
-#4 x 3 
\ No newline at end of file
+#4 x 3 
+
+#sem_details_indirect_effect_all_lakes<- DIE
+#mysave(sem_details_indirect_effect_all_lakes, dir = "outputs/SEM", overwrite = TRUE)