Run an Extrapolation Detection (NT2) analysis and calculate MIC

In this vignette, I will run an ExDet analysis for each regional-scale model and calculate the most important covariate based on each model. The ExDet will be used for two main purposes:

This metric will allow us to quantify the level of extrapolation that is happening with our models as they are being projected to new locations and time periods. An ExDet value below 0 or above 1 indicates extrapolation. Below 0 indicates a type 1 extrapolation within a single covariate, while a value above 1 indicates a type 2 extrapolation between multiple covariates.
We will “weight” each regional-scale model using the ExDet measurement in preparation for an ensemble model. This way, our predictions are weighted by our level of confidence and the level of climate analogy.

I will specifically use the background points that were used to train each model to perform the ExDet analysis. By doing this, I am effectively testing whether the covariates in other locations are significantly different from our model. This analysis will be performed for both the historical and climate change versions of the rasters.

In addition to the ExDet analysis, I will be performing an MIC (Most Influential Covariate) analysis. This analysis reveals which covariate is causing the extrapolation in the ExDet, if present.

Both of these analyses will be performed using the dsmextra package.

Setup

library(tidyverse)

library(here) 
# here() is set at the root folder of this package

# spatial data handling
library(terra)
library(raster)
library(dsmextra) # ExDet and MIC

# aes
library(viridis)
library(ggnewscale)
library(patchwork)
library(scales)

library(scari)

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.

This is a styling object for the exdet plots.

map_style_exdet <- list(
  theme_classic(),
  theme(
    legend.position = "bottom",
    panel.background = element_rect(fill = "azure3",colour = "azure3"),
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()    
  ),
  scale_x_continuous(expand = c(0, 0)),
  scale_y_continuous(expand = c(0, 0)),
  coord_equal()
)

# color scheme taken from: https://colorbrewer2.org/?type=diverging&scheme=RdBu&n=11

First, I will outline this helper function to format all figure numbers to no more than 2 decimal places.

format_decimal <- function(x) as.numeric(sprintf("%.*f", 2, x))

This next helper function will make mapping of our exdet results more consistent, as the scale is rather complex. This function includes the above helper.

# raster should be the raster resulting from the exdet analysis
# raster.df should be the dataframe version of this raster.

map_exdet <- function(raster, raster.df) {
  
  # rename plot layer (3rd column)
  raster.df <- dplyr::rename(raster.df, ExDet = 3) 
  
  # plot
  ggplot() +
    # default plot style
    map_style_exdet +
    # plot regular raster of values first
    geom_raster(data = raster.df, aes(x = x, y = y, fill = ExDet)) +
    # fill scale
    scale_fill_gradientn(
      name = "ExDet Score",
      # the color scale used- must be 1 color per break point, with at least 1 extra color between (1, 3, 5 and 7 are the break point colors- see breaks)
      colors = c("#67001f", "#f4a582", "white", "white", "white", "#92c5de", "#053061"),
      # the scale values themselves
      values = scales::rescale(x = c(
        min(values(raster), na.rm = TRUE), 
        0, 1, 
        max(values(raster), na.rm = TRUE)
        )),
      # the breaks indicating a change in color
      breaks = c(
        format_decimal(min(values(raster), na.rm = TRUE)), 
        0, 1, 
        format_decimal(max(values(raster), na.rm = TRUE))
        ),
      # the limits of the color scale
      limits = c(
        format_decimal(min(values(raster), na.rm = TRUE)), 
        format_decimal(max(values(raster), na.rm = TRUE))
        ),
      guide = guide_colorbar(frame.colour = "black", ticks.colour = "black", barwidth = 20, order = 1)
      ) 
  
}

Finally, the style for the MIC maps.

map_style_MIC <- list(
  theme_classic(),
  theme(
    legend.position = "bottom",
    panel.background = element_rect(fill = "azure3",colour = "azure3"),
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()    
  ),
  scale_x_continuous(expand = c(0, 0)),
  scale_y_continuous(expand = c(0, 0)),
  coord_equal()
)

# color scheme taken from: https://colorbrewer2.org/?type=diverging&scheme=RdBu&n=11

Import datasets

First, I will import and tidy the global bioclim layers for 1981-2010 and 2041-2070.

# path to directory
mypath <- file.path(here::here() %>% 
                       dirname(),
                     "maxent/historical_climate_rasters/chelsa2.1_30arcsec")

