cleaning up repo for publication. updated readme to walk through the …

…steps to recreate figures

mc_sim.R → 1_mc_sim.R

@@ -28,7 +28,7 @@ set.seed(82072)
 
 disp_rates <- 10^seq(-5, 0, length.out = 50)
 germ_fracs <- seq(.1,1, length.out = 10)
-surv_fracs <- c(.1, .5, 1)
+surv_fracs <- c(.5, .99)
 
 params <- expand.grid(disp_rates, germ_fracs, surv_fracs)
 
@@ -84,9 +84,7 @@ for(x in conditions){
                                                               env_niche_optima = "even")
 
                                disp_array <- generate_dispersal_matrices(landscape, species, patches, species_traits, torus = FALSE)
-                               # int_mat <- species_int_mat(species = species, intra = intra,
-                               #                            min_inter = min_inter, max_inter = max_inter,
-                               #                            comp_scaler = comp_scaler, plot = TRUE)
+
 
 
                                N <- init_community(initialization = initialization, species = species, patches = patches)
@@ -107,8 +105,6 @@ for(x in conditions){
                                  # compute r
                                  r <- compute_r_xt(species_traits, env = env, species = species)
 
-                                 # compute growth
-                                 #N_hat <- N*r / (1 + N%*%int_mat)
 
                                  # who germinates? Binomial distributed
                                  N_germ <- germination(N + D, species_traits)
@@ -131,15 +127,7 @@ for(x in conditions){
 
                                  dispSP <- colSums(E)
 
-                                 # determine immigrants to each patch
-                                 # I_hat_raw <- disp_array[,,1]%*%E
-                                 # I_hat <- t(t(I_hat_raw)/colSums(I_hat_raw))
-                                 # I_hat[is.nan(I_hat)] <- 1
-                                 # I <- sapply(1:species, function(x) {
-                                 #   if(dispSP[x]>0){
-                                 #     table(factor(sample(x = patches, size = dispSP[x], replace = TRUE, prob = I_hat[,x]), levels = 1:patches))
-                                 #   } else {rep(0, patches)}
-                                 # })
+
 
                                  I_hat_raw <- matrix(nrow = patches, ncol = species)
 
@@ -181,14 +169,7 @@ for(x in conditions){
                                                        dynamics_i)
                                }
 
-                               # every 20 tsteps?
-                               # dynamics_subset <- dynamics_out %>% 
-                               #   filter(time %in% seq(0, timesteps, by = 20))
-
-                               #dynamics_total <- rbind(dynamics_total, dynamics_subset)
-
-                               #saveRDS(dynamics_out, file = paste0("sim_output/sim_disp",disp,"_germ_",germ,"_surv_",surv,"_maxinter_",max_inter,"_mininter_",min_inter,".rds"))
-
+
 
                                return(dynamics_out)
                              }

mc_sim.sh → 1_mc_sim.sh

@@ -7,10 +7,10 @@
 #SBATCH --mem=124GB
 #SBATCH --mail-type=ALL
 #SBATCH [email protected]
-#SBATCH --job-name=metacom_99
+#SBATCH --job-name=metacom
 
 cd /gscratch/nwisnosk/GitHub/metacom-coexistence
 
 module load gcc/7.3.0 r/3.6.1
 
-R CMD BATCH --no-restore --no-save mc_sim_0.99.R mc_sim99.log
+R CMD BATCH --no-restore --no-save 1_mc_sim.R mc_sim.log

mc_sim_sensitivity.R → 2_mc_sim_sensitivity.R

@@ -92,9 +92,6 @@ for(e in 1:nenvs){
                                      env_niche_optima = "even")
 
       disp_array <- generate_dispersal_matrices(landscape, species, patches, species_traits, torus = FALSE)
-      # int_mat <- species_int_mat(species = species, intra = intra,
-      #                            min_inter = min_inter, max_inter = max_inter,
-      #                            comp_scaler = comp_scaler, plot = TRUE)
 
 
       N <- init_community(initialization = initialization, species = species, patches = patches)
