In karinorman/biodivTS: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(echo = TRUE)

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

library(ggmap)
library(ggplot2)
library(dplyr)
library(ggthemes)
library(mapproj)
library(stringr)
library(cowplot)
library(tidyverse)
library(sp)
library(rworldmap)

load(here::here("data/model_data.rda"))
load(here::here("data/rarefyID_cell_centre.rda"))

Define color palette

my_pal <- c("#d9be3b", "#77976e", "#086788", "#cc8214", "#bfac88", "#99adbf")

#library(forcats)
source(here::here("R/coords2continent.R"))

#get rarefyid centers
included_ids <- model_data %>%
  filter(taxa != "Amphibians") %>%
  select(rarefyID, taxa, climate, realm, duration) %>%
  distinct() %>% 
  mutate(rarefyID = str_replace_all(rarefyID, pattern = "_bird", ""), 
  rarefyID = str_replace_all(rarefyID, pattern = "_mamm", ""))

study_data <- rarefyID_cell_centre %>%
  inner_join(included_ids) #%>%
  #mutate(taxa = as.factor(taxa), taxa = fct_relevel(taxa, "Birds"))

study_data <- study_data %>% 
  mutate(country = coords2continent(study_data %>% select(rarefyID_x, rarefyID_y)))

study_map <- map_data("world") %>%
  ggplot(aes(long, lat, group = group)) +
  geom_polygon(fill = "white", colour = "grey") + 
  coord_map(xlim=c(-180,180), ylim=c(-57, 90), projection = "mercator") +
  geom_point(data = study_data %>% filter(taxa == "Birds"),
             aes(rarefyID_x, rarefyID_y, color = taxa), alpha = 0.4, inherit.aes = FALSE) +
  geom_point(data = study_data %>% filter(taxa == "Mammals"),
             aes(rarefyID_x, rarefyID_y, color = taxa), alpha = 0.4, inherit.aes = FALSE) +
  #plot the US amphibian point on bottom so we can see it
  # geom_point(data = study_data %>% 
  #              #filter(taxa %in% c("Amphibians")), 
  #              filter(taxa %in% c("Amphibians"), rarefyID_y > 0),
  #            aes(rarefyID_x, rarefyID_y ), color = "#d9be3b", alpha = 0.8, inherit.aes = FALSE) +
  ggthemes::theme_map() +
  scale_color_manual(values = my_pal[2:3]) +
  theme(legend.title  = element_blank()) +
  #theme(legend.title  = element_blank(), legend.text = element_text(size = 12)) +
  guides(color = guide_legend(override.aes= list(alpha = 0.8)))

ggsave(here::here("figures/study_map.jpeg"), bg = "white", dpi = 800)

duration_hist <- study_data %>% 
  #filter(taxa == "Birds") %>%
  ggplot(aes(duration)) + 
  geom_histogram(aes(fill = taxa)) + 
  theme_classic() + 
  #theme_classic(base_size = 12) + 
  facet_wrap(~taxa, scales = "free") +
  xlab("Time Series Duration") +
  theme(strip.background = element_blank(), 
        strip.text = element_blank(),
        legend.position = "none", 
        axis.title = element_text(size=8),
        ) +
  scale_fill_manual(values = my_pal[2:3])

plot_grid(study_map, duration_hist, labels = "AUTO", nrow = 2, align = 'vh', hjust = -1, axis = "l", scale = c(1,0.9), rel_heights = c(2, 1.25))

ggsave2(here::here("figures/study_map_hist.jpeg"), bg = "white", dpi = 1000, width = 180, units = "mm")

Hexbin map

coast <- rnaturalearth::ne_coastline(scale = "small", returnclass = "sf")

basemap <- ggplot(data = coast) + geom_sf(color = "grey") + 
  ggthemes::theme_map() #+
  #coord_sf(crs = 3857)  #+
  #ggthemes::theme_map() 

#xlim=c(-180,180), ylim=c(-57, 90), 
bird_hexmap <- basemap +
  geom_hex(data = study_data %>% select(rarefyID, rarefyID_x, rarefyID_y, taxa) %>% 
             distinct() %>%
             filter(taxa == "Birds"),
             aes(x = rarefyID_x, y = rarefyID_y), color = "grey",# bins = 50, color = my_pal[2]) 
                 bins = 35) +
  scale_fill_gradient2()
             scale_fill_gradient(low = "#bbe9ff", high = "#086788")


