Analyze biological realism of model response curves

Overview

In this vignette, I will plot the marginal and univariate response curves for the global and each of the regional-scale models. These plots will help me to explore the biological realism of our modeled scales.

I hypothesized that the regional models would have more nuanced responses to the environmental variables than the global model, so these plots will allow me to explore the differences in response curves between these models. These regional models should more closely is because the regional models are trained on a smaller subset of the global data, and thus may be more sensitive to the range of our covariates present in each region.

Note for unsmoothed versions of the same curves, simply change the geom_smooth functions to geom_line and remove the se and method arguments.

Setup

# general tools
library(tidyverse)  #data manipulation
library(here) #making directory pathways easier on different instances
# here::here() starts at the root folder of this package.
library(devtools)

# SDMtune and dependencies
library(SDMtune) # main package used to run SDMs
library(dismo) # package underneath SDMtune
library(rJava)
library(plotROC) # plots ROCs

# spatial data handling
library(raster) 
library(terra) 
library(sf)

library(viridis)
library(patchwork)

Note: I will be setting the global options of this document so that only certain code chunks are rendered in the final .html file. I will set the eval = FALSE so that none of the code is re-run (preventing files from being overwritten during knitting) and will simply overwrite this in chunks with plots.

These plots depict the activity of the variables put into the MaxEnt models. I will create a combined figure for the 3 models in the regional_ensemble

mypath <- file.path(here::here() %>% 
                       dirname(),
                     "maxent/models")

ensemble_colors <- c(
  "Rn (native)" = "#4daf4a",
  "Ri.NAmerica" =  "#e41a1c",
  "Ri.Asia" = "#377eb8"
)

marginal_plots_style <- list(
  theme_bw(),
  # labs
  labs(
    title = "Marginal response of 'regional_ensemble' models",
    color = "Model",
    caption = "rug plots indicate presences"
    ),
  ylab("cloglog_output"),
  # aes
  scale_color_manual(
    name = "model",
    values = ensemble_colors,
    aesthetics = "color"
  ),
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0, 1)),
  theme(legend.position = "bottom"),
  geom_hline(aes(yintercept = 0), color = "black", linewidth = 0.5) 
)

regional_native_model <- read_rds(file = file.path(mypath, "slf_regional_native_v3", "regional_native_model.rds"))

regional_invaded_model <- read_rds(file = file.path(mypath, "slf_regional_invaded_v7", "regional_invaded_model.rds"))

regional_invaded_asian_model <- read_rds(file = file.path(mypath, "slf_regional_invaded_asian_v2", "regional_invaded_asian_model.rds"))

Bio 11