# the env covariates used to train the model global
x_global_hist_env_covariates_list <- list.files(path = file.path(mypath, "v1_maxent_10km"), pattern = "\\.asc$", full.names = TRUE) %>%
    grep("bio2_1981-2010_global.asc|bio11_1981-2010_global.asc|bio12_1981-2010_global.asc|bio15_1981-2010_global.asc", ., value = TRUE)

# path to directory
  mypath <- file.path(here::here() %>% 
                       dirname(),
                     "maxent/future_climate_rasters/chelsa2.1_30arcsec")

# the env covariates for performing xy predictions for global slf and IVR points

# SSP126
x_global_126_env_covariates_list <- list.files(path = file.path(mypath, "2041-2070_ssp126_GFDL", "v1_maxent_10km"), pattern = "\\_global.asc$", full.names = TRUE) %>%
  # dont include Access to cities
  grep(pattern = "atc_2015", value = TRUE, invert = TRUE)

# SSP370
x_global_370_env_covariates_list <- list.files(path = file.path(mypath, "2041-2070_ssp370_GFDL", "v1_maxent_10km"), pattern = "\\_global.asc$", full.names = TRUE) %>%
  # dont include Access to cities
  grep(pattern = "atc_2015", value = TRUE, invert = TRUE)

# SSP585
x_global_585_env_covariates_list <- list.files(path = file.path(mypath, "2041-2070_ssp585_GFDL", "v1_maxent_10km"), pattern = "\\_global.asc$", full.names = TRUE) %>%
  # dont include Access to cities
  grep(pattern = "atc_2015", value = TRUE, invert = TRUE)

I will unify the variable naming and import the rasters.

# layer name object. Check order of layers first
env_layer_names <- c("bio11", "bio12", "bio15", "bio2")

# stack env covariates
x_global_hist_env_covariates <- terra::rast(x = x_global_hist_env_covariates_list)

# attributes
nlyr(x_global_hist_env_covariates)
names(x_global_hist_env_covariates)
minmax(x_global_hist_env_covariates)
# ext(x_global_hist_env_covariates)
# crs(x_global_hist_env_covariates)


# I will change the name of the variables because they are throwing errors in SDMtune
names(x_global_hist_env_covariates) <- env_layer_names
# confirmed- SDMtune doesnt like dashes in column names (it is read as a mathematical operation)

# SSP126
x_global_126_env_covariates <- terra::rast(x = x_global_126_env_covariates_list)

# attributes
nlyr(x_global_126_env_covariates)
names(x_global_126_env_covariates)
minmax(x_global_126_env_covariates)
# ext(x_global_126_env_covariates)
# crs(x_global_126_env_covariates)

names(x_global_126_env_covariates) <- env_layer_names


# SSP370
x_global_370_env_covariates <- terra::rast(x = x_global_370_env_covariates_list)

# attributes
nlyr(x_global_370_env_covariates)
names(x_global_370_env_covariates)
minmax(x_global_370_env_covariates)
# ext(x_global_370_env_covariates)
# crs(x_global_370_env_covariates)

names(x_global_370_env_covariates) <- env_layer_names


# SSP585
x_global_585_env_covariates <- terra::rast(x = x_global_585_env_covariates_list)

# attributes
nlyr(x_global_585_env_covariates)
names(x_global_585_env_covariates)
minmax(x_global_585_env_covariates)
# ext(x_global_585_env_covariates)
# crs(x_global_585_env_covariates)

names(x_global_585_env_covariates) <- env_layer_names

rm(x_global_hist_env_covariates_list)
rm(x_global_126_env_covariates_list)
rm(x_global_370_env_covariates_list)
rm(x_global_585_env_covariates_list)

I also need the background point sets that were used to train each regional model. These copies need to include the values of each raster at those particular locations.

# regional_native
regional_native_background <- read_csv(
  file.path(here::here(), "vignette-outputs", "data-tables", "regional_native_background_points_with_data_v4.csv")
  ) %>%
  as.data.frame()

## Rows: 10000 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Species
## dbl (6): X, Y, bio11, bio12, bio15, bio2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# regional_invaded
regional_invaded_background <- read_csv(
  file.path(here::here(), "vignette-outputs", "data-tables", "regional_invaded_background_points_with_data_v8.csv")
  ) %>%
  as.data.frame()

## Rows: 10000 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Species
## dbl (6): X, Y, bio11, bio12, bio15, bio2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# regional_invaded_asian
regional_invaded_asian_background <- read_csv(
  file.path(here::here(), "vignette-outputs", "data-tables", "regional_invaded_asian_background_points_with_data_v3.csv")
  ) %>%
  as.data.frame()

