Skip to contents

This vignette will plot the variable permutation importance for all models. Variable importance is critical for me to explore because it helps me to identify which environmental factors are most influential to the presence of SLF based on each model. Variable permutation importance is especially important for discerning fundamental differences between the regional-scale models. It will help us to discern whether SLF occupies a fundamentally different climate in different regions. We can also compare the permutation importance across models to discern whether SLF presence in our different regions is driven by the same or different environmental factors.

Finally, the variable permutation importance will tell us if our ensemble of regional-scale models is actually highlighting a different range of climatic factors that is traditionally missed by a global-scale model. The global-scale model is essentially the mean of all data, but this ignores invasion history (SLF has occupied different regions at different times) and the fact that SLF may actually occupy fundamentally different climates in different regions. By modeling these regions separately and ensembling them into a “regional_ensemble” model, we hypothesized this model would emphasize different climatic variables than the global mean model. We also believe that bio 11, which represents winter temperature minimums, will significantly drive our models, and therefore SLF distribution, across regions.

We will plot the variable importance for all models on the same plot for direct comparison.

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) # for running MaxEnt
library(plotROC) # plots ROCs

# html tools
library(kableExtra)
library(webshot)
library(webshot2)
mypath <- file.path(here::here() %>% 
                       dirname(),
                     "maxent/models")
ensemble_colors <- c(
  "Rn (native)" = "#4daf4a",
  "Ri.NAmerica" =  "#e41a1c",
  "Ri.Asia" = "#377eb8"
)

Variable importance

First, import the variable importance output from SDMtune.

regional_native_var_imp <- read.csv(file = file.path(mypath, "slf_regional_native_v5", "regional_native_variable_importance.csv"))

regional_invaded_n_american_var_imp <- read.csv(file = file.path(mypath, "slf_regional_invaded_n_american_v9", "regional_invaded_n_american_variable_importance.csv"))

regional_invaded_asian_var_imp <- read.csv(file = file.path(mypath, "slf_regional_invaded_asian_v4", "regional_invaded_asian_variable_importance.csv"))

Next, use a built-in function from SDMtune to get the base plot for each model.

I had to manually manipulate the ggplot objects according to the order of the bars from least to greatest. Essentially, I reordered the elements of plot_object[["data"]][[1]][["x"]] so that they were in order that would put the amounts of [["data"]][[1]][["y"]], which represent the proportion of each bar, in order from Bio 2 -> bio 15. Note that these values must be changed with every new model. So when I re-order them this way, the numeric amount in [["data"]][[1]][["y"]] should assigned the correct number value, where: * 1 = Bio 1 * 2 = Bio 11 * 3 = Bio 12 * 4 = Bio 15

regional_native_var_imp_plot <- SDMtune::plotVarImp(
  df = regional_native_var_imp
) %>%
  ggplot_build()

# change groups
regional_native_var_imp_plot[["data"]][[1]][["x"]] <- c(1, 4, 3, 2)



regional_invaded_n_american_var_imp_plot <- SDMtune::plotVarImp(
  df = regional_invaded_n_american_var_imp
) %>%
  ggplot_build()

# change groups
 regional_invaded_n_american_var_imp_plot[["data"]][[1]][["x"]] <- c(3, 1, 2, 4)



regional_invaded_asian_var_imp_plot <- SDMtune::plotVarImp(
  df = regional_invaded_asian_var_imp
) %>%
  ggplot_build()

# change groups
regional_invaded_asian_var_imp_plot[["data"]][[1]][["x"]] <- c(3, 1, 4, 2)

Finally, plot the regional ensemble models.

var_imp_ensemble <- ggplot() +
  # native model data
  geom_col(data = regional_native_var_imp_plot$data[[1]], aes(x = x + 0.2, y = y, fill = "Rn (native)"), color = "black", width = 0.2) +
  # invaded_n_american model data
  geom_col(data = regional_invaded_n_american_var_imp_plot$data[[1]], aes(x = x, y = y, fill = "Ri.NAmerica"), color = "black", width = 0.2) +
  # invaded_asian model data
  geom_col(data = regional_invaded_asian_var_imp_plot$data[[1]], aes(x = x - 0.2, y = y, fill = "Ri.Asia"), color = "black", width = 0.2) +
  labs(
    title = "Variable Importance for 'regional_ensemble' component models",
    x = "",
    y = "Permutation importance"
    ) +
  scale_x_continuous(
    breaks = c(1, 2, 3, 4),
    labels = c("bio 2", "bio 11", "bio 12", "bio 15")
  ) +
  scale_y_continuous(labels = scales::percent) +
    # aes
  theme_bw() +
  scale_fill_manual(
    name = "model",
    values = ensemble_colors,
    aesthetics = "fill"
  ) +
  theme(legend.position = "bottom") +
  coord_flip()

