Plotting the simulation results

Load packages and results

packages <- c("dplyr", "purrr", "ggplot2", "reshape2", "Cairo", "knitr",  
              "latex2exp", "pander", "grid", "gridExtra", "ggthemes", "data.table",
              "readr", "tibble", "biglm", "kableExtra", "future", "future.apply")
shh <- suppressMessages(lapply(packages, library, character.only = TRUE, quietly = TRUE))

if(!file.exists("data/all_results.rds")){
  
  # Open all the results files and bind them together into a tibble
  results <- lapply(list.files(path = file.path(getwd(), "data"), pattern = "results_", full.names = TRUE), 
                    function(focal_file) {
                      readRDS(focal_file) %>% filter(cost_A == 0 & cost_B == 0)}) %>% 
    rbindlist() %>% as_tibble()
  
  # Make sure each unique parameter space is run exactly once
  results <- results %>% 
    distinct(release_strategy, W_shredding_rate, Z_conversion_rate,
             Zr_creation_rate, Zr_mutation_rate, Wr_mutation_rate,
             cost_Zdrive_female, cost_Zdrive_male, male_migration_prob,
             female_migration_prob, migration_type, n_patches, softness,
             male_weighting, density_dependence_shape, max_fecundity, realisations,
             initial_A, initial_B, .keep_all = TRUE)
  
  saveRDS(results, file = "data/all_results.rds")
} else results <- read_rds("data/all_results.rds")

results$went_extinct <- ifelse(results$outcome == "Population extinct", 1, 0)
results$migration_type[results$migration_type == "local"] <- "Local"
results$migration_type[results$migration_type == "global"] <- "Global"

# Command to upload the results file to Spartan - needed for Spartan to queue up all the parameter spaces not already finished
# scp /Users/lholman/Rprojects/W_shredder/data/all_results.rds lukeholman@spartan:/data/projects/punim0243/W_shredder/data/all_results.rds

# Helper function to get a list of all the model parameters that vary between runs
find_variable_parameters <- function(dat){
  dat %>% 
    select(-id, -realisations, -generation_extinct, -generation_Zd_extinct, 
           -generation_W_extinct, -generation_Zd_fixed, -outcome, -went_extinct, 
           -mating_table, -initial_Wr, -initial_Zr) %>%
    sapply(function(x) length(unique(x))) %>% 
    keep(~.x > 1) %>% names()
}
variable_parameters <- find_variable_parameters(results)
combinations <- apply(combn(variable_parameters, 2), 2, paste0, collapse = " x ") # all 2-way combos of parameters

# Make a data frame to convert R-friendly names to figure-friendly names
nice_names <- data.frame(
  original = c(variable_parameters, combinations),
  new = gsub("_", " ", c(variable_parameters, combinations)),
  stringsAsFactors = FALSE) %>%
  mutate(
    new = gsub("rel", "Rel", new),
    new = gsub("W shredding rate", "Gene drive in females ($p_{shred}$)", new),
    new = gsub("Z conversion rate", "Gene drive in males ($p_{conv}$)", new),
    new = gsub("Zr creation rate", "Rate of NHEJ ($p_{nhej}$)", new),
    new = gsub("Zr mutation rate", "Z-linked resistance ($\\mu{_Z}$)", new),
    new = gsub("Wr mutation rate", "W-linked resistance ($\\mu{_W}$)", new),
    new = gsub("cost Zdrive female", "Cost of Z* to females ($c_f$)", new),
    new = gsub("cost Zdrive male", "Cost of Z* to males ($c_m$)", new),
    new = gsub("male migration prob", "Male dispersal frequency", new),
    new = gsub("feMale dispersal frequency", "Female dispersal frequency (xf)", new),
    new = gsub("Male dispersal frequency", "Male dispersal frequency (xm)", new),
    new = gsub("xm", "$x_m$", new), 
    new = gsub("xf", "$x_f$", new),
    new = gsub("migration type", "Dispersal type", new),
    new = gsub("n patches", "Number of patches ($k$)", new),
    new = gsub("softness", "Importance of local density ($\\psi$)", new),
    new = gsub("male weighting", "Male effect on density ($\\delta$)", new),
    new = gsub("density dependence shape", "Shape of density dependence ($\\alpha$)", new),
    new = gsub("max fecundity", "Maximum fecundity ($r$)", new),
    new = gsub("initial A", "Initial freq. W-shredding resistance allele A", new),
    new = gsub("initial B", "Initial freq. gene conversion resistance allele B", new),
    new = as.character(TeX(new))) %>% 
  mutate(new = gsub("mu", "\\mu", new))