regional_native_marginal_bio11 <- SDMtune::plotResponse(
  model = regional_native_model,
  var = "bio11",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()


regional_invaded_marginal_bio11 <- SDMtune::plotResponse(
  model = regional_invaded_model,
  var = "bio11",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

regional_invaded_asian_marginal_bio11 <- SDMtune::plotResponse(
  model = regional_invaded_asian_model,
  var = "bio11",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio11_max_activity_values <- c(
  native_bio_11 = regional_native_marginal_bio11$data[[1]][which.max(regional_native_marginal_bio11$data[[1]]$y), ]$x,
  invaded_bio_11 = regional_invaded_marginal_bio11$data[[1]][which.max(regional_invaded_marginal_bio11$data[[1]]$y), ]$x,
  invaded_asian_bio_11 = regional_invaded_asian_marginal_bio11$data[[1]][which.max(regional_invaded_asian_marginal_bio11$data[[1]]$y), ]$x
)

bio11_marginal_combined <- ggplot() +
  # native model data
  geom_smooth(data = regional_native_marginal_bio11$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_marginal_bio11$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # invaded_asian model data
  geom_smooth(data = regional_invaded_asian_marginal_bio11$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # Rn (native)
  geom_vline(aes(xintercept = bio11_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio11_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio11_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_marginal_bio11$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_marginal_bio11$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_marginal_bio11$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  marginal_plots_style +
  # labels
  labs(
    subtitle = "Bio 11 (mean temperature of the coldest quarter)",
    x = "mean temperature of the coldest quarter (C)"
    ) +
  # aesthetics
  xlim(-20, 20)

ggsave(
  bio11_marginal_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_bio11_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 12

regional_native_marginal_bio12 <- SDMtune::plotResponse(
  model = regional_native_model,
  var = "bio12",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()


regional_invaded_marginal_bio12 <- SDMtune::plotResponse(
  model = regional_invaded_model,
  var = "bio12",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

regional_invaded_asian_marginal_bio12 <- SDMtune::plotResponse(
  model = regional_invaded_asian_model,
  var = "bio12",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio12_max_activity_values <- c(
  native_bio_12 = regional_native_marginal_bio12$data[[1]][which.max(regional_native_marginal_bio12$data[[1]]$y), ]$x,
  invaded_bio_12 = regional_invaded_marginal_bio12$data[[1]][which.max(regional_invaded_marginal_bio12$data[[1]]$y), ]$x,
  invaded_asian_bio_12 = regional_invaded_asian_marginal_bio12$data[[1]][which.max(regional_invaded_asian_marginal_bio12$data[[1]]$y), ]$x
)

bio12_marginal_combined <- ggplot() +
  # native model data
  geom_smooth(data = regional_native_marginal_bio12$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_marginal_bio12$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_marginal_bio12$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # Rn (native)
  geom_vline(aes(xintercept = bio12_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio12_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio12_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_marginal_bio12$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_marginal_bio12$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_marginal_bio12$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  marginal_plots_style +
  # labels
  labs(
    subtitle = "Bio 12 (Annual Precipitation)",
    x = "annual precipitation (mm)"
    ) +
  # aesthetics
  xlim(0, 4250)

ggsave(
  bio12_marginal_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_bio12_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 15

regional_native_marginal_bio15 <- SDMtune::plotResponse(
  model = regional_native_model,
  var = "bio15",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()


regional_invaded_marginal_bio15 <- SDMtune::plotResponse(
  model = regional_invaded_model,
  var = "bio15",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

regional_invaded_asian_marginal_bio15 <- SDMtune::plotResponse(
  model = regional_invaded_asian_model,
  var = "bio15",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio15_max_activity_values <- c(
  native_bio_15 = regional_native_marginal_bio15$data[[1]][which.max(regional_native_marginal_bio15$data[[1]]$y), ]$x,
  invaded_bio_15 = regional_invaded_marginal_bio15$data[[1]][which.max(regional_invaded_marginal_bio15$data[[1]]$y), ]$x,
  invaded_asian_bio_15 = regional_invaded_asian_marginal_bio15$data[[1]][which.max(regional_invaded_asian_marginal_bio15$data[[1]]$y), ]$x
)

bio15_marginal_combined <- ggplot() +
  # native model data
  geom_smooth(data = regional_native_marginal_bio15$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_marginal_bio15$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # invaded_asian model data
  geom_smooth(data = regional_invaded_asian_marginal_bio15$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # Rn (native)
  geom_vline(aes(xintercept = bio15_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio15_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio15_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_marginal_bio15$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_marginal_bio15$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_marginal_bio15$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  marginal_plots_style +
  # labels
  labs(
    subtitle = "Bio 15 (Precipitation Seasonality, CV)",
    x = "precipitation seasonality (%)"
    ) +
  # aesthetics
  xlim(0, 120)

ggsave(
  bio15_marginal_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_bio15_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 2

regional_native_marginal_bio2 <- SDMtune::plotResponse(
  model = regional_native_model,
  var = "bio2",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()


regional_invaded_marginal_bio2 <- SDMtune::plotResponse(
  model = regional_invaded_model,
  var = "bio2",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

regional_invaded_asian_marginal_bio2 <- SDMtune::plotResponse(
  model = regional_invaded_asian_model,
  var = "bio2",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio2_max_activity_values <- c(
  native_bio_2 = regional_native_marginal_bio2$data[[1]][which.max(regional_native_marginal_bio2$data[[1]]$y), ]$x,
  invaded_bio_2 = regional_invaded_marginal_bio2$data[[1]][which.max(regional_invaded_marginal_bio2$data[[1]]$y), ]$x,
  invaded_asian_bio_2 = regional_invaded_asian_marginal_bio2$data[[1]][which.max(regional_invaded_asian_marginal_bio2$data[[1]]$y), ]$x
)

bio2_marginal_combined <- ggplot() +
  # native model data
  geom_smooth(data = regional_native_marginal_bio2$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_marginal_bio2$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # invaded_asian model data
  geom_smooth(data = regional_invaded_asian_marginal_bio2$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # Rn (native)
  geom_vline(aes(xintercept = bio2_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio2_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio2_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_marginal_bio2$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_marginal_bio2$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_marginal_bio2$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  marginal_plots_style +
  # labels
  labs(
    subtitle = "Bio 2 (Mean Diurnal Range)",
    x = "mean diurnal range (C)"
    ) +
  # aesthetics
  xlim(0, 15)

ggsave(
  bio2_marginal_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_bio2_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

create univariate response curves

I will smooth these plots using a GAM formula, after the methods of Elith et al, 2011.

A note that these curves were created using a Beta multiplier of 1.5.

ensemble_colors <- c(
  "Rn (native)" = "#4daf4a",
  "Ri.NAmerica" =  "#e41a1c",
  "Ri.Asia" = "#377eb8"
)

univar_plots_style <- list(
  theme_bw(),
  # labs
  labs(
    title = "Univariate response in 'regional_ensemble' models",
    color = "Model",
    caption = "rug plots indicate presences"
    ),
  ylab("cloglog_output"),
  # aes
  scale_color_manual(
    name = "model",
    values = ensemble_colors,
    aesthetics = "color"
  ),
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0, 1)),
  theme(legend.position = "bottom"),
  geom_hline(aes(yintercept = 0), color = "black", linewidth = 0.5) 
)

Bio 11

regional_native_univar_bio11 <- SDMtune::plotResponse(
  model = regional_native_model,
  var = "bio11",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()


regional_invaded_univar_bio11 <- SDMtune::plotResponse(
  model = regional_invaded_model,
  var = "bio11",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

regional_invaded_asian_univar_bio11 <- SDMtune::plotResponse(
  model = regional_invaded_asian_model,
  var = "bio11",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio11_univar_max_activity_values <- c(
  native_bio11 = regional_native_univar_bio11$data[[1]][which.max(regional_native_univar_bio11$data[[1]]$y), ]$x,
  invaded_bio11 = regional_invaded_univar_bio11$data[[1]][which.max(regional_invaded_univar_bio11$data[[1]]$y), ]$x,
  invaded_asian_bio11 = regional_invaded_asian_univar_bio11$data[[1]][which.max(regional_invaded_asian_univar_bio11$data[[1]]$y), ]$x
)

bio11_univar_combined <- ggplot() +
  # native model data
  geom_smooth(data = regional_native_univar_bio11$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_univar_bio11$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_univar_bio11$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # Rn (native)
  geom_vline(aes(xintercept = bio11_univar_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio11_univar_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio11_univar_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_univar_bio11$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_univar_bio11$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_univar_bio11$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  univar_plots_style +
  # labels
  labs(
    subtitle = "Bio 11 (mean temperature of the coldest quarter)",
    x = "mean temperature of the coldest quarter (C)"
    ) +
  # aesthetics
  xlim(-20, 20)

ggsave(
  bio11_univar_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_bio11_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 26 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_smooth()`).

## Warning: Removed 31 rows containing missing values or values outside the scale range
## (`geom_smooth()`).

## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_smooth()`).

Bio 12

regional_native_univar_bio12 <- SDMtune::plotResponse(
  model = regional_native_model,
  var = "bio12",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()


regional_invaded_univar_bio12 <- SDMtune::plotResponse(
  model = regional_invaded_model,
  var = "bio12",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

regional_invaded_asian_univar_bio12 <- SDMtune::plotResponse(
  model = regional_invaded_asian_model,
  var = "bio12",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio12_univar_max_activity_values <- c(
  native_bio12 = regional_native_univar_bio12$data[[1]][which.max(regional_native_univar_bio12$data[[1]]$y), ]$x,
  invaded_bio12 = regional_invaded_univar_bio12$data[[1]][which.max(regional_invaded_univar_bio12$data[[1]]$y), ]$x,
  invaded_asian_bio12 = regional_invaded_asian_univar_bio12$data[[1]][which.max(regional_invaded_asian_univar_bio12$data[[1]]$y), ]$x
)

bio12_univar_combined <- ggplot() +
  # native model data
  geom_smooth(data = regional_native_univar_bio12$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_univar_bio12$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # invaded_asian model data
  geom_smooth(data = regional_invaded_asian_univar_bio12$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # Rn (native)
  geom_vline(aes(xintercept = bio12_univar_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio12_univar_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio12_univar_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_univar_bio12$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_univar_bio12$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_univar_bio12$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  univar_plots_style +
  # labels
  labs(
    subtitle = "Bio 12 (Annual Precipitation)",
    x = "annual precipitation (mm)"
    ) +
  # aesthetics
  xlim(0, 4250)

ggsave(
  bio12_univar_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_bio12_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 15

regional_native_univar_bio15 <- SDMtune::plotResponse(
  model = regional_native_model,
  var = "bio15",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()


regional_invaded_univar_bio15 <- SDMtune::plotResponse(
  model = regional_invaded_model,
  var = "bio15",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

regional_invaded_asian_univar_bio15 <- SDMtune::plotResponse(
  model = regional_invaded_asian_model,
  var = "bio15",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio15_univar_max_activity_values <- c(
  native_bio15 = regional_native_univar_bio15$data[[1]][which.max(regional_native_univar_bio15$data[[1]]$y), ]$x,
  invaded_bio15 = regional_invaded_univar_bio15$data[[1]][which.max(regional_invaded_univar_bio15$data[[1]]$y), ]$x,
  invaded_asian_bio15 = regional_invaded_asian_univar_bio15$data[[1]][which.max(regional_invaded_asian_univar_bio15$data[[1]]$y), ]$x
)

bio15_univar_combined <- ggplot() +
  # native model data
  geom_smooth(data = regional_native_univar_bio15$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_univar_bio15$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # invaded_asian model data
  geom_smooth(data = regional_invaded_asian_univar_bio15$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # Rn (native)
  geom_vline(aes(xintercept = bio15_univar_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio15_univar_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio15_univar_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_univar_bio15$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_univar_bio15$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_univar_bio15$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  univar_plots_style +
  # labels
  labs(
    subtitle = "Bio 15 (Precipitation Seasonality, CV)",
    x = "precipitation seasonality (%)"
    ) +
  # aesthetics
  xlim(0, 120)

ggsave(
  bio15_univar_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_bio15_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 2

regional_native_univar_bio2 <- SDMtune::plotResponse(
  model = regional_native_model,
  var = "bio2",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()


regional_invaded_univar_bio2 <- SDMtune::plotResponse(
  model = regional_invaded_model,
  var = "bio2",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

regional_invaded_asian_univar_bio2 <- SDMtune::plotResponse(
  model = regional_invaded_asian_model,
  var = "bio2",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio2_univar_max_activity_values <- c(
  native_bio2 = regional_native_univar_bio2$data[[1]][which.max(regional_native_univar_bio2$data[[1]]$y), ]$x,
  invaded_bio2 = regional_invaded_univar_bio2$data[[1]][which.max(regional_invaded_univar_bio2$data[[1]]$y), ]$x,
  invaded_asian_bio2 = regional_invaded_asian_univar_bio2$data[[1]][which.max(regional_invaded_asian_univar_bio2$data[[1]]$y), ]$x
)

bio2_univar_combined <- ggplot() +
  # native model data
  geom_smooth(data = regional_native_univar_bio2$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_univar_bio2$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_univar_bio2$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # Rn (native)
  geom_vline(aes(xintercept = bio2_univar_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio2_univar_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio2_univar_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_univar_bio2$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_univar_bio2$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_univar_bio2$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  univar_plots_style +
  # labels
  labs(
    subtitle = "Bio 2 (Mean Diurnal Range)",
    x = "mean diurnal range (C)"
    ) +
  # aesthetics
  xlim(0, 15)

ggsave(
  bio2_univar_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_bio2_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Global Model

Load in objects and setup

univar_plots_style_global <- list(
  theme_bw(),
  # labs
  labs(
    title = "Univariate response of 'global' model",
    caption = "rug plot indicates presences"
    ),
  ylab("cloglog_output"),
  # aes
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0, 1)),
  geom_hline(aes(yintercept = 0), color = "black", linewidth = 0.5) 
)

global_model <- read_rds(file = file.path(mypath, "slf_global_v3", "global_model.rds"))

marginal_plots_style_global <- list(
  theme_bw(),
  # labs
  labs(
    title = "Marginal response of 'global' model",
    caption = "rug plot indicates presences"
    ),
  ylab("cloglog_output"),
  # aes
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0, 1)),
  geom_hline(aes(yintercept = 0), color = "black", linewidth = 0.5) 
)

global_marginal_bio11 <- SDMtune::plotResponse(
  model = global_model,
  var = "bio11",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()


global_marginal_bio12 <- SDMtune::plotResponse(
  model = global_model,
  var = "bio12",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

global_marginal_bio15 <- SDMtune::plotResponse(
  model = global_model,
  var = "bio15",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()
           

global_marginal_bio2 <- SDMtune::plotResponse(
  model = global_model,
  var = "bio2",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = TRUE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio11_marginal_max_activity_global <- global_marginal_bio11$data[[1]][which.max(global_marginal_bio11$data[[1]]$y), ]$x

bio12_marginal_max_activity_global <- global_marginal_bio12$data[[1]][which.max(global_marginal_bio12$data[[1]]$y), ]$x

bio15_marginal_max_activity_global <- global_marginal_bio15$data[[1]][which.max(global_marginal_bio15$data[[1]]$y), ]$x

bio2_marginal_max_activity_global <- global_marginal_bio2$data[[1]][which.max(global_marginal_bio2$data[[1]]$y), ]$x

plots

bio11_marginal_global <- ggplot() +
  geom_smooth(data = global_marginal_bio11$data[[1]], aes(x = x, y = y), color = "orangered2", se = FALSE, method = "gam") +
  # plot peak trendline and labels
  geom_vline(aes(xintercept = bio11_marginal_max_activity_global), color = "orangered2", linetype = "dashed", linewidth = 0.5) +
  # plot rug
  geom_rug(aes(x = global_marginal_bio11$data[[3]]$x), color = "orangered2", sides = "t") +
  # default style
  marginal_plots_style_global +
  # labels
  # labels
  labs(
    subtitle = "Bio 11 (mean temperature of the coldest quarter)",
    x = "mean temperature of the coldest quarter (C)"
    ) +
  # aesthetics
  xlim(-20, 20)

ggsave(
  bio11_marginal_global, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "global_bio11_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

bio12_marginal_global <- ggplot() +
  geom_smooth(data = global_marginal_bio12$data[[1]], aes(x = x, y = y), color = "orangered2", se = FALSE, method = "gam") +
  # plot peak trendline and labels
  geom_vline(aes(xintercept = bio12_marginal_max_activity_global), color = "orangered2", linetype = "dashed", linewidth = 0.5) +
  # plot rug
  geom_rug(aes(x = global_marginal_bio12$data[[3]]$x), color = "orangered2", sides = "t") +
  # default style
  marginal_plots_style_global +
  # labels
  # labels
  labs(
    subtitle = "Bio 12 (Annual Precipitation)",
    x = "annual precipitation (mm)"
    ) +
  # aesthetics
  xlim(0, 4250)

ggsave(
  bio12_marginal_global, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "global_bio12_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

bio15_marginal_global <- ggplot() +
  geom_smooth(data = global_marginal_bio15$data[[1]], aes(x = x, y = y), color = "orangered2", se = FALSE, method = "gam") +
  # plot peak trendline and labels
  geom_vline(aes(xintercept = bio15_marginal_max_activity_global), color = "orangered2", linetype = "dashed", linewidth = 0.5) +
  # plot rug
  geom_rug(aes(x = global_marginal_bio15$data[[3]]$x), color = "orangered2", sides = "t") +
  # default style
  marginal_plots_style_global +
  # labels
  labs(
    subtitle = "Bio 15 (Precipitation Seasonality, CV)",
    x = "precipitation seasonality (%)"
    ) +
  # aesthetics
  xlim(0, 120)

ggsave(
  bio15_marginal_global, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "global_bio15_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

bio2_marginal_global <- ggplot() +
  geom_smooth(data = global_marginal_bio2$data[[1]], aes(x = x, y = y), color = "orangered2", se = FALSE, method = "gam") +
  # plot peak trendline and labels
  geom_vline(aes(xintercept = bio2_marginal_max_activity_global), color = "orangered2", linetype = "dashed", linewidth = 0.5) +
  # plot rug
  geom_rug(aes(x = global_marginal_bio2$data[[3]]$x), color = "orangered2", sides = "t") +
  # default style
  marginal_plots_style_global +
  # labels
  labs(
    subtitle = "Bio 2 (Mean Diurnal Range)",
    x = "mean diurnal range (C)"
    ) +
  # aesthetics
  xlim(0, 15)

ggsave(
  bio2_marginal_global, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "global_bio2_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Create univariate plots

global_univar_bio11 <- SDMtune::plotResponse(
  model = global_model,
  var = "bio11",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()


global_univar_bio12 <- SDMtune::plotResponse(
  model = global_model,
  var = "bio12",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

global_univar_bio15 <- SDMtune::plotResponse(
  model = global_model,
  var = "bio15",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()
           

global_univar_bio2 <- SDMtune::plotResponse(
  model = global_model,
  var = "bio2",
  type = "cloglog",
  fun = mean,
  only_presence = TRUE,
  marginal = FALSE,
  rug = TRUE
) %>%
  ggplot_build()

# peak y trendline object 
bio11_univar_max_activity_global <- global_univar_bio11$data[[1]][which.max(global_univar_bio11$data[[1]]$y), ]$x

bio12_univar_max_activity_global <- global_univar_bio12$data[[1]][which.max(global_univar_bio12$data[[1]]$y), ]$x

bio15_univar_max_activity_global <- global_univar_bio15$data[[1]][which.max(global_univar_bio15$data[[1]]$y), ]$x

bio2_univar_max_activity_global <- global_univar_bio2$data[[1]][which.max(global_univar_bio2$data[[1]]$y), ]$x

plots

bio11_univar_global <- ggplot() +
  geom_smooth(data = global_univar_bio11$data[[1]], aes(x = x, y = y), color = "orangered2", se = FALSE, method = "gam") +
  # plot peak trendline and labels
  geom_vline(aes(xintercept = bio11_univar_max_activity_global), color = "orangered2", linetype = "dashed", linewidth = 0.5) +
  # plot rug
  geom_rug(aes(x = global_univar_bio11$data[[3]]$x), color = "orangered2", sides = "t") +
  # default style
  univar_plots_style_global +
  # labels
  # labels
  labs(
    subtitle = "Bio 11 (mean temperature of the coldest quarter)",
    x = "mean temperature of the coldest quarter (C)"
    ) +
  # aesthetics
  xlim(-20, 20)

ggsave(
  bio11_univar_global, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "global_bio11_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

bio12_univar_global <- ggplot() +
  geom_smooth(data = global_univar_bio12$data[[1]], aes(x = x, y = y), color = "orangered2", se = FALSE, method = "gam") +
  # plot peak trendline and labels
  geom_vline(aes(xintercept = bio12_univar_max_activity_global), color = "orangered2", linetype = "dashed", linewidth = 0.5) +
  # plot rug
  geom_rug(aes(x = global_univar_bio12$data[[3]]$x), color = "orangered2", sides = "t") +
  # default style
  univar_plots_style_global +
  # labels
  # labels
  labs(
    subtitle = "Bio 12 (Annual Precipitation)",
    x = "annual precipitation (mm)"
    ) +
  # aesthetics
  xlim(0, 4250)

ggsave(
  bio12_univar_global, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "global_bio12_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

bio15_univar_global <- ggplot() +
  geom_smooth(data = global_univar_bio15$data[[1]], aes(x = x, y = y), color = "orangered2", se = FALSE, method = "gam") +
  # plot peak trendline and labels
  geom_vline(aes(xintercept = bio15_univar_max_activity_global), color = "orangered2", linetype = "dashed", linewidth = 0.5) +
  # plot rug
  geom_rug(aes(x = global_univar_bio15$data[[3]]$x), color = "orangered2", sides = "t") +
  # default style
  univar_plots_style_global +
  # labels
  labs(
    subtitle = "Bio 15 (Precipitation Seasonality, CV)",
    x = "precipitation seasonality (%)"
    ) +
  # aesthetics
  xlim(0, 120)

ggsave(
  bio15_univar_global, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "global_bio15_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

bio2_univar_global <- ggplot() +
  geom_smooth(data = global_univar_bio2$data[[1]], aes(x = x, y = y), color = "orangered2", se = FALSE, method = "gam") +
  # plot peak trendline and labels
  geom_vline(aes(xintercept = bio2_univar_max_activity_global), color = "orangered2", linetype = "dashed", linewidth = 0.5) +
  # plot rug
  geom_rug(aes(x = global_univar_bio2$data[[3]]$x), color = "orangered2", sides = "t") +
  # default style
  univar_plots_style_global +
  # labels
  labs(
    subtitle = "Bio 2 (Mean Diurnal Range)",
    x = "mean diurnal range (C)"
    ) +
  # aesthetics
  xlim(0, 15)

ggsave(
  bio2_univar_global, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "global_bio2_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Combine ensemble and global models

ensemble_global_colors <- c(
  "Rn (native)" = "#4daf4a",
  "Ri.NAmerica" =  "#e41a1c",
  "Ri.Asia" = "#377eb8",
  "global" = "black"
)

marginal_plots_style <- list(
  theme_bw(),
  # labs
  labs(
    title = "Marginal response of 'regional_ensemble' and global models",
    color = "Model",
    caption = "rug plots indicate presences"
    ),
  ylab("cloglog_output"),
  # aes
  scale_color_manual(
    name = "model",
    values = ensemble_global_colors,
    aesthetics = "color"
  ),
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0, 1)),
  theme(legend.position = "bottom"),
  geom_hline(aes(yintercept = 0), color = "black", linewidth = 0.5) 
)

Bio 11

bio11_marginal_combined <- ggplot() +
  # global model
  geom_smooth(data = global_marginal_bio11$data[[1]], aes(x = x, y = y, color = "global"), se = FALSE, method = "gam") +
  # native model data
  geom_smooth(data = regional_native_marginal_bio11$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_marginal_bio11$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_marginal_bio11$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # global
  geom_vline(aes(xintercept = bio11_marginal_max_activity_global, color = "global"), linetype = "dashed", linewidth = 0.5) +
  # Rn (native)
  geom_vline(aes(xintercept = bio11_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio11_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio11_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # global is just all presences so it does not need a rug plot
  # Rn (native)
  geom_rug(aes(x = regional_native_marginal_bio11$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_marginal_bio11$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_marginal_bio11$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  marginal_plots_style +
  # labels
  labs(
    subtitle = "Bio 11 (mean temperature of the coldest quarter)",
    x = "Bio 11: mean temperature of the coldest quarter (C)"
    ) +
  # aesthetics
  xlim(-20, 20)

ggsave(
  bio11_marginal_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_global_bio11_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 12

bio12_marginal_combined <- ggplot() +
  # global model
  geom_smooth(data = global_marginal_bio12$data[[1]], aes(x = x, y = y, color = "global"), se = FALSE, method = "gam") +
  # Rn (native) model data
  geom_smooth(data = regional_native_marginal_bio12$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_marginal_bio12$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_marginal_bio12$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # global
  geom_vline(aes(xintercept = bio12_marginal_max_activity_global, color = "global"), linetype = "dashed", linewidth = 0.5) +
  # Rn (native)
  geom_vline(aes(xintercept = bio12_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio12_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio12_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_marginal_bio12$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_marginal_bio12$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_marginal_bio12$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  marginal_plots_style +
  # labels
  labs(
    subtitle = "Bio 12 (Annual Precipitation)",
    x = "Bio 12: annual precipitation (mm)"
    ) +
  # aesthetics
  xlim(0, 4250)

ggsave(
  bio12_marginal_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_global_bio12_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 15

bio15_marginal_combined <- ggplot() +
  # global model
  geom_smooth(data = global_marginal_bio15$data[[1]], aes(x = x, y = y, color = "global"), se = FALSE, method = "gam") +
  # Rn (native) model data
  geom_smooth(data = regional_native_marginal_bio15$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_marginal_bio15$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_marginal_bio15$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # global
  geom_vline(aes(xintercept = bio15_marginal_max_activity_global, color = "global"), linetype = "dashed", linewidth = 0.5) +
  # Rn (native)
  geom_vline(aes(xintercept = bio15_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio15_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio15_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_marginal_bio15$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_marginal_bio15$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_marginal_bio15$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  marginal_plots_style +
  # labels
  labs(
    subtitle = "Bio 15 (Precipitation Seasonality, CV)",
    x = "Bio 15: precipitation seasonality (%)"
    ) +
  # aesthetics
  xlim(0, 120)

ggsave(
  bio15_marginal_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_global_bio15_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 2

bio2_marginal_combined <- ggplot() +
  # global model
  geom_smooth(data = global_marginal_bio2$data[[1]], aes(x = x, y = y, color = "global"), se = FALSE, method = "gam") +
  # Rn (native) model data
  geom_smooth(data = regional_native_marginal_bio2$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_marginal_bio2$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_marginal_bio2$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # global
  geom_vline(aes(xintercept = bio2_marginal_max_activity_global, color = "global"), linetype = "dashed", linewidth = 0.5) +
  # Rn (native)
  geom_vline(aes(xintercept = bio2_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio2_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio2_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_marginal_bio2$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_marginal_bio2$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_marginal_bio2$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  marginal_plots_style +
  # labels
  labs(
    subtitle = "Bio 2 (Mean Diurnal Range)",
    x = "Bio 2: mean diurnal range (C)"
    ) +
  # aesthetics
  xlim(0, 15)

ggsave(
  bio2_marginal_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_global_bio2_marginal_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

# bio 11
write_rds(
  bio11_marginal_combined,
  file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio11_marginal_resp_curve_cloglog.rds")
)
# bio 12
write_rds(
  bio12_marginal_combined,
  file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio12_marginal_resp_curve_cloglog.rds")
)
# bio 15
write_rds(
  bio15_marginal_combined,
  file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio15_marginal_resp_curve_cloglog.rds")
)
# bio 2
write_rds(
  bio2_marginal_combined,
  file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio2_marginal_resp_curve_cloglog.rds")
)

create univariate response curves

ensemble_colors <- c(
  "Rn (native)" = "#4daf4a",
  "Ri.NAmerica" =  "#e41a1c",
  "Ri.Asia" = "#377eb8",
  "global" = "black"
)

univar_plots_style <- list(
  theme_bw(),
  # labs
  labs(
    title = "Univariate response in 'regional_ensemble' and 'global' models",
    color = "Model",
    caption = "rug plots indicate presences"
    ),
  ylab("cloglog_output"),
  # aes
  scale_color_manual(
    name = "model",
    values = ensemble_colors,
    aesthetics = "color"
  ),
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0, 1)),
  theme(legend.position = "bottom"),
  geom_hline(aes(yintercept = 0), color = "black", linewidth = 0.5) 
)

Bio 11

bio11_univar_combined <- ggplot() +
  # global
  geom_smooth(data = global_univar_bio11$data[[1]], aes(x = x, y = y, color = "global"), se = FALSE, method = "gam") +
  # Rn (native) model data
  geom_smooth(data = regional_native_univar_bio11$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_univar_bio11$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_univar_bio11$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # global
   geom_vline(aes(xintercept = bio11_univar_max_activity_global, color = "global"), linetype = "dashed", linewidth = 0.5) +
  # Rn (native)
  geom_vline(aes(xintercept = bio11_univar_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio11_univar_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio11_univar_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_univar_bio11$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_univar_bio11$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_univar_bio11$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  univar_plots_style +
  # labels
  labs(
    subtitle = "Bio 11 (mean temperature of the coldest quarter)",
    x = "Bio 11: mean temperature of the coldest quarter (C)"
    ) +
  # aesthetics
  xlim(-20, 20)

ggsave(
  bio11_univar_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_global_bio11_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 12

bio12_univar_combined <- ggplot() +
  # global
  geom_smooth(data = global_univar_bio12$data[[1]], aes(x = x, y = y, color = "global"), se = FALSE, method = "gam") +
  # Rn (native) model data
  geom_smooth(data = regional_native_univar_bio12$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_univar_bio12$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_univar_bio12$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # global
  geom_vline(aes(xintercept = bio12_univar_max_activity_global, color = "global"), linetype = "dashed", linewidth = 0.5) +
  # Rn (native)
  geom_vline(aes(xintercept = bio12_univar_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio12_univar_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio12_univar_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_univar_bio12$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_univar_bio12$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_univar_bio12$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  univar_plots_style +
  # labels
  labs(
    subtitle = "Bio 12 (Annual Precipitation)",
    x = "Bio 12: annual precipitation (mm)"
    ) +
  # aesthetics
  xlim(0, 4250)

ggsave(
  bio12_univar_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_global_bio12_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 15

bio15_univar_combined <- ggplot() +
  # global
  geom_smooth(data = global_univar_bio15$data[[1]], aes(x = x, y = y, color = "global"), se = FALSE, method = "gam") +
  # native model data
  geom_smooth(data = regional_native_univar_bio15$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_univar_bio15$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_univar_bio15$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # global
  geom_vline(aes(xintercept = bio15_univar_max_activity_global, color = "global"), linetype = "dashed", linewidth = 0.5) +
  # Rn (native)
  geom_vline(aes(xintercept = bio15_univar_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio15_univar_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio15_univar_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_univar_bio15$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_univar_bio15$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_univar_bio15$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  univar_plots_style +
  # labels
  labs(
    subtitle = "Bio 15 (Precipitation Seasonality, CV)",
    x = "Bio 15: precipitation seasonality (%)"
    ) +
  # aesthetics
  xlim(0, 120)

ggsave(
  bio15_univar_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_global_bio15_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

Bio 2

bio2_univar_combined <- ggplot() +
  # global
  geom_smooth(data = global_univar_bio2$data[[1]], aes(x = x, y = y, color = "global"), se = FALSE, method = "gam") +
  # Rn (native) model data
  geom_smooth(data = regional_native_univar_bio2$data[[1]], aes(x = x, y = y, color = "Rn (native)"), se = FALSE, method = "gam") +
  # invaded model data
  geom_smooth(data = regional_invaded_univar_bio2$data[[1]], aes(x = x, y = y, color = "Ri.NAmerica"), se = FALSE, method = "gam") +
  # Ri.Asia model data
  geom_smooth(data = regional_invaded_asian_univar_bio2$data[[1]], aes(x = x, y = y, color = "Ri.Asia"), se = FALSE, method = "gam") +
  # plot peak trendline and labels
  # global
  geom_vline(aes(xintercept = bio15_univar_max_activity_global, color = "global"), linetype = "dashed", linewidth = 0.5) +
  # Rn (native)
  geom_vline(aes(xintercept = bio2_univar_max_activity_values[[1]], color = "Rn (native)"), linetype = "dashed", linewidth = 0.5) +
  # invaded
  geom_vline(aes(xintercept = bio2_univar_max_activity_values[[2]], color = "Ri.NAmerica"), linetype = "dashed", linewidth = 0.5) +
  # Ri.Asia
  geom_vline(aes(xintercept = bio2_univar_max_activity_values[[3]], color = "Ri.Asia"), linetype = "dashed", linewidth = 0.5) +
  # plot rug plots
  # Rn (native)
  geom_rug(aes(x = regional_native_univar_bio2$data[[2]]$x, color = "Rn (native)"), sides = "t") +
  # invaded
  geom_rug(aes(x = regional_invaded_univar_bio2$data[[2]]$x, color = "Ri.NAmerica"), sides = "t", alpha = 0.8) +
  # Ri.Asia
  geom_rug(aes(x = regional_invaded_asian_univar_bio2$data[[2]]$x, color = "Ri.Asia"), sides = "t", alpha = 0.6) +
  # default style
  univar_plots_style +
  # labels
  labs(
    subtitle = "Bio 2 (Mean Diurnal Range)",
    x = "Bio 2: mean diurnal range (C)"
    ) +
  # aesthetics
  xlim(0, 15)

ggsave(
  bio2_univar_combined, 
  filename = file.path(
    here::here(), "vignette-outputs", "figures", "regional_ensemble_global_bio2_univar_resp_curve_cloglog.jpg"
    ),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

# bio 11
write_rds(
  bio11_univar_combined,
  file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio11_univar_resp_curve_cloglog.rds")
)
# bio 12
write_rds(
  bio12_univar_combined,
  file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio12_univar_resp_curve_cloglog.rds")
)
# bio 15
write_rds(
  bio15_univar_combined,
  file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio15_univar_resp_curve_cloglog.rds")
)
# bio 2
write_rds(
  bio2_univar_combined,
  file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio2_univar_resp_curve_cloglog.rds")
)

univariate response curves

# bio 11
bio11_univar_combined <- read_rds(file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio11_univar_resp_curve_cloglog.rds")
)
# bio 12
bio12_univar_combined <- read_rds(file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio12_univar_resp_curve_cloglog.rds")
)
# bio 15
bio15_univar_combined <- read_rds(file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio15_univar_resp_curve_cloglog.rds")
)
# bio 2
bio2_univar_combined <- read_rds(file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio2_univar_resp_curve_cloglog.rds")
)

bio2_univar_combined  <- bio2_univar_combined +
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold", size = 25)
    ) +
  labs(tag = "A") 



bio11_univar_combined <- bio11_univar_combined + 
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold", size = 25)
    ) +
  labs(tag = "B") 



bio12_univar_combined <- bio12_univar_combined + 
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold", size = 25)
    ) +
  labs(tag = "C") 


  
bio15_univar_combined <- bio15_univar_combined + 
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold", size = 25)
    ) +
  labs(tag = "D") 




univar_combined_patchwork <- (bio2_univar_combined + bio11_univar_combined + bio12_univar_combined + bio15_univar_combined) +
  plot_annotation(title = "Univariate response of 'regional_ensemble' and 'global' models") +
  plot_layout(axes = "collect_y", axis_titles = "collect_y", guides = "collect") & # collect y axis and legends
  theme(legend.position = "bottom")

ggsave(
  univar_combined_patchwork, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "regional_ensemble_global_combined_univar_resp_curve_cloglog.jpg"),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

univar_combined_patchwork_pt1 <- (bio2_univar_combined + bio11_univar_combined) +
  plot_annotation(title = "Univariate response of 'regional_ensemble' and 'global' models") +
  plot_layout(axes = "collect_y", axis_titles = "collect_y", guides = "collect", nrow = 2) & # collect y axis and legends
  theme(legend.position = "bottom") 

# save
ggsave(
  univar_combined_patchwork_pt1, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "regional_ensemble_global_combined_univar_resp_curve_cloglog_half1.jpg"),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

univar_combined_patchwork_pt2 <- (bio12_univar_combined + bio15_univar_combined) +
  plot_annotation(title = "Univariate response of 'regional_ensemble' and 'global' models") +
  plot_layout(axes = "collect_y", axis_titles = "collect_y", guides = "collect", nrow = 2) & # collect y axis and legends
  theme(legend.position = "bottom") 

# save
ggsave(
  univar_combined_patchwork_pt2, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "regional_ensemble_global_combined_univar_resp_curve_cloglog_half2.jpg"),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

bio11_univar_combined <- bio11_univar_combined + 
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold", size = 20)
    ) +
  labs(tag = "A") 



bio15_univar_combined  <- bio15_univar_combined +
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold", size = 20)
    ) +
  labs(tag = "B")  




univar_combined_patchwork_main <- (bio11_univar_combined / bio15_univar_combined) +
  plot_annotation(title = "Univariate response of 'regional_ensemble' and 'global' models") +
  plot_layout(axes = "collect_y", axis_titles = "collect_y", guides = "collect") & # collect y axis and legends
  theme(legend.position = "bottom")

ggsave(
  univar_combined_patchwork_main, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "regional_ensemble_global_combined_univar_resp_curve_cloglog_main.jpg"),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

bio12_univar_combined <- bio12_univar_combined + 
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold", size = 20)
    ) +
  labs(tag = "A") 



bio2_univar_combined  <- bio2_univar_combined +
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold", size = 20)
    ) +
  labs(tag = "B") 




univar_combined_patchwork_supp <- (bio12_univar_combined / bio2_univar_combined) +
  plot_annotation(title = "Univariate response of 'regional_ensemble' and 'global' models") +
  plot_layout(axes = "collect_y", axis_titles = "collect_y", guides = "collect") & # collect y axis and legends
  theme(legend.position = "bottom")

ggsave(
  univar_combined_patchwork_supp, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "regional_ensemble_global_combined_univar_resp_curve_cloglog_supp.jpg"),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

marginal

# bio 11
bio11_marginal_combined <- read_rds(file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio11_marginal_resp_curve_cloglog.rds")
)
# bio 12
bio12_marginal_combined <- read_rds(file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio12_marginal_resp_curve_cloglog.rds")
)
# bio 15
bio15_marginal_combined <- read_rds(file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio15_marginal_resp_curve_cloglog.rds")
)
# bio 2
bio2_marginal_combined <- read_rds(file = file.path(here::here(), "vignette-outputs", "figures", "figures-rds", "regional_ensemble_global_bio2_marginal_resp_curve_cloglog.rds")
)

bio11_marginal_combined <- bio11_marginal_combined + 
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold")
    ) +
  labs(tag = "A") 



bio12_marginal_combined <- bio12_marginal_combined + 
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold")
    ) +
  labs(tag = "B") 


  
bio15_marginal_combined <- bio15_marginal_combined + 
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold")
    ) +
  labs(tag = "C") 


bio2_marginal_combined  <- bio2_marginal_combined +
  theme(
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank(),
    plot.tag.location = "panel",
    plot.tag.position = c(0.05, 0.9),
    plot.tag = element_text(face = "bold")
    ) +
  labs(tag = "D") 




marginal_combined_patchwork <- (bio11_marginal_combined + bio12_marginal_combined + bio15_marginal_combined + bio2_marginal_combined) +
  plot_annotation(title = "Marginal response of 'regional_ensemble' and 'global' models") +
  plot_layout(axes = "collect_y", axis_titles = "collect_y", guides = "collect") & # collect y axis and legends
  theme(legend.position = "bottom")

ggsave(
  marginal_combined_patchwork, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "regional_ensemble_global_combined_marginal_resp_curve_cloglog.jpg"),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = "retina"
  )

References

Elith, J., Kearney, M., & Phillips, S. (2010). The art of modelling range-shifting species. Methods in Ecology and Evolution, 1(4), 330–342. https://doi.org/10.1111/j.2041-210X.2010.00036.x

Samuel M. Owens¹

2024-08-23

Overview

Setup

create marginal response plots

Bio 11

Bio 12

Bio 15

Bio 2

create univariate response curves

Bio 11

Bio 12

Bio 15

Bio 2

Global Model

Load in objects and setup

Create marginal plots

plots

Create univariate plots

plots

Combine ensemble and global models

marginal response

Bio 11

Bio 12

Bio 15

Bio 2

create univariate response curves

Bio 11

Bio 12

Bio 15

Bio 2

Patchwork together plots

univariate response curves

marginal

References

Analyze biological realism of model response curves

Samuel M. Owens1

2024-08-23

Overview

Setup

create marginal response plots

Bio 11

Bio 12

Bio 15

Bio 2

create univariate response curves

Bio 11

Bio 12

Bio 15

Bio 2

Global Model

Load in objects and setup

Create marginal plots

plots

Create univariate plots

plots

Combine ensemble and global models

marginal response

Bio 11

Bio 12

Bio 15

Bio 2

create univariate response curves

Bio 11

Bio 12

Bio 15

Bio 2

Patchwork together plots

univariate response curves

marginal

References

Samuel M. Owens¹