Now, repeat the same process, but add the global-scale model to the same plot for direct comparison.

global_var_imp <- read.csv(file = file.path(mypath, "slf_global_v5", "global_variable_importance.csv"))
global_var_imp_plot <- SDMtune::plotVarImp(
  df = global_var_imp
) %>%
  ggplot_build()

global_var_imp_plot[["data"]][[1]][["x"]] <- c(1, 4, 3, 2)
# re-define ensemble colors
ensemble_global_colors <- c(
  "Rn (native)" = "#4daf4a",
  "Ri.NAmerica" =  "#e41a1c",
  "Ri.Asia" = "#377eb8",
  "Global" = "darkgrey"
)

# plot all together
var_imp_ensemble_global <- ggplot() +
  # native model data
  geom_col(data = global_var_imp_plot$data[[1]], aes(x = x - 0.4, y = y, fill = "Global"), color = "black", width = 0.2) +
  # native model data
  geom_col(data = regional_native_var_imp_plot$data[[1]], aes(x = x - 0.2, y = y, fill = "Rn (native)"), color = "black", width = 0.2) +
  # invaded_n_american model data
  geom_col(data = regional_invaded_n_american_var_imp_plot$data[[1]], aes(x = x, y = y, fill = "Ri.NAmerica"), color = "black", width = 0.2) +
  # invaded_asian model data
  geom_col(data = regional_invaded_asian_var_imp_plot$data[[1]], aes(x = x + 0.2, y = y, fill = "Ri.Asia"), color = "black", width = 0.2) +
  labs(
    title = "Variable Importance for all models of SLF establishment potential",
    x = "",
    y = "Permutation importance"
    ) +
  scale_x_continuous(
    breaks = c(1, 2, 3, 4),
    labels = c("bio 2", "bio 11", "bio 12", "bio 15")
  ) +
  scale_y_continuous(labels = scales::percent) +
    # aes
  theme_bw() +
  scale_fill_manual(
    name = "model",
    values = ensemble_global_colors,
    aesthetics = "fill"
  ) +
  theme(legend.position = "bottom") +
  coord_flip()

Right off, we can see important differences in the variable importance between the global and regional ensemble models.

The global mean model emphasizes Bio 11 (minimum winter temperature) as the most important variable by far; we believe this is rooted in the biology of the species (Lee et al, 2011). The other regional-scale models except for the model of the invaded region in North America agree on this point. The invaded North American model places the highest importance on Bio 15 (precipitation seasonality). This is unexpected because this model closely mirrors the global mean in this variable, but it does depart significantly from the invaded_asian and native regional models (see below). The different variable emphasis may have resulted from the loss of data in those ranges. The model of the invaded region in asia also considers Bio 15 to be the most important after Bio 11, which reinforces its importance. It is likely that part of the importance of this variable is related to the monsoonal differences in precipitation seasonality, which is also reflected in the response curves below.

These findings are consistent with our hypothesis that SLF presence is driven by different environmental factors in different regions. Another of our initial hypotheses was that winter temperature minimums would drive the distribution of SLF, which we are seeing from most of our models, including the global-scale model. However, the global-scale model considers Bio 15 to be the second least important, while two of the regional-scale models clearly emphasize its importance. Indeed, our regional_ensemble is clarifying important variables we would not have seen with a global mean model.

References

  1. Lee, J.-S., Kim, I.-K., Koh, S.-H., Cho, S. J., Jang, S.-J., Pyo, S.-H., & Choi, W. I. (2011). Impact of minimum winter temperature on Lycorma delicatula (Hemiptera: Fulgoridae) egg mortality. Journal of Asia-Pacific Entomology, 14(1), 123–125. https://doi.org/10.1016/j.aspen.2010.09.004