mammal_hexmap <- basemap +
  geom_hex(data = study_data %>% select(rarefyID, rarefyID_x, rarefyID_y, taxa) %>% 
             distinct() %>%
             filter(taxa == "Mammals"),
             aes(x = rarefyID_x, y = rarefyID_y), color = "grey",# bins = 50, color = my_pal[2]) 
                 bins = 35) +
  scale_fill_gradient2()

Metric plots with trend lines

load(here::here("data/metric_model_table.rda"))

plot_metric <- function(metric_name, ylabel, yval, xlabel, data, model_coef, color_list, legend){

  if(is.na(xlabel)) xlabel <- NULL

  data <- data %>%
    filter(metric == metric_name, !is.na(value), !is.infinite(value))

  model_coef <- model_coef %>%
    filter(metric == metric_name) %>%
    select(term, estimate)

  intercept <- model_coef %>% filter(term == "(Intercept)") %>% pull(estimate)
  slope <- model_coef %>% filter(term == "year_scaled") %>% pull(estimate)

  pred <- data %>%
    select(year_scaled) %>%
    distinct() %>%
    mutate(pred = intercept + year_scaled * slope)

  trend_plot <- data %>%
    left_join(pred, by = "year_scaled") %>%
    #left_join(metadata %>% select(study_id, climate)) %>%
    ggplot(aes(x = year, y = !!sym(yval))) +
    # geom_smooth(aes(color = climate, group = rarefyID), method = "glm", se = FALSE) +
    # geom_point(aes(color = climate, group = rarefyID)) +
    geom_point(aes(group = rarefyID) ,color = "grey", size = 0.25) +
    geom_smooth(aes(color = climate, group = study_id), method = "glm", se = FALSE, size = 0.7) +
    #geom_smooth(color = "black", method = "glm", se = FALSE) +
    geom_line(mapping = aes(y = pred)) +
    theme_classic() + 
    #theme_classic(base_size = 12) + 
    #theme(legend.position = "none") +
    ylab(ylabel) +
    xlab(xlabel) +
    labs(color = "Climate") +
    scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) +
    scale_colour_manual(values = color_list, drop = FALSE)

  if (legend == TRUE){
    trend_plot <- trend_plot + theme(axis.title = element_text(size=7)) 
  } else{
    trend_plot <- trend_plot + theme(axis.title = element_text(size=7), legend.position = "none")
  }

  ggsave(paste0(here::here("figures"), "/",metric_name, "_trend.png"), trend_plot)
  return(trend_plot)
}

group_colors <- c("Global" = my_pal[1], "Polar/Temperate" = my_pal[6], "Temperate/Tropical" = my_pal[3], "Tropical" = my_pal[4], "Temperate" = my_pal[2], "Overall Mean" = "#000000")

plot_map <- tibble(metric_name = c("S", "Jaccard_base", "SES_FRic", "SES_FDiv", "SES_FEve"), 
           ylabel = c("log(S)", "log(Jaccard + 1)", "Functional Richness SES", "Functional Divergence SES", "Functional Evenness SES"),
           #ylabel = c("log(S)", "Jaccard", "Functional\nRichness SES", "Functional\nDivergence SES", "Functional\nEvenness SES"),
           yval = c("log(value)", "logvalue", "value", "value", "value"),
           xlabel = c(NA, NA, NA, "year", "year"))
m_plots <- pmap(plot_map, plot_metric, data = model_data %>% mutate("log(value)" = log(value)), 
                model_coef = metric_model_table %>% filter(model_id == "base"), 
                color_list = group_colors, 
                legend = FALSE)

plot_join <- plot_grid(plotlist=m_plots, labels = "AUTO", nrow = 3, align = 'vh',
                       hjust = -1, axis = "l", scale = 0.9)

leg_plot <- model_data %>%
  bind_rows(data.frame(climate = c("Overall Mean", "Overall Mean"), year = c(2019, 2020), value = c(1, 1))) %>%
  ggplot(aes(x = year, y = value)) +
  geom_point(color = "grey", size = 0.50) +
  geom_smooth(aes(color = climate, group = study_id), method = "glm", se = FALSE, size = 0.7) +
  theme_classic() +
  theme(axis.title = element_text(size=7)) +
  #theme(legend.title = element_text(size = 12), legend.text = element_text(size = 12)) +
  scale_colour_manual(values = group_colors, drop = FALSE, name = "Climate")

