update the simulation scripts

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,36 +57,69 @@ 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
 
 # replicate N_REP times
 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")

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)

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