# Sort the results between those that model a W-shredder, versus those that model a Z* that sterilises females
W_shredder <- results %>% filter(cost_Zdrive_female != 1)
female_sterilising <- results %>% filter(cost_Zdrive_female == 1)

pretty <- function(x) prettyNum(x, big.mark = ",", scientific = FALSE)

# Count the number of runs and save for entry into the manuscript
data.frame(total_runs = pretty(nrow(results)),
           W_shredder_runs = pretty(nrow(W_shredder)),
           female_sterilising_runs = pretty(nrow(female_sterilising)),
           W_shredder_runs_extinct = pretty(sum(W_shredder$went_extinct)),
           female_sterilising_runs_extinct = pretty(sum(female_sterilising$went_extinct)))  %>%
  saveRDS(file = "data/number_of_runs.rds")
  
rm(results)

Table showing the frequencies of each possible outcome

Table S1: The number and percentage of simulation runs that ended with the five possible outcomes, for the subset of simulation runs focusing on a W-shredder gene drive.

make_tally_table <- function(dat){
  outcomes <- dat$outcome %>% table() %>% melt() %>% arrange(-value) %>% mutate(p = round(100 * value / sum(value), 1))
  names(outcomes) <- c("Outcome", "Number of simulations", "%")
  outcomes %>% mutate(Outcome = gsub("without", "without causing", Outcome),
                      Outcome = gsub("n extinct", "n went extinct", Outcome),
                      Outcome = gsub("d extinct", "d went extinct", Outcome),
                      Outcome = gsub("Zd", "Z*", Outcome)) 
}

table_S1 <- make_tally_table(W_shredder) 
saveRDS(table_S1, file = "figures/tableS1.rds")
table_S1 %>% pander()

Table S2: The number and percentage of simulation runs that ended with the five possible outcomes, for the subset of simulation runs focusing on a female-sterilising Z-linked gene drive.

table_S2 <- make_tally_table(female_sterilising)
saveRDS(table_S2, file = "figures/tableS2.rds")
table_S2 %>% pander()

Make Figure 2

plot_run <- function(){
  
  get_data <- function(model_id, allele_freqs_list, result_df, label){
    
    df <- pluck_allele_freqs(model_id, allele_freqs_list) 
    
    alleles_to_plot <- group_by(df, allele) %>% 
      summarise(uniques = length(unique(frequency))) %>% 
      filter(uniques > 1) %>% # Don't plot alleles that stay at 0 whole time
      pull(allele)
    
    alleles_to_plot <- c(alleles_to_plot, "N")
    df <- df[df$allele %in% alleles_to_plot, ]
    df$frequency[df$allele == "N"] <- df$frequency[df$allele == "N"] /
      max(df$frequency[df$allele == "N"])
     
    paras <- result_df %>% mutate(id = as.character(id)) %>% filter(id == model_id)
    last_generation <- tail(df,12) %>% select(allele, frequency)
    print(paras); print(last_generation)

    df %>% mutate(facet = label,
                  allele = replace(allele, allele == "Zd", "Z*"),
                  allele = replace(allele, allele == "females", "Females"))
  }
  
  allele_freqs <- read_rds("data/allele_freqs_1.rds")
  
  rbind(
    get_data("10024119427447", allele_freqs, W_shredder, "A. Z* causes rapid extinction"),
    get_data("1002596940286", allele_freqs, W_shredder, "B. Z* fails to cause extinction"),
    get_data("10024112958792", allele_freqs, W_shredder, "C. Resistance to Z* prevents extinction")) %>%
    ggplot(aes(x = generation, y = frequency, colour = allele, group = allele)) + 
    geom_vline(xintercept = 50, linetype = 3, colour = "grey10", size = 0.9) + 
    geom_line(size = 0.9, alpha = 0.9) + 
    facet_wrap(~facet, scales = "free_x") +
    theme_hc() + scale_colour_hc(name = "") +
    theme(strip.text = element_text(hjust = 0), 
          strip.background = element_rect(fill = "grey90"),
          legend.position = "top",
          axis.ticks.y = element_blank()) + 
    xlab("Generation") + ylab("Frequency")
}