# extract the legend from one of the plots
legend <- get_legend(leg_plot)

# add the legend to the row we made earlier. Give it one-third of 
# the width of one plot (via rel_widths).
#plot_join + draw_grob(legend,  2/3.3, -0.3, .3/3.3, 1, scale = 0.6)
plot_join + draw_grob(legend,  2/3.3, -0.3,.3/3.3, 1, scale = 0.01)

ggsave2(here::here("figures", "3met_long.jpeg"), bg = "white", width = 180, unit = "mm", dpi = 1000)

Get hindcast Jaccard plot for supplement

jaccard_hind_plot <- plot_metric("Jaccard_hind", ylabel = "log(Jaccard + 1)", yval = "logvalue", xlabel = "year", data = model_data %>% mutate("log(value)" = log(value)), 
            model_coef = metric_model_table %>% filter(model_id == "base"), 
            color_list = group_colors,
            legend = TRUE)

jaccard_hind_plot + draw_grob(legend,  2/3.3, -0.3,.3/3.3, 1, scale = 0.01)

ggsave(here::here("figures", "jaccard_hind_plot.jpeg"), bg = "white", width = 6, height = 4)

Effect size plots

load(here::here("data/model_output.rda"))

metric_names <- unique(model_output$metric)
raw_metrics <- c("SES_FRic", "SES_FDiv", "SES_FEve", 
                 "CWM_diet_inv", "CWM_forstrat_ground", "CWM_forstrat_watbelowsurf", "FEve")
log_metrics <- c("FRic", "S", "FDiv")
logone_metrics <- metric_names[!metric_names %in% c(raw_metrics, log_metrics)]

durations <- model_data %>%
  select(study_id, duration) %>%
  group_by(study_id) %>%
  mutate(duration = max(duration)) %>% 
  distinct()

slopes <- model_output %>% 
  filter(term == "year_scaled", model_id == "base",
         !metric %in% c("FDiv", "FRic", "FEve", "Jaccard_next", "Jaccard_hind")) %>% 
  select(metric, estimate, p.adjust, study_id) %>%
  mutate(estimate_backtrans = case_when(
           metric %in% log_metrics ~ exp(estimate),
           metric %in% logone_metrics ~ 1 + exp(estimate),
           .default = estimate
         )) %>%
  group_by(metric) %>%
  mutate(estimate_scale = scale(estimate_backtrans)[,1]) %>%
  ungroup() %>%
    mutate(ind = case_when(
    estimate > 0 & p.adjust < 0.05 ~ "blue", 
    estimate < 0 & p.adjust < 0.05 ~ "red", 
    .default = "grey"
  )) %>%
  left_join(durations) %>%
  mutate(duration_fac = as.factor(duration)) %>%
  filter(estimate_scale < 3.5, estimate_scale > -3.5)

colgraid <- colorRampPalette(c("white", "#77976e"))
graid_pal <- colgraid(26)


plot_slopes <- function(metric_name, plot_label, keep_x_axis, color){

  if (metric_name == "empty") return(NULL)

  colgraid <- colorRampPalette(c("white", color))
  graid_pal <- colgraid(26)

  plot_data <- slopes %>%
    filter(metric == metric_name) 

  # plotting parameters
  lim_buffer <- plot_data %>% filter(!is.na(estimate)) %>% pull(estimate) %>% sd()/4
  lim <- max(abs(plot_data$estimate))

  gen_est <- plot_data %>% filter(is.na(study_id))

  metric_plot <- plot_data %>% 
    filter(!is.na(study_id)) %>% 
    ggplot() +
    geom_dotplot(aes(x = estimate, fill = duration_fac),
                 stackgroups = TRUE,
                 method = "histodot"#,
                 #binwidth = 0.1, dotsize = 1.4
    ) +
    scale_fill_manual(values = graid_pal, guide = "none", drop = FALSE) +
    geom_vline(aes(xintercept = gen_est$estimate), linetype = "longdash") +
    #coord_fixed(ratio=5)+#, ylim = c(0,5)) +
    #scale_y_continuous(limits=c(0, 1), expand = c(0, 0), breaks = seq(0, 1,1/10), labels=seq(0,10)) +
    theme_minimal() +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          axis.title.y = element_blank(), 
          axis.text.y=element_blank()#, 
          #axis.title=element_text(size=14)
          ) +
    xlim(-lim - lim_buffer, 
         lim + lim_buffer) +
    xlab("temporal slope estimate") +
    geom_hline(aes(yintercept = 0)) +
    annotate('text', x = (-lim + (2*lim_buffer)), y = 0.8, label = plot_label, size = unit(4, "pt")) 

  if(isFALSE(keep_x_axis)) {
    #metric_plot <- metric_plot + theme(axis.text.x = element_blank(), axis.title.x = element_blank())
  } else{
    #metric_plot <- metric_plot + xlim(-3.5, 3.5) +
      #geom_hline(aes(yintercept = 0))
  }

  if(gen_est$p.adjust < 0.05){
    metric_plot <- metric_plot + annotate('text', x = (gen_est$estimate + (lim_buffer*2)), y = 0.8, label = "*", size = unit(16, "pt"))
  }

  return(metric_plot)
}  