## Rows: 10000 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Species
## dbl (6): X, Y, bio11, bio12, bio15, bio2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# I will also perform this operation for the global model

# global
global_background <- read_csv(
  file.path(here::here(), "vignette-outputs", "data-tables", "global_background_points_with_data_v4.csv")
  ) %>%
  as.data.frame()

## Rows: 20000 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Species
## dbl (6): X, Y, bio11, bio12, bio15, bio2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

I need to import the rasters of the area from which the background points were chosen so these can be displayed on each ExDet and MIC plot

# path to directory
mypath <- file.path(here::here() %>% 
                       dirname(),
                     "maxent/historical_climate_rasters/chelsa2.1_30arcsec")


# native
regional_native_background_area <- terra::rast(x = file.path(mypath, "v1_maxent_10km", "bio11_1981-2010_regional_native_KG.asc"))
# invaded
regional_invaded_background_area <- terra::rast(x = file.path(mypath, "v1_maxent_10km", "bio11_1981-2010_regional_invaded_KG.asc"))
# invaded_asian
regional_invaded_asian_background_area <- terra::rast(x = file.path(mypath, "v1_maxent_10km", "bio11_1981-2010_regional_invaded_asian_KG.asc"))


# convert to df
regional_native_background_area_df <- terra::as.data.frame(regional_native_background_area, xy = TRUE)
regional_invaded_background_area_df <- terra::as.data.frame(regional_invaded_background_area, xy = TRUE)
regional_invaded_asian_background_area_df <- terra::as.data.frame(regional_invaded_asian_background_area, xy = TRUE)

1. Extract variables for climate data

First, I will load the background points dataset used to train the model and that will be used to train the ExDet analysis. I have extracted the climate dataset from 1981-2010 for these points, but not for our projected time period of 2041-2070 and our 3 ssp scenarios. I will extract the values for the background points from each model and SSP scenario.

# regional_native model
# extract values from 2060 projected raster
regional_native_background_126 <- regional_native_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_126_env_covariates, y = ., method = "simple", ID = FALSE)




# regional_invaded model
regional_invaded_background_126 <- regional_invaded_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_126_env_covariates, y = ., method = "simple", ID = FALSE)





# regional_invaded_asian model
regional_invaded_asian_background_126 <- regional_invaded_asian_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_126_env_covariates, y = ., method = "simple", ID = FALSE)



# global
global_background_126 <- global_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_126_env_covariates, y = ., method = "simple", ID = FALSE)

# regional_native model
# extract values from 2060 projected raster
regional_native_background_370 <- regional_native_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_370_env_covariates, y = ., method = "simple", ID = FALSE)




# regional_invaded model
regional_invaded_background_370 <- regional_invaded_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_370_env_covariates, y = ., method = "simple", ID = FALSE)





# regional_invaded_asian model
regional_invaded_asian_background_370 <- regional_invaded_asian_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_370_env_covariates, y = ., method = "simple", ID = FALSE)



# global
global_background_370 <- global_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_370_env_covariates, y = ., method = "simple", ID = FALSE)

# regional_native model
# extract values from 2060 projected raster
regional_native_background_585 <- regional_native_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_585_env_covariates, y = ., method = "simple", ID = FALSE)




# regional_invaded model
regional_invaded_background_585 <- regional_invaded_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_585_env_covariates, y = ., method = "simple", ID = FALSE)





# regional_invaded_asian model
regional_invaded_asian_background_585 <- regional_invaded_asian_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_585_env_covariates, y = ., method = "simple", ID = FALSE)



# global
global_background_585 <- global_background %>%
  dplyr::select("X", "Y") %>%
  # extract new values
  terra::extract(x_global_585_env_covariates, y = ., method = "simple", ID = FALSE)

The dsmextra::compute_extrapolation() function only takes data columns (no descriptors), so I will isolate those from each background point set that was loaded in.

# native
regional_native_background <- dplyr::select(regional_native_background, 4:7)
# invaded
regional_invaded_background <- dplyr::select(regional_invaded_background, 4:7)
# regional_invaded_asian
regional_invaded_asian_background <- dplyr::select(regional_invaded_asian_background, 4:7)
# global
global_background <- dplyr::select(global_background, 4:7)

Now I will convert the global-scale climate rasters to df. These will be used for projection, to compare the exdet and MIC analyses with. The background point sets are used to train the analyses, and the global rasters are used for projection.

# historical
global_projection_hist_df <- terra::as.data.frame(x_global_hist_env_covariates, xy = TRUE)

