Youtube tutorial on implementing the bias test from Warren et al. 2021

I've just posted a Youtube tutorial on how to implement the bias test from last year's paper "The effects of climate change on Australia’s only endemic Pokémon: Measuring bias in species distribution models".  It's a really neat way to see what sort of methodological biases might exist in a given study design, and can even tell you where in geographic space to trust the predictions your models make!

https://www.youtube.com/watch?v=HtURN7HaSLg 

Code is here:

library(ENMTools)

present.files <- list.files("~/Dropbox/Ongoing Projects/Brazilian Ants and Bias/wc/",

                            pattern = ".gri", full.names = TRUE)

env.present <- stack(present.files)

future.files <- list.files("~/Dropbox/Ongoing Projects/Brazilian Ants and Bias/wc.future/",

                           pattern = ".gri", full.names = TRUE)

env.future <- stack(future.files)

adlerzi <- enmtools.species()

adlerzi$presence.points <- read.csv("~/Dropbox/Ongoing Projects/Brazilian Ants and Bias/Real Data Analyses/ant_spp/pt_data/Procryptocerus.adlerzi.csv")

adlerzi$species.name <- "Procryptocerus adlerzi"

adlerzi <- check.species(adlerzi)

adlerzi$range <- background.buffer(adlerzi$presence.points, 500000, mask = env.present[[1]])

adlerzi <- check.species(adlerzi)

interactive.plot(adlerzi)

raster.cor.plot(env.present)

# Based on that I'm going to go with bio1, bio12, and bio15

env.present <- env.present[[c("bio1", "bio12", "bio15")]]

env.future <- env.future[[c("bio1", "bio12", "bio15")]]

env.future

env.present

# Different resolutions so we need to resample

env.present <- resample(env.present, env.future)

adlerzi.gam <- enmtools.gam(adlerzi, env.present, test.prop = 0.2, rts.reps = 100)

adlerzi.future <- predict(adlerzi.gam, env.future)

plot(adlerzi.future$suitability.plot)

plot(adlerzi.gam)

predicted.change <- adlerzi.future$suitability - adlerzi.gam$suitability

plot(predicted.change)

bias.rasters <- lapply(adlerzi.gam$rts.test$rts.models,

                       function(x) predict(env.future, x$model, type = "response") - 

                         predict(env.present, x$model, type = "response"))

bias.stack <- stack(bias.rasters)

mean.pred <- calc(bias.stack, fun = mean)

plot(mean.pred)

sd.pred <- calc(bias.stack, fun = sd)

plot(sd.pred)

# Get upper and lower bounds of 95% CI

pred.95.upper <- calc(bias.stack, fun = function(x) quantile(x, probs = 0.975, na.rm= TRUE))

plot(pred.95.upper)

pred.95.lower <- calc(bias.stack, fun = function(x) quantile(x, probs = 0.025, na.rm= TRUE))

plot(pred.95.lower)

# Find which predictions fall within the CI of expected predictions under bias-only

plot(predicted.change > pred.95.lower & predicted.change < pred.95.upper)