muted <- khroma::color("muted")(5)

slope_plot_df <- data.frame(metric_name = c("S", "Jaccard_base", 
                     "SES_FRic", "SES_FEve", "SES_FDiv", 
                     "CWM_diet_inv", "CWM_diet_scav", "CWM_diet_vfish", "CWM_diet_vunk", "CWM_diet_planto",
                     "CWM_diet_seed", "CWM_diet_vend", "CWM_diet_nect", "CWM_diet_fruit", "CWM_diet_vect", 
                     "CWM_forstrat_ground", "CWM_forstrat_watbelowsurf", "CWM_forstrat_wataroundsurf", 
                     "CWM_forstrat_understory", "CWM_forstrat_aerial", 
                     "CWM_forstrat_canopy", "CWM_forstrat_midhigh", "CWM_bodymass_value", rep("empty", 2)),
                     plot_label = c("Species\nRichness", "Jaccard\nSimilarity", 
                     "FRic SES", "FEve SES", "FDiv SES", 
                     "Invertebrates", "Scavenge", "Fish", "Other\nVertebrates", "Other\nPlants",
                     "Seeds", "Mammals\n& Birds", "Nectar", "Fruit", "Reptiles &\nAmphibians", 
                     "Ground", "Below Water\nSurface", "Around\nWater\nSurface", 
                     "Understory", "Aerial", "Canopy", "Midhigh", 
                     "Body Mass", rep("empty", 2)), 
                     keep_x_axis = c(FALSE, FALSE, 
                     FALSE, FALSE, TRUE, 
                     FALSE, FALSE, FALSE, FALSE, TRUE,
                     FALSE, FALSE, FALSE, FALSE, TRUE, 
                     FALSE, FALSE, FALSE, 
                     FALSE, TRUE, 
                     FALSE, FALSE, TRUE, rep(FALSE, 2)),
                     color = c(rep(muted[1], 2), rep(muted[2], 3), rep(muted[3], 10), rep(muted[4], 7), rep(muted[5], 3))) %>%
  mutate(across(keep_x_axis, as.logical))

slope_plots <- pmap(slope_plot_df, plot_slopes)

slope_plot_grid <- cowplot::plot_grid(plotlist=slope_plots, nrow = 5, align = 'vh',
                       hjust = -1, axis = "l", byrow = FALSE)

cowplot::save_plot(here::here("figures/slope_plot_grid_test.png"), slope_plot_grid, ncol = 7, nrow = 5, bg = "white", base_asp = .7, dpi = 1000)

get_legends <- function(color){
  duration_legend_plot <- slopes %>%
    ggplot(aes(y = metric, x = estimate, color = duration)) +
    geom_point() +
    scale_color_gradient(low = "white", high = color) +
    theme(legend.title = element_blank(), legend.text = element_blank(), 
          legend.background = element_rect(fill='transparent'))

  duration_legend <- get_legend(duration_legend_plot)
  leg <- gridExtra::grid.arrange(duration_legend)
  ggsave(here::here(paste0("figures/slope_legends/", color, "_legend.png")), leg, width = .6, height = 1.5,
         bg = "transparent", dpi = 800)
}

map(muted, get_legends)

Let's try dot plots of each categorical slope

library(patchwork)