fig2 <- plot_run()

fig2

Past versions of unnamed-chunk-5-1.png

Figure 2: Three illustrative runs of the simulation, showing evolution in response to the introduction of 20 males carrying a W-shredder at Generation 50 (marked by the dotted line). In panel A, the driving Z* allele fixed very quickly, causing population extinction as the number of females dropped to zero. In panel B, the Z* allele spread up until the point that its fitness costs began to negate its transmission advantage, causing the population to halve in size but not to go extinct. In panel C, the Z* allele invaded, selecting for the resistance alleles A and Zr, and causing the Z* allele to reverse course and go extinct. The population size N is shown as a fraction of its maximum value of 10,000. Table S3 gives the parameter spaces used for these three runs.

tabl <- t(rbind(
  data.frame(Panel = "A", W_shredder %>% filter(id == "10024119427447")),
  data.frame(Panel = "B", W_shredder %>% filter(id == "1002596940286")),
  data.frame(Panel = "C", W_shredder %>% filter(id == "10024112958792"))) %>%
  select(release_strategy, W_shredding_rate, 
         Z_conversion_rate, Zr_creation_rate, Zr_mutation_rate,
         Wr_mutation_rate, cost_Zdrive_female, cost_Zdrive_male,
         male_migration_prob, female_migration_prob,
         migration_type, n_patches, softness, male_weighting,
         density_dependence_shape, max_fecundity, initial_A, initial_B))

nice_names2 <- data.frame(original = variable_parameters,
                          new = gsub("_", " ", variable_parameters),
                          stringsAsFactors = FALSE) %>%
  mutate(
    new = gsub("rel", "Rel", new),
    new = gsub("W shredding rate", "Gene drive in females (p_shred)", new),
    new = gsub("Z conversion rate", "Gene drive in males (p_conv)", new),
    new = gsub("Zr creation rate", "Rate of NHEJ (p_nhej)", new),
    new = gsub("Zr mutation rate", "Z-linked resistance (μ_Z)", new),
    new = gsub("Wr mutation rate", "W-linked resistance (μ_W)", new),
    new = gsub("cost Zdrive female", "Cost of Z* to females (c_f)", new),
    new = gsub("cost Zdrive male", "Cost of Z* to males (c_m)", new),
    new = gsub("male migration prob", "Male dispersal frequency", new),
    new = gsub("feMale dispersal frequency", "Female dispersal frequency (x_f)", new),
    new = gsub("Male dispersal frequency", "Male dispersal frequency (x_m)", new),
    new = gsub("migration type", "Dispersal type", new),
    new = gsub("n patches", "Number of patches (k)", new),
    new = gsub("softness", "Importance of local density (ψ)", new),
    new = gsub("male weighting", "Male effect on density (𝛿)", new),
    new = gsub("density dependence shape", "Shape of density dependence (α)", new),
    new = gsub("max fecundity", "Maximum fecundity (r)", new),
    new = gsub("initial A", "Intial freq. W-shredding resistance allele A", new),
    new = gsub("initial B", "Intial freq. gene conversion resistance allele B", new)) 

rownames(tabl) <- nice_names2$new[match(rownames(tabl), nice_names2$original)]
rownames(tabl) <- gsub("[.]", "", rownames(tabl))
colnames(tabl) <- c("Panel A", "Panel B", "Panel C")

Table S3: List of the parameter values used to generate the simulation runs shown in Figure 1.

as.data.frame(tabl) %>% kable(format = "html", escape = FALSE) %>% kable_styling()

Plot the effects of each parameter on extinction