# CMIP6
global_projection_126_df <- terra::as.data.frame(x_global_126_env_covariates, xy = TRUE)
global_projection_370_df <- terra::as.data.frame(x_global_370_env_covariates, xy = TRUE)
global_projection_585_df <- terra::as.data.frame(x_global_585_env_covariates, xy = TRUE)

2. Run ExDet / MIC Function

First, I will run the function that calculates the ExDet and MIC indices. I will use the background points datasets to train each analysis and will project to the entire global raster. This will be done for each of the historical and climate change time periods. This function outputs rasters and data frames that can be used for plotting.

To save length in this vignette, I will only display the sections creating the ExDet and MIC plots based on the historical data, as these sections are essentially copy-pasted.

# output directory
mypath <- file.path(here::here() %>% 
                       dirname(),
                     "maxent")

2.1 Historical Climate

exdet_MIC_regional_native_hist <- dsmextra::compute_extrapolation(
  samples = regional_native_background, # training set of background points
  covariate.names = env_layer_names, # names of covariates
  prediction.grid = global_projection_hist_df, # projection area for extrapolation test
  coordinate.system = "ESRI:54017", # CRS
  verbose = TRUE
)

exdet_MIC_regional_invaded_hist <- dsmextra::compute_extrapolation(
  samples = regional_invaded_background, # training set of background points
  covariate.names = env_layer_names, # names of covariates
  prediction.grid = global_projection_hist_df, # projection area for extrapolation test
  coordinate.system = "ESRI:54017", # CRS
  verbose = TRUE
)

exdet_MIC_regional_invaded_asian_hist <- dsmextra::compute_extrapolation(
  samples = regional_invaded_asian_background, # training set of background points
  covariate.names = env_layer_names, # names of covariates
  prediction.grid = global_projection_hist_df, # projection area for extrapolation test
  coordinate.system = "ESRI:54017", # CRS
  verbose = TRUE
)

also do global model

exdet_MIC_global_hist <- dsmextra::compute_extrapolation(
  samples = global_background, # training set of background points
  covariate.names = env_layer_names, # names of covariates
  prediction.grid = global_projection_hist_df, # projection area for extrapolation test
  coordinate.system = "ESRI:54017", # CRS
  verbose = TRUE
)

3. regional_native- Save/plot ExDet/MIC rasters

Now I need to pull the rasters from the above functions and save them to my working directory folder. Then, I will create plots of each raster and save those as well.

3.1 regional_native hist

First, I will save the summary files from the analysis to each model folder.

Save outputs

# summary of data
write_csv(
  x = exdet_MIC_regional_native_hist[["data"]][["all"]],
  file = file.path(mypath, "models", "slf_regional_native_v4", "exdet_mic_analysis_results.csv")
  )

# exdet
write_csv(
  x = as.data.frame(exdet_MIC_regional_native_hist[["summary"]][["extrapolation"]]),
  file = file.path(mypath, "models", "slf_regional_native_v4", "exdet_analysis_summary.csv")
  )

# mic univariate
write_csv(
  x = as.data.frame(exdet_MIC_regional_native_hist[["summary"]][["mic"]][[1]]),
  file = file.path(mypath, "models", "slf_regional_native_v4", "mic_analysis_univariate_summary.csv")
  )

# mic univariate
write_csv(
  x = as.data.frame(exdet_MIC_regional_native_hist[["summary"]][["mic"]][[2]]),
  file = file.path(mypath, "models", "slf_regional_native_v4", "mic_analysis_combinatorial_summary.csv")
  )

Next, I will save the output rasters.

# exdet
raster::writeRaster(
  x = exdet_MIC_regional_native_hist[["rasters"]][["ExDet"]][["all"]], 
  filename = file.path(mypath, "models", "slf_regional_native_v4", "ExDet_slf_regional_native_v4_background_projGlobe_1981-2010.tif"),
  filetype = "GTiff",
  overwrite = FALSE
  )

raster::writeRaster(
  x = exdet_MIC_regional_native_hist[["rasters"]][["mic"]][["all"]], 
  filename = file.path(mypath, "models", "slf_regional_native_v4", "MIC_slf_regional_native_v4_background_projGlobe_1981-2010.tif"),
  filetype = "GTiff",
  overwrite = FALSE
  )

Plot ExDet

Finally, I will create plots of each of the ExDet and MIC analyses.