group_slopes <- model_output %>% 
    filter(effect == "group_slope", !metric %in% c("Jaccard_next", "FRic", "FEve", "FDiv"), 
           !term %in% c("Global", "Polar/Temperate")) %>% 
  #get the general model estimates for foraging CWM's, since they contain only bird communities
  bind_rows(model_output %>% filter(model_id == "base", term == "year_scaled", is.na(study_id),
                        metric %in% c("CWM_forstrat_ground", "CWM_forstrat_watbelowsurf", "CWM_forstrat_wataroundsurf", 
                                                                                 "CWM_forstrat_understory", "CWM_forstrat_aerial", 
                                                                                 "CWM_forstrat_canopy", "CWM_forstrat_midhigh")) %>%
              mutate(model_id = "taxa", term = "Birds")) %>%
  select(model_id, term, estimate, p.adjust, metric) %>%
  mutate(sig = ifelse(p.adjust < 0.05, "yes", "no")) %>%
  left_join(data.frame(metric = c("S", "Jaccard_base", 
                     "SES_FRic", "SES_FEve", "SES_FDiv", 
                     "CWM_diet_inv", "CWM_diet_scav", "CWM_diet_vfish", "CWM_diet_vunk", "CWM_diet_planto",
                     "CWM_diet_seed", "CWM_diet_vend", "CWM_diet_nect", "CWM_diet_fruit", "CWM_diet_vect", 
                     "CWM_forstrat_ground", "CWM_forstrat_watbelowsurf", "CWM_forstrat_wataroundsurf", 
                     "CWM_forstrat_understory", "CWM_forstrat_aerial", 
                     "CWM_forstrat_canopy", "CWM_forstrat_midhigh", "CWM_bodymass_value"),
                     plot_label = c("Species Richness", "Jaccard Similarity", 
                     "FRic SES", "FEve SES", "FDiv SES", 
                     "CWM Invertebrate Diet", "CWM Scavenge Diet", "CWM Fish Diet", "CWM Other Vertebrates Diet", "CWM Other Plants Diet",
                     "CWM Seeds Diet", "CWM Mammals & Birds Diet", "CWM Nectar Diet", "CWM Fruit Diet", "CWM Reptiles & Amphibians Diet", 
                     "CWM Ground Foraging", "CWM Below Water Surface Foraging", "CWM Around Water Surface Foraging", 
                     "CWM Understory Foraging", "CWM Aerial Foraging", "CWM Canopy Foraging", "CWM Midhigh Foraging", 
                     "CWM Body Mass"))) %>%
  mutate(estimate_backtrans = case_when(
           metric %in% log_metrics ~ exp(estimate),
           metric %in% logone_metrics ~ 1 + exp(estimate),
           .default = estimate
         )) %>%
  group_by(metric) %>%
  mutate(numeric_metric = as.factor(cur_group_id()),
         estimate_scale = scale(estimate_backtrans)[,1]) %>%
  ungroup()

sig_contrasts <- model_output %>%
  filter(effect == "contrast", p.adjust < 0.05)

get_estimate <- function(metric_name, term_name){
  group_slopes %>% 
    filter(metric == metric_name, term == term_name) %>%
    pull(estimate_scale)
}

climate_plot <- group_slopes %>%
  filter(model_id == "climate") %>%
  ggplot(aes(y = numeric_metric, x = estimate_scale, color = term)) +
  geom_point(aes(shape = sig), size = 3) +
  geom_segment(aes(x = (get_estimate("CWM_diet_vunk", "Tropical") + 0.1), y = 11, xend = (get_estimate("CWM_diet_vunk", "Temperate") - .095)), color = "black") +
  #ylim(-7,7) +
  #scale_x_discrete(expand=c(0.001, 4)) +
  theme_bw() +
  theme(axis.title.y = element_blank(), legend.title = element_blank(),
        legend.position = "bottom", 
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.border = element_blank(),
        axis.line = element_line(colour = "black")
        ) +
  scale_shape_manual(values = c(1, 19), guide = "none") +
  #ggpubr::geom_bracket(data = group_slopes %>% filter(model_id == "climate"), xmin = unlist(get_estimate("CWM_diet_vunk", "Marine")), xmax = unlist(get_estimate("CWM_diet_vunk", "Terrestrial")), y.position = 11, label = "") +
  scale_y_discrete(labels = group_slopes %>% select(numeric_metric, plot_label) %>% distinct() %>% tibble::deframe(), drop = FALSE) +
  scale_color_manual(values = group_colors) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  guides(color=guide_legend(nrow=2, byrow=TRUE)) +
  xlab("temporal slope estimate")