@@ -115,9 +112,6 @@ for(e in 1:nenvs){
         # compute r
         r <- compute_r_xt(species_traits, env = env, species = species)
 
-        # compute growth
-        #N_hat <- N*r / (1 + N%*%int_mat)
-
         # who germinates? Binomial distributed
         N_germ <- germination(N + D, species_traits)
 
@@ -139,15 +133,6 @@ for(e in 1:nenvs){
 
         dispSP <- colSums(E)
 
-        # determine immigrants to each patch
-        # I_hat_raw <- disp_array[,,1]%*%E
-        # I_hat <- t(t(I_hat_raw)/colSums(I_hat_raw))
-        # I_hat[is.nan(I_hat)] <- 1
-        # I <- sapply(1:species, function(x) {
-        #   if(dispSP[x]>0){
-        #     table(factor(sample(x = patches, size = dispSP[x], replace = TRUE, prob = I_hat[,x]), levels = 1:patches))
-        #   } else {rep(0, patches)}
-        # })
 
         I_hat_raw <- matrix(nrow = patches, ncol = species)
 

mc_sim_sensitivity.sh → 2_mc_sim_sensitivity.sh

@@ -14,4 +14,4 @@ cd /gscratch/nwisnosk/GitHub/metacom-coexistence
 
 module load gcc/7.3.0 r/3.6.1
 
-R CMD BATCH --no-restore --no-save mc_sim_sensitivity.R mc_sim_sensitivity.log
+R CMD BATCH --no-restore --no-save 2_mc_sim_sensitivity.R mc_sim_sensitivity.log

3_make_diversity_plots.R

@@ -0,0 +1,130 @@
+library("tidyverse")
+
+theme_set(theme_light() + 
+            theme(strip.text = element_text(color = "black", size = 14),
+                  strip.background = element_rect(fill = "gray90"),
+                  axis.text = element_text(size = 14),
+                  axis.title = element_text(size = 16),
+                  axis.ticks.length = unit(x = .2, units = "cm"),
+                  legend.title = element_text(size = 14),
+                  legend.text = element_text(size = 12)))
+
+# competition scenarios
+equal_above <- "sim_output/final_equal_above.csv"
+stable_above <- "sim_output/final_stable_above.csv"
+equal_below <- "sim_output/final_equal_below.csv"
+stable_below <- "sim_output/final_stable_below.csv"
+
+
+file_list <- c(equal_above, stable_above, equal_below, stable_below)
+competition = file_list[1]
+make_plots <- function(competition){
+  condition_type <- str_remove(competition, ".csv") %>% 
+    str_remove("sim_output/final_") %>% 
+    str_split("_") %>% .[[1]] %>% .[1]
+
+  position <- str_remove(competition, ".csv") %>% 
+    str_remove("sim_output/final_") %>% 
+    str_split("_") %>% .[[1]] %>% .[2]
+
+
+  # read data and ignore s=0.1 condition
+  last_t_out <- read_csv(competition) %>% 
+    filter(survival %in% c(0.5, 0.99))
+
+
+  disp_rates <- sort(unique(last_t_out$dispersal))
+  germ_fracs <- sort(unique(last_t_out$germination))
+  surv_fracs <- sort(unique(last_t_out$survival))
+
+  alpha_out <- last_t_out %>%
+    group_by(dispersal, germination, survival, patch, rep) %>% 
+    mutate(N_pa = 1*(N > 0)) %>%
+    summarize(alpha = sum(N_pa)) %>%
+    ungroup()
+
+  alpha_out <- alpha_out %>%
+    group_by(dispersal, germination, survival, rep) %>%
+    summarize(mean_alpha = mean(alpha))
+
+  gamma_out <- last_t_out %>%
+    ungroup() %>%
+    mutate(N_pa = 1*(N > 0)) %>%
+    group_by(dispersal, germination, survival, rep) %>%
+    summarize(gamma = length(unique(species)))
+
+  div_part <- alpha_out %>%
+    left_join(gamma_out) %>%
+    mutate(beta = gamma / mean_alpha)
+
+  no_sb <- div_part %>% 
+    filter(germination == 1) %>% 
+    mutate(survival = "No seed bank")
+
+  div_part <- div_part %>% 
+    mutate(
+      survival = case_when(
+        survival == 0.5 ~ "Survival = 0.5",
+        survival == 0.99 ~ "Survival = 0.99"
+      )
+    ) %>% bind_rows(no_sb)
+
+
+  germ_plot <- div_part %>%
+
+    gather(mean_alpha, beta, gamma, key = "partition", value = "diversity") %>%
+    mutate(partition = factor(partition, levels = c("mean_alpha", "beta", "gamma"),
+                              labels = c("mean alpha", "beta", "gamma"))) %>%
+    ggplot(aes(x = dispersal, y = diversity, color = germination, fill = germination, group = germination)) +
+    geom_point(alpha = 0.05) +
+    #geom_line(mapping = aes(group = paste(germination, rep)), stat = "smooth", alpha = 0.2, se = F) +
+    geom_smooth(alpha = 0.25) +
+    facet_grid(rows = vars(partition), cols = vars(survival), drop = FALSE, scales = "free_y") +
+    scale_x_log10() +
+    scale_y_continuous(breaks = scales::pretty_breaks(4)) +
+    scale_color_viridis_c("Germination", option = "B", end = .8) +
+    scale_fill_viridis_c("Germination", option = "B", end = .8) +
+    labs(y = "Diversity", x = "Dispersal") +
+    theme(legend.position = "right")
+  ggsave(filename = paste0("figures/germination_", position, "_", condition_type,".pdf"),
+         plot = germ_plot, width = 12, height = 8)
+
+  surv_plot <- div_part %>%
+    filter(germination <= 0.5) %>% 
+    mutate(germination = factor(germination, levels = germ_fracs, labels = paste("Germ. =", round(germ_fracs,2)))) %>%
+    gather(mean_alpha, beta, gamma, key = "partition", value = "diversity") %>%
+    mutate(partition = factor(partition, levels = c("mean_alpha", "beta", "gamma"),
+                              labels = c("mean alpha", "beta", "gamma"))) %>%
+    ggplot(aes(x = dispersal, y = diversity, color = survival, fill = survival, group = survival)) +
+    geom_point(alpha = 0.05) +
+    #geom_line(mapping = aes(group = paste(survival, rep)), stat = "smooth", alpha = 0.2, se = F) +
+    geom_smooth(alpha = 0.25) +
+    facet_grid(rows = vars(partition), cols = vars(germination), drop = TRUE, scales = "free_y") +
+    scale_x_log10() +
+    scale_y_continuous(breaks = scales::pretty_breaks(4)) +
+    scale_color_viridis_d("Survival", option = "D", end = .9) +
+    scale_fill_viridis_d("Survival", option = "D", end = .9) +
+    labs(y = "Diversity", x = "Dispersal") +
+    theme(legend.position = "right") 
+  ggsave(filename = paste0("figures/survival_", position, "_", condition_type,".pdf"), 
+         plot = surv_plot, width = 12, height = 8)
+
+  germ_comp <- div_part %>% 
+    filter(germination == 1 | germination == 0.1) %>% 
+    filter(survival == 1) %>% 
+    mutate(germination = factor(germination, levels = germ_fracs, labels = paste("Germ. =", round(germ_fracs,2))),
+           survival = factor(survival, levels = surv_fracs, labels = paste("Surv. =", round(surv_fracs,2)))) %>%
+    gather(mean_alpha, beta, gamma, key = "partition", value = "diversity") %>%
+    mutate(partition = factor(partition, levels = c("mean_alpha", "beta", "gamma"), 
+                              labels = c("mean alpha", "beta", "gamma"))) 
+
+
+}
+
+### Loop through files
+for(competition in file_list){
+  make_plots(competition)
+}
+
+
+

sensitivity_analysis.R → 4_plot_sensitivity_analysis.R

@@ -2,10 +2,7 @@ library("tidyverse")
 library("broom")
 library("colorspace")
 
-normalize = function(x){
-  normalized = (x-min(x))/(max(x)-min(x))
-  return(normalized)
-}
+
 
 theme_set(theme_light() + 
             theme(strip.text = element_text(color = "black", size = 14),
@@ -128,16 +125,7 @@ x_axis_labs <- c("dispersal",
                  "germination : survival", "dispersal : germination : survival")
 
 int_sensitivity_plot <- sens_out %>% 
-  # mutate(term = factor(term, levels = c(
-  #   "Intercept", "dispersal", "germination", "survival",
-  #   "dispersal:germination", "dispersal:survival", "germination:survival",
-  #   "dispersal:germination:survival"
-  # ), labels = c(
-  #   "intercept", "dispersal", "germination", "survival",
-  #   "dispersal x \n germination", "dispersal x \n survival", "germination x \n survival",
-  #   "dispersal x \n germination x \n survival"
-  # ))) %>%
-
+
   ggplot(aes(y = estimate, ymin = (estimate - std.error), ymax = (estimate + std.error),
              color = response, fill = response, x = term)) +
   geom_vline(xintercept = c(1:9+0.5), alpha = 0.5, color = "gray50", size = 0.2) +

mc_sensitivity.sh → 4_plot_sensitivity_analysis.sh

@@ -14,4 +14,4 @@ cd /gscratch/nwisnosk/GitHub/metacom-coexistence
 
 module load gcc/7.3.0 r/3.6.1
 
-R CMD BATCH --no-restore --no-save sensitivity_analysis.R sensitivity_analysis.log
+R CMD BATCH --no-restore --no-save 4_plot_sensitivity_analysis.R plot_sensitivity_analysis.log

README.md

@@ -1,15 +1,17 @@
-# metacom-coexistence
-Examining coexistence in metacommunities with dormancy.
+Examining diversity in metacommunities with dormant seed banks.
 
 This repository contains the code to reproduce the manuscript: 
 
-Wisnoski NI, and LG Shoemaker. "Seed banks alter metacommunity diversity: the interactive effects of competition, germination, and survival"
+**Wisnoski NI, and LG Shoemaker. "Seed banks alter metacommunity diversity: the interactive effects of competition, dispersal, and dormancy"**
 
 To reproduce the analysis and figures from the manuscript, run the following:
 
-1. `mc_sim.R` = this file includes the main simulation across parameter space in the paper.
-2. `sim_output/extract_sim_output.sh` = commands included here parse the simulation output and split it into files for above and belowground data for equal and stable interactions. Paths and filenames here may be user-specific and likely need tweaking for your system.
-3.  `make_plots.R` = this file makes figures from the output of the simulations
-4.  `mc_sim_tradeoff.R` = simulation for condition where there's interspecific trait variation
-5.  `analyze_tradeoff.R` = making tradeoff related figures.
+1. `1_mc_sim.R` = this file includes the main simulation across parameter space in the paper. It generates output files that are used in figures 2 and 3.
+2. `sim_output/extract_sim_output.sh` = commands included here parse the simulation output and split it into files for above and belowground data for equal and stable interactions. Paths and filenames will be user-specific and need to be specified for your system.
+3.  `2_mc_sim_sensitivity.R` = this file includes the simulation model for the sensitivity analysis. It generates output files that are used to make figure 4. 
+4.  `3_make_diversity_plots.R` = this file makes figures 2 and 3 from the output of the simulations reformated by the `extract_sim_output.sh` script.
+5.  `4_plot_sensitivity_analysis.R` = this script creates figure 4 from the output of the sensitivity analysis simulation script.
+6.  `s1_mc_sim_tradeoff.R` = simulation model for the scenario where species exhibit interspecific trait variation in dispersal, germination, and survival. 
+7.  `s2_analyze_tradeoff.R` = this script makes the tradeoff related figures included in the supplement.
+8.  `s3_visualize_example_tradeoffs.R` = this script demonstrates the differences in none, weak, and strong tradeoff scenarios explored in the supplement.