# re-import raster
exdet_regional_native_hist <- terra::rast(x = file.path(mypath, "models", "slf_regional_native_v4", "ExDet_slf_regional_native_v4_background_projGlobe_1981-2010.tif"))
# convert to df
exdet_regional_native_hist_df <- terra::as.data.frame(exdet_regional_native_hist, xy = TRUE)
# mapping helper
exdet_regional_native_hist_plot <- map_exdet(
  raster = exdet_regional_native_hist, 
  raster.df = exdet_regional_native_hist_df
  ) 



# add new scale for bg area
exdet_regional_native_hist_plot <- exdet_regional_native_hist_plot +
  # new scale
  ggnewscale::new_scale_fill() +
  # add layer for bg point area
  geom_raster(data = regional_native_background_area_df, aes(x = x, y = y, fill = "Background")) +
  # scale 2
  scale_fill_manual(name = "", values = c("Background" = "azure4"), guide = guide_legend(order = 2)) +
  theme(legend.key = element_rect(color = "black")) +
  # labs
  labs(
    title = "ExDet analysis | bio2, bio11, bio12 and bio15 | 1981-2010",
    subtitle = "Global climate divergence from 'Rn' model training set",
    fill = "ExDet Analysis"
    )

ggsave(
  exdet_regional_native_hist_plot, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "ExDet_slf_regional_native_v4_background_projGlobe_1981-2010.jpg"),
  height = 8, 
  width = 12,
  device = jpeg,
  dpi = "retina"
  )

Plot MIC

# re-import
mic_regional_native_hist_df <- terra::rast(x = file.path(mypath, "models", "slf_regional_native_v4", "MIC_slf_regional_native_v4_background_projGlobe_1981-2010.tif")) %>%
  terra::as.data.frame(., xy = TRUE) %>%
  dplyr::rename(mic = 3) # rename mic column

# plot mean raster first
mic_regional_native_hist_plot <- ggplot() +
  # default plot style
  map_style_MIC +
  # plot regular raster of values first
  geom_raster(data = mic_regional_native_hist_df, aes(x = x, y = y, fill = as.factor(mic))) +
  # fill scale 1
  viridis::scale_fill_viridis(name = "MIC", discrete = TRUE, labels = c("analogous", "bio11", "bio12", "bio15", "bio2"), guide = guide_legend(order = 1)) +
  # new scale
  ggnewscale::new_scale_fill() +
  # add layer for bg point area
  geom_raster(data = regional_native_background_area_df, aes(x = x, y = y, fill = "Background")) +
  # scale 2
  scale_fill_manual(name = "", values = c("Background" = "azure4"), guide = guide_legend(order = 2)) +
  theme(legend.key = element_rect(color = "black")) +
  # aes
  labs(
    title = "Most Influential Covariate (MIC) analysis | 1981-2010",
    subtitle = "Covariate causing most extrapolation from 'Rn (native)' model training set"
    )

ggsave(
  mic_regional_native_hist_plot, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "MIC_slf_regional_native_v4_background_projGlobe_1981-2010.jpg"),
  height = 8, 
  width = 12,
  device = jpeg,
  dpi = "retina"
  )

I will repeat this workflow for each of the models and their CMIP6 version climate data, so I will not display the workflow for the other models and ssp scenarios.

In the next vignette, I will utilize the ExDet analysis in my ensemble techniques for the regional-scale models. This ensembling technique should give me a unified prediction with which to compare the global-scale model and make predictions.

References

Bouchet, P. J., Miller, D. L., Roberts, J. J., Mannocci, L., Harris, C. M., & Thomas, L. (2020). dsmextra: Extrapolation assessment tools for density surface models. Methods in Ecology and Evolution, 11(11), 1464–1469. https://doi.org/10.1111/2041-210X.13469
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
Mesgaran, M. B., Cousens, R. D., & Webber, B. L. (2014). Here be dragons: A tool for quantifying novelty due to covariate range and correlation change when projecting species distribution models. Diversity and Distributions, 20(10), 1147–1159. https://doi.org/10.1111/ddi.12209

Samuel M. Owens¹

2025-07-28

Setup

Import datasets

1. Extract variables for climate data

2. Run ExDet / MIC Function

2.1 Historical Climate

3. regional_native- Save/plot ExDet/MIC rasters

3.1 regional_native hist

Save outputs

Plot ExDet

Plot MIC

References

Run an Extrapolation Detection (NT2) analysis and calculate MIC

Samuel M. Owens1

2025-07-28

Setup

Import datasets

1. Extract variables for climate data

2. Run ExDet / MIC Function

2.1 Historical Climate

3. regional_native- Save/plot ExDet/MIC rasters

3.1 regional_native hist

Save outputs

Plot ExDet

Plot MIC

References

Samuel M. Owens¹