climate_plot

taxa_plot <- group_slopes %>%
  filter(model_id == "taxa") %>%
  ggplot(aes(y = numeric_metric, x = estimate_scale, color = term)) +
  geom_point(aes(shape = sig), size = 3) +
  #ylim(-7,7) +
  #scale_x_discrete(expand=c(0.001, 4)) +
  theme_bw() +
  theme(axis.title.y = element_blank(), legend.title = element_blank(),
        axis.line.y = element_blank(), axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.border = element_blank(),
        axis.line = element_line(colour = "black")
        ) +
  scale_shape_manual(values = c(1, 19), guide = "none") +
  scale_color_manual(values = my_pal) +
  geom_vline(xintercept = 0, linetype = "dashed")+
  guides(color=guide_legend(nrow=2, ncol = 2, byrow=TRUE)) +
  scale_y_discrete(labels = group_slopes %>% select(numeric_metric, plot_label) %>% distinct() %>% tibble::deframe(), drop = FALSE) +
  xlab("temporal slope estimate")

taxa_plot

pa_plot <- group_slopes %>%
  filter(model_id == "protected_area") %>%
  mutate(term = ifelse(term == TRUE, "Protected", "Unprotected")) %>%
  ggplot(aes(y = numeric_metric, x = estimate_scale, color = term)) +
  geom_point(aes(shape = sig), size = 3) +
  #ylim(-7,7) +
  #scale_x_discrete(expand=c(0.001, 4)) +
  theme_bw() +
  theme(axis.title.y = element_blank(), legend.title = element_blank(),
        legend.position = "bottom",
        axis.line.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        panel.border = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.line = element_line(colour = "black")
        ) +
  scale_shape_manual(values = c(1, 19), guide = "none")+
  scale_color_manual(values = my_pal) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  guides(color=guide_legend(nrow=2, ncol = 2, byrow=TRUE)) +
  scale_y_discrete(labels = group_slopes %>% select(numeric_metric, plot_label) %>% distinct() %>% tibble::deframe(), drop = FALSE) +
  xlab("temporal slope estimate")

pa_plot

realm_plot <- group_slopes %>%
  filter(model_id == "realm") %>%
  ggplot(aes(y = numeric_metric, x = estimate_scale, color = term)) +
  geom_point(aes(shape = sig), size = 3) +
   geom_segment(aes(x = (get_estimate("CWM_diet_vunk", "Marine") + 0.075), y = 11, xend = (get_estimate("CWM_diet_vunk", "Terrestrial")) - 0.075), color = "black") +
  #ylim(-7,7) +
  #scale_x_discrete(expand=c(0.001, 4)) +
  theme_bw() +
  theme(axis.title.y = element_blank(), legend.title = element_blank(),
        axis.line.y = element_blank(), axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.border = element_blank(),
        axis.line = element_line(colour = "black")
        ) +
  scale_shape_manual(values = c(1, 19), guide = "none") +
  scale_color_manual(values = my_pal) +
  geom_vline(xintercept = 0, linetype = "dashed")+
  guides(color=guide_legend(nrow=2, ncol = 2, byrow=TRUE)) +
  scale_y_discrete(labels = group_slopes %>% select(numeric_metric, plot_label) %>% distinct() %>% tibble::deframe(), drop = FALSE) +
  xlab("temporal slope estimate")

realm_plot

# contrasts_plot_grid <- cowplot::plot_grid(plotlist = list(climate_plot, taxa_plot, pa_plot, realm_plot), nrow = 1, rel_widths = c(1.7, 1, 1, 1), rel_heights = c(1, .6, .6, .6))
# cowplot::save_plot(here::here("figures/contrasts_plot_grid.png"), contrasts_plot_grid, ncol = 4, nrow = 1, bg = "white", base_asp = 1, base_height = 5)

contrasts_plot_grid <- climate_plot + taxa_plot + pa_plot + realm_plot + plot_layout(ncol = 4, axis_titles = "collect")
ggsave(here::here("figures/contrasts_plot_grid.jpg"), contrasts_plot_grid, width = 15, heigh = 6, dpi = 1000)