diff --git a/scripts/02_simulation/02_generate_environments.R b/scripts/02_simulation/02_generate_environments.R index 0322ccd..19df983 100644 --- a/scripts/02_simulation/02_generate_environments.R +++ b/scripts/02_simulation/02_generate_environments.R @@ -11,6 +11,12 @@ library(dplyr) library(ggplot2) +# load relevant functions +source("scripts/Function_plotting_theme.R") + +# get a colour palette +col_pal <- wesanderson::wes_palette(name = "Darjeeling1", 5) + # function to standardise a variable between two values range_stand <- function(x, min, max) { @@ -26,7 +32,7 @@ patches <- 5 timesteps <- 300 # the scale of environmental variation -env1Scale = 4 +env1Scale <- 5 # environmental variation for the spatial insurance effect @@ -51,11 +57,25 @@ for(j in 1:N_REP) { } -# check an example -ggplot(data = env_SI[[1]], - mapping = aes(x = time, y = env1, colour = as.character(patch) )) + - geom_line() - +# plot out 5 examples of the spatially fluctuating environment +p1 <- + bind_rows(env_SI[1], .id = "rep") %>% + as_tibble %>% + mutate(patch = as.character(patch), + id = paste(rep, patch, sep = "_") ) %>% + ggplot(data = ., + mapping = aes(x = time, y = env1, colour = patch, group = id )) + + geom_line(alpha = 1) + + labs(colour = "Patch") + + ylab("Environment (0-1)") + + xlab("Time") + + scale_colour_manual(values = col_pal) + + guides(colour = guide_legend(override.aes = list(alpha = 1))) + + scale_y_continuous(limits = c(-0.05, 1.05), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + + theme_meta() + + theme(legend.position = "top", + legend.key = element_rect(fill = NA)) +plot(p1) # environmental variation for the temporal insurance effect @@ -63,24 +83,43 @@ ggplot(data = env_SI[[1]], env_TI <- vector("list", length = N_REP) for(j in 1:N_REP) { + # get a fluctuating environment + env1 <- synchrony::phase.partnered(n = timesteps, gamma = env1Scale, mu = 0.5, sigma = 0.25)$timeseries[,1] + env_df <- data.frame() for(i in 1:patches){ - env1 <- synchrony::phase.partnered(n = timesteps, gamma = env1Scale, mu = 0.5, sigma = 0.25)$timeseries[,1] - env1 <- range_stand(env1, 0.1, 0.9) + env1 <- env1 + rnorm(n = 1, mean = 0, sd = 0.02) env_df <- rbind(env_df, data.frame(env1 = env1, patch = i, time = 1:timesteps)) } + env_df$env1 <- range_stand(env_df$env1, 0.1, 0.9) + # add to output list env_TI[[j]] <- env_df } -# plot an example -ggplot(data = env_TI[[1]], - mapping = aes(x = time, y = env1, colour = as.character(patch) )) + - geom_line() - +# plot out 5 examples of the temporally fluctuating environment +p2 <- + bind_rows(env_TI[1], .id = "rep") %>% + as_tibble %>% + mutate(patch = as.character(patch), + id = paste(rep, patch, sep = "_") ) %>% + ggplot(data = ., + mapping = aes(x = time, y = env1, colour = patch, group = id )) + + geom_line(alpha = 0.5) + + labs(colour = "Patch") + + ylab("") + + xlab("Time") + + scale_colour_manual(values = col_pal) + + guides(colour = guide_legend(override.aes = list(alpha = 1))) + + scale_y_continuous(limits = c(-0.05, 1.05), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + + theme_meta() + + theme(legend.position = "top", + legend.key = element_rect(fill = NA)) +plot(p2) + # environmental variation with spatial and temporal insurance @@ -101,10 +140,25 @@ for(j in 1:N_REP) { } -# plot an example -ggplot(data = env_com[[1]], - mapping = aes(x = time, y = env1, colour = as.character(patch) )) + - geom_line() +# plot out 5 examples of the temporally fluctuating environment +p3 <- + bind_rows(env_com[1], .id = "rep") %>% + as_tibble %>% + mutate(patch = as.character(patch), + id = paste(rep, patch, sep = "_") ) %>% + ggplot(data = ., + mapping = aes(x = time, y = env1, colour = patch, group = id )) + + geom_line(alpha = 0.5) + + labs(colour = "Patch") + + ylab("") + + xlab("Time") + + scale_colour_manual(values = col_pal) + + guides(colour = guide_legend(override.aes = list(alpha = 1))) + + scale_y_continuous(limits = c(-0.05, 1.05), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + + theme_meta() + + theme(legend.position = "top", + legend.key = element_rect(fill = NA)) +plot(p3) # quantify the variation in space and time @@ -148,6 +202,13 @@ bind_rows(env_com, .id = "rep") %>% # comparing the variances confirms the simulations +# make a plot for the supplementary +ggpubr::ggarrange(p1, p2, p3, + labels = c("a", "b", "c"), + font.label = list(size = 11, face = "plain"), + common.legend = TRUE, + nrow = 1) + # export these data for use in the simulations saveRDS(env_SI, "scripts/02_simulation/02_generated_environment_SI.rds") saveRDS(env_TI, "scripts/02_simulation/02_generated_environment_TI.rds") diff --git a/scripts/02_simulation/03_simulate_metacommunities.R b/scripts/02_simulation/03_simulate_metacommunities.R index 5f4c1df..12841fc 100644 --- a/scripts/02_simulation/03_simulate_metacommunities.R +++ b/scripts/02_simulation/03_simulate_metacommunities.R @@ -20,7 +20,7 @@ source("scripts/01_partition_functions/01_isbell_2018_partition.R") # load the simulated environments env_SI <- readRDS("scripts/02_simulation/02_generated_environment_SI.rds") -env_TI <- readRDS("scripts/02_simulation/02_generated_environment_SI.rds") +env_TI <- readRDS("scripts/02_simulation/02_generated_environment_TI.rds") env_com <- readRDS("scripts/02_simulation/02_generated_environment_com.rds") # pull the simulated environments into a list @@ -30,7 +30,7 @@ env_list <- list(env_SI, env_TI, env_com) set.seed(54258748) # number of replicate simulations -N_REP <- 500 +N_REP <- 10 assertthat::assert_that( N_REP == length(env_SI) && N_REP == length(env_SI) && N_REP == length(env_com) diff --git a/scripts/02_simulation/05_analyse_accuracy.R b/scripts/02_simulation/05_analyse_accuracy.R index 6557470..be298fa 100644 --- a/scripts/02_simulation/05_analyse_accuracy.R +++ b/scripts/02_simulation/05_analyse_accuracy.R @@ -26,17 +26,16 @@ output_rye <- output_df %>% filter( !(Beff %in% c( "LS", "LC", "TC", "NO" )) ) -# subset out the first 500 spatial insurance effects -output_rye <- - output_rye %>% - filter(sim_rep %in% 501:1000) +# add a variable describing the different levels of variation +output_rye$env_var <- rep(c("Spatial env. variation", "Temporal env. variation", "Combination"), each = 10*7) # get a table of the range of the different BEF effects output_rye %>% - group_by(Beff) %>% + group_by(env_var, Beff) %>% summarise(mean_BEF = mean(BEF_obs), min_BEF = min(BEF_obs), - max_BEF = max(BEF_obs)) + max_BEF = max(BEF_obs)) %>% + filter(Beff %in% c("NBE", "SI", "TI")) # does the observed value lie in the 90% percentile interval output_rye <- @@ -112,8 +111,4 @@ plot(p1) ggsave(filename = "figures/SI2_fig_S5.svg", p1, unit = "cm", width = 21, height = 21) - - - - ### END