# Find the proportion of runs that went extinct for each parameter value
get_percent_extinct <- function(dat, parameter){
  dat %>%
    group_by(!! sym(parameter)) %>%
    summarise(extinct = sum(went_extinct),
              not_extinct = n() - sum(went_extinct),
              prob = list(binom.test(extinct, extinct + not_extinct))) %>%
    rowwise() %>%
    mutate(response_variable = 100 * extinct/(extinct + not_extinct),
           lower_95_CI = 100 * prob$conf.int[1],
           upper_95_CI = 100 * prob$conf.int[2],
           parameter =  parameter) %>% rename(value = !! sym(parameter)) %>%
    select(parameter, value, everything()) %>% select(-prob) %>% ungroup()
} 

# Find the median time to extinction, among just those runs in which extinction occurred, for each parameter value
get_extinction_time <- function(dat, parameter){
  dat %>%
    filter(!is.na(generation_extinct)) %>% # only include runs where extinction occurred
    group_by(!! sym(parameter)) %>%
    summarise(response_variable = median(generation_extinct),
              lower_95_CI = quantile(generation_extinct, probs = 0.025),
              upper_95_CI = quantile(generation_extinct, probs = 0.975),
              n = n(),
              parameter =  parameter) %>% rename(value = !! sym(parameter)) %>%
    select(parameter, value, everything()) %>% ungroup()
} 

# To calculate separately for the 3 levels of Z_conversion_rate (male gene drive)
split_by_male_drive <- function(dat, fun){
  rbind(
    lapply(
      variable_parameters[variable_parameters != "Z_conversion_rate"], 
      get(fun), 
      dat = dat %>%
        filter(Z_conversion_rate == 0)) %>% do.call("rbind", .) %>% mutate(Z_conversion_rate = 0),
    lapply(
      variable_parameters[variable_parameters != "Z_conversion_rate"], 
      get(fun), 
      dat = dat %>%
        filter(Z_conversion_rate == 0.5)) %>% do.call("rbind", .) %>% mutate(Z_conversion_rate = 0.5),
    lapply(
      variable_parameters[variable_parameters != "Z_conversion_rate"], 
      get(fun), 
      dat = dat %>%
        filter(Z_conversion_rate == 0.95)) %>% do.call("rbind", .) %>% mutate(Z_conversion_rate = 0.95)
  ) %>% mutate(Z_conversion_rate = factor(Z_conversion_rate, levels = c(0, 0.5, 0.95)))
}

# % extinct, calculated separately for the three levels of Z_conversion_rate
W_shredder_split_by_male_drive_percent_extinct <- split_by_male_drive(W_shredder, "get_percent_extinct")
female_sterilising_split_by_male_drive_percent_extinct <- split_by_male_drive(female_sterilising, "get_percent_extinct") %>%
  filter(!(parameter %in% c("cost_Zdrive_female", "W_shredding_rate")))

# Time until extinction, calculated separately for the three levels of Z_conversion_rate
W_shredder_split_by_male_drive_time_to_extinct <- split_by_male_drive(W_shredder, "get_extinction_time")
female_sterilising_split_by_male_drive_time_to_extinct <- split_by_male_drive(female_sterilising, "get_extinction_time") %>%
  filter(!(parameter %in% c("cost_Zdrive_female", "W_shredding_rate")))

# Function to make Figure 3 (and similar supplementary figures)
make_multipanel_figure <- function(dat, ylab){
  parameter_importance <- dat %>%
    mutate(parameter = nice_names$new[match(parameter, nice_names$original)]) %>%
    group_by(parameter) %>% summarise(range = max(response_variable) - min(response_variable)) %>%
    arrange(-range) %>% pull(parameter)
  
  levels <- levels(factor(dat$value))
  levels <- c(levels[!(levels %in% c("all_patches", 100))], 100, "all_patches")
  levels <- replace(levels, levels == "all_patches", "Release in\nall patches")
  levels <- replace(levels, levels == "one_patch", "Release in\na single patch")
  
  pd <- position_dodge(0.2)

  dat %>% 
    mutate(parameter = nice_names$new[match(parameter, nice_names$original)],
           parameter = factor(parameter, parameter_importance),
           value = replace(value, value == "all_patches", "Release in\nall patches"),
           value = replace(value, value == "one_patch", "Release in\na single patch"),
           value = factor(value, levels)) %>%
    ggplot(aes(value, response_variable, colour = Z_conversion_rate)) + 
    geom_errorbar(aes(ymin = lower_95_CI, ymax = upper_95_CI), width = 0, position = pd) +
    geom_point(position = pd) + 
    scale_colour_brewer(palette = "Set2", name = "Gene drive in males\n(p_conv)") + 
    xlab("Parameter value") +
    ylab(ylab) + 
    theme_bw() + 
    theme(strip.background = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.border = element_rect(fill = NA, size = 1))
  
}

