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)