fig3 <- make_multipanel_figure(W_shredder_split_by_male_drive_percent_extinct, 
                      ylab = "% simulation runs resulting in extinction (\u00B1 95% CIs)")

cairo_pdf(file = "figures/figure3.pdf", width = 9, height = 12)
fig3 + facet_wrap(~parameter, scales = "free", labeller = label_parsed, ncol = 4) + theme(legend.position = c(.9, .1))
dev.off()

quartz_off_screen 
                2 

fig3 + facet_wrap(~parameter, scales = "free", labeller = label_parsed, ncol = 3) 

S1_fig <- make_multipanel_figure(filter(W_shredder_split_by_male_drive_time_to_extinct, n >= 40), 
                        ylab = "Median generations until extinction (\u00B1 95% quantiles)") + 
  facet_wrap(~parameter, scales = "free", labeller = label_parsed, ncol = 3) + 
  coord_cartesian(ylim = c(0, 150))

ggsave(S1_fig, filename = "figures/S1_fig.pdf", height = 12, width = 19.5)
saveRDS(S1_fig, file = "figures/S1_fig.rds")
S1_fig

S3_fig <- make_multipanel_figure(female_sterilising_split_by_male_drive_percent_extinct,
                        ylab = "% simulation runs resulting in extinction (\u00B1 95% CIs)") + 
  facet_wrap(~parameter, scales = "free", labeller = label_parsed, ncol = 4)

ggsave(S3_fig, filename = "figures/S3_fig.pdf", height = 9.6, width = 10)
saveRDS(S3_fig, file = "figures/S3_fig.rds")
S3_fig

S4_fig <- make_multipanel_figure(filter(female_sterilising_split_by_male_drive_time_to_extinct, n >= 40), 
                        ylab = "Median generations until extinction (\u00B1 95% quantiles)") + 
  facet_wrap(~parameter, scales = "free", labeller = label_parsed, ncol = 3) + 
  coord_cartesian(ylim = c(0, 100))

ggsave(S4_fig, filename = "figures/S4_fig.pdf", height = 12, width = 10)
saveRDS(S4_fig, file = "figures/S4_fig.rds")
S4_fig

Finding variables with interacting effects on the extinction probability

run_big_glm <- function(dat, sterilising = FALSE){
  
  if(sterilising){
    variable_parameters <- variable_parameters[!(variable_parameters %in% c("W_shredding_rate", "cost_Zdrive_female"))] 
  }
  
  formula <- paste("went_extinct ~ (",
                   paste0(variable_parameters, collapse = " + "), ")^2",
                   sep = "")
  
  if(!sterilising) formula <- paste(formula, "- W_shredding_rate:cost_Zdrive_female")
  
  bigglm(as.formula(formula), 
         data = dat %>% mutate_if(is.numeric, function(x) as.numeric(scale(x))), 
         chunksize = 10000)
}

# Helper to get the fixed effects table out of the biglm object (no 'coefficients' slot, apparently)
tidy_model <- function(model){
  
  out <- capture.output(model %>% summary)[-(1:5)] %>% lapply(function(x) {
    x <- strsplit(x, split = " ")[[1]]
    x[x != ""]
    }) %>% do.call("rbind", .) %>% as.data.frame(stringsAsFactors = FALSE) %>% as_tibble() %>% 
    mutate(V1 = as.factor(gsub(":", " x ", V1))) %>% mutate_if(is.character, as.numeric) 
  
  names(out) <- c("Parameter", "Estimate", "Lower_95_CI", "Upper_95_CI", "SE", "p")
  out %>% arrange(-abs(Estimate))
}

if(!file.exists("data/W_shredder_model.rds")) {
  W_shredder_model <- run_big_glm(W_shredder)
  saveRDS(tidy_model(W_shredder_model), "data/W_shredder_model.rds")
  female_sterilising_model <- run_big_glm(female_sterilising, sterilising = TRUE)
  saveRDS(tidy_model(female_sterilising_model), "data/female_sterilising_model.rds")
} else {
  W_shredder_model <- readRDS("data/W_shredder_model.rds")
  female_sterilising_model <- readRDS("data/female_sterilising_model.rds")
}

importance_plot <- function(dat, n = "all"){
  
  if(n != "all"){
    dat <- dat %>%
      mutate(row=1:n()) %>% filter(row %in% 1:n)
  } 
  
  for(i in 1:nrow(dat)){
    if(dat$Estimate[i] < 0){
      lower <- dat$Lower_95_CI[i]
      dat$Lower_95_CI[i] <- dat$Upper_95_CI[i]
      dat$Upper_95_CI[i] <- lower
    }
  }
  
  dat <- dat %>%
    mutate(Parameter = gsub("one_patch", "", as.character(Parameter)),
           Parameter = gsub("Local", "", Parameter),
           Parameter = factor(Parameter, rev(Parameter)),
           tick_label = nice_names$new[match(Parameter, nice_names$original)]) 
  
  dat %>%
    ggplot(aes(Parameter, abs(Estimate))) + 
    geom_bar(stat = "identity", fill = "tomato") +
    geom_errorbar(aes(ymin = abs(Lower_95_CI), ymax = abs(Upper_95_CI)), width = 0) + 
    coord_flip() + 
    scale_y_continuous(expand = c(0, 0)) + 
    scale_x_discrete(labels = parse(text = rev(dat$tick_label))) + 
    xlab("Predictor variable") + ylab("Absolute effect size (\u00B1 95% CIs)") + 
    theme_classic() + 
    theme(axis.ticks.y = element_blank())
}

W_shred        <- importance_plot(W_shredder_model, n = 25)
female_sterile <- importance_plot(female_sterilising_model, n = 25)

W_shred

female_sterile

get_percent_extinct_double <- function(dat, p1, p2){
  output <-  dat %>%
    group_by(!! sym(p1), !! sym(p2)) %>%
    summarise(extinct = sum(went_extinct),
              not_extinct = n() - sum(went_extinct),
              prob = list(binom.test(extinct, extinct + not_extinct))) %>%
    rowwise() %>%
    mutate(percent_extinct = 100 * extinct / (extinct + not_extinct),
           lower_95_CI = 100 * prob$conf.int[1],
           upper_95_CI = 100 * prob$conf.int[2],
           parameter_1 =  p1, 
           parameter_2 =  p2) 
  names(output)[1:2] <- paste("value", 1:2, sep = "_")
  output %>%
    mutate(value_1 = as.character(value_1), value_2 = as.character(value_2)) %>%
    select(parameter_1, parameter_2, value_1, value_2, everything()) %>% 
    select(-prob) %>% ungroup()
} 

top_12_shredder    <- (W_shredder_model %>% mutate(Parameter = as.character(Parameter)) %>% filter(grepl(" x ", Parameter)) %>% pull(Parameter))[1:12]
top_12_sterilising <- (female_sterilising_model %>% mutate(Parameter = as.character(Parameter)) %>% filter(grepl(" x ", Parameter)) %>% pull(Parameter))[1:12]

two_way_W_shredder <- map2_df(str_split(top_12_shredder, " x ", simplify = TRUE)[, 1], 
                              str_split(top_12_shredder, " x ", simplify = TRUE)[, 2], 
                              get_percent_extinct_double, dat = W_shredder)

two_way_sterilising <- map2_df(str_split(top_12_sterilising, " x ", simplify = TRUE)[, 1], 
                               str_split(top_12_sterilising, " x ", simplify = TRUE)[, 2], 
                               get_percent_extinct_double, dat = W_shredder)



two_way_plot <- function(two_way_data, pairwise_combos){
  
  two_way_data <- two_way_data %>%
    mutate(pasted = paste(parameter_1, parameter_2, sep = " x "),
           parameter_1 = nice_names$new[match(parameter_1, nice_names$original)],
           parameter_2 = nice_names$new[match(parameter_2, nice_names$original)]) %>%
    filter(pasted %in% pairwise_combos) 
  z_range <- range(two_way_data$percent_extinct) * 1.1
  
  make_one <- function(df){
    
    df$value_1 <- factor(df$value_1, unique(df$value_1))
    df$value_2 <- factor(df$value_2, unique(df$value_2))
    
    ggplot(df, aes(value_1, value_2, fill = percent_extinct)) + 
      geom_tile(colour = "black", linetype = 3) + 
      scale_fill_distiller(palette = "PuBuGn", direction = 1, limits = z_range, name = "% extinct") + 
      theme_bw() + 
      theme() + 
      scale_x_discrete(expand = c(0,0), name = parse(text = df$parameter_1[1])) + 
      scale_y_discrete(expand = c(0,0), name = parse(text = df$parameter_2[1])) 
      
  }
  
  grid_arrange_shared_legend <- function(...) {
    plots <- list(...)
    g <- ggplotGrob(plots[[1]] + theme(legend.position="bottom"))$grobs
    legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
    lheight <- sum(legend$height)
    plots <- lapply(plots, function(x) ggplotGrob(x + theme(legend.position = "none")))
    
    p <- gtable_cbind(gtable_rbind(plots[[1]], plots[[4]], plots[[7]]),
                      gtable_rbind(plots[[2]], plots[[5]], plots[[8]]),
                      gtable_rbind(plots[[3]], plots[[6]], plots[[9]]))
    
    grid.arrange(p, legend, 
                 nrow = 2, 
                 heights = unit.c(unit(1, "npc") - lheight, lheight))
  }

  dat_list <- vector(mode = "list", length = length(pairwise_combos))
  for(i in 1:length(pairwise_combos)){
    dat_list[[i]] <- two_way_data %>% filter(pasted == pairwise_combos[i])
  }
  
  lapply(dat_list, make_one) %>% do.call("grid_arrange_shared_legend", .)
}

two_way_plot_shredder <- two_way_plot(two_way_W_shredder, top_12_shredder)

two_way_plot_sterilising <- two_way_plot(two_way_sterilising, top_12_sterilising)
saveRDS(two_way_plot_shredder, file = "figures/S2_fig.rds")
saveRDS(two_way_plot_sterilising, file = "figures/S6_fig.rds")
grid.draw(two_way_plot_shredder)

grid.draw(two_way_plot_sterilising)

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] grid      parallel  stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] future.apply_1.0.1 future_1.11.1.1    kableExtra_0.9.0  
 [4] biglm_0.9-1        DBI_1.0.0          data.table_1.12.2 
 [7] ggthemes_4.0.1     gridExtra_2.3      pander_0.6.2      
[10] latex2exp_0.4.0    knitr_1.22         Cairo_1.5-9       
[13] ggplot2_3.1.0      tibble_2.0.99.9000 readr_1.1.1       
[16] rslurm_0.4.0       Rcpp_1.0.1         reshape2_1.4.3    
[19] stringr_1.4.0      tidyr_0.8.2        purrr_0.3.2       
[22] dplyr_0.8.0.1     

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5   xfun_0.6           listenv_0.7.0     
 [4] colorspace_1.3-2   htmltools_0.3.6    viridisLite_0.3.0 
 [7] yaml_2.2.0         utf8_1.1.4         rlang_0.3.4       
[10] pillar_1.3.1.9000  glue_1.3.1.9000    withr_2.1.2       
[13] RColorBrewer_1.1-2 plyr_1.8.4         munsell_0.5.0     
[16] gtable_0.2.0       workflowr_1.3.0    rvest_0.3.2       
[19] codetools_0.2-15   evaluate_0.13      labeling_0.3      
[22] fansi_0.4.0        scales_1.0.0       backports_1.1.2   
[25] fs_1.3.1           hms_0.4.2          digest_0.6.18     
[28] stringi_1.4.3      rprojroot_1.3-2    cli_1.1.0         
[31] tools_3.5.1        magrittr_1.5       lazyeval_0.2.2    
[34] crayon_1.3.4       whisker_0.3-2      pkgconfig_2.0.2   
[37] xml2_1.2.0         assertthat_0.2.1   rmarkdown_1.10    
[40] httr_1.4.0         rstudioapi_0.10    R6_2.4.0          
[43] globals_0.12.4     git2r_0.23.0       compiler_3.5.1