modTools: simple tools for simpler modelling


This website provides a series of free software tools, mostly written in R, that can be used for various tasks related to the analysis of species’ distribution and ecology (although many of them can be applied to other purposes). Most of them are now being included in R packages — please cite the appropriate package (mentioned in the blog post featuring each function) if you use these tools, which are released under a Creative Commons Attribution license. Please note that these tools are experimental (like R, they come with absolutely no warranty) and occasionally edited with corrections or improvements, so preferably check back for the latest version before you use a function.

Leave a comment here if you use these functions in publications — this way I can cite you in the blog, track how people are using the tools, and justify the time and effort dedicated to them. Thanks are due to Arpat Ozgul and François Guilhaumon, who helped me through the steep part of the R learning curve (to me it was actually more like an overhang). Comments and suggestions are very welcome, although a quick response cannot be guaranteed!

Plot shared favourability for two competing species

The sharedFav function below implements the graphical analyses of Acevedo et al. (2010, 2012) on biogeographical interactions. It takes two vectors of favourability values at different localities for, respectively, a stronger and a weaker competing species (or two equally strong competitors), and plots their favourableness or shared favourability to assess potential competitive interactions.

sharedFav <- function(strong_F, weak_F, conf = 0.95, main = "Shared favourability") {
stopifnot(length(strong_F) == length(weak_F))
opar <- par(no.readonly = T)
par(mar = c(4, 4, 2, 4.5))
F_intersection <- fuzzyOverlay(cbind(strong_F, weak_F), op = "intersection")
F_union <- fuzzyOverlay(cbind(strong_F, weak_F), op = "union")
Fovl <- sum(F_intersection, na.rm = TRUE) / sum(F_union, na.rm = TRUE)
brks <- seq(0, 1, by = 0.1)
bins <- 1:10
bin <- cut(F_intersection, breaks = brks, labels = bins)
strong_mean <- tapply(strong_F, INDEX = bin, FUN = mean)
weak_mean <- tapply(weak_F, INDEX = bin, FUN = mean)
strong_ci <- tapply(strong_F, INDEX = bin, FUN = function(x) t.test(x, conf.level = conf, na.action = na.pass)$
weak_ci <- tapply(weak_F, INDEX = bin, FUN = function(x) t.test(x, conf.level = conf)$
strong_ci[names(Filter(is.null, strong_ci))] <- NA
weak_ci[names(Filter(is.null, weak_ci))] <- NA
strong_ci_up <- unlist(lapply(strong_ci, `[`, 2))
strong_ci_dn <- unlist(lapply(strong_ci, `[`, 1))
weak_ci_up <- unlist(lapply(weak_ci, `[`, 2))
weak_ci_dn <- unlist(lapply(weak_ci, `[`, 1))
bin_size <- table(bin)
names(bin_size) <- names(bins)
props <- bin_size / length(bin)
bin <- as.integer(bin)
bar_plot <- barplot(rep(NA, length(bins)), ylim = c(0, 1), xlab = "Favourability intersection", ylab = "Mean favourability", names.arg = brks[-1], main = main)
col_bar <- "grey50"
col_ci <- "grey"
poly_left <- mean(bar_plot[2:3])
poly_right <- mean(bar_plot[8:9])
polygon(x = c(poly_left, poly_left, poly_right, poly_right), y = c(0, 1, 1, 0), col = "lightgrey", border = NA, density = 25, angle = -45)
par(new = TRUE)
barplot(props, col = col_bar, border = FALSE, xaxt = "n", yaxt = "n", add = TRUE)
axis(side = 4, col = col_bar, col.axis = col_bar, col.ticks = col_bar, col.lab = col_bar)
mtext(side = 4, line = 3, "Proportion of localities", col = col_bar)
abline(h = 0.8, col = "grey", lty = 3)
strong_x <- bar_plot - 0.1
weak_x <- bar_plot + 0.1
arrows(strong_x, strong_ci_dn, strong_x, strong_ci_up, code = 3, length = 0.03, angle = 90, col = col_ci)
arrows(weak_x, weak_ci_dn, weak_x, weak_ci_up, code = 3, length = 0.03, angle = 90, col = col_ci)
lines(x = strong_x, y = strong_mean, lwd = 2, lty = 1)
lines(x = weak_x, y = weak_mean, lwd = 2, lty = 2)
points(x = strong_x, y = strong_mean, pch = 20)
points(x = weak_x, y = weak_mean, pch = 15, cex = 0.8)

The function returns the fuzzy overlap index and a shared favourability plot similar to those of Acevedo et al. (2010, 2012). For more details, examples of use and additional references, see the help files of these functions in package fuzzySim >= 2.2, currently available on R-Forge. Special thanks to Pelayo Acevedo and Darío Chamorro for clarifications on how to build the plot!



Acevedo P., Ward A.I., Real R. & Smith G.C. (2010) Assessing biogeographical relationships of ecologically related species using favourability functions: a case study on British deer. Diversity and Distributions, 16: 515-528

Acevedo P., Jiménez-Valverde A., Melo-Ferreira J., Real R. & Alves, P.C. (2012) Parapatric species and the implications for climate change studies: a case study on hares in Europe. Global Change Biology, 18: 1509-1519

Richerson P.J. & Lum K. (1980) Patterns of plant species and diversity in California: relation to weather and topography. American Naturalist 116: 504-536


Calculate biotic threat of a stronger over a weaker species

The bioThreat function below takes two vectors of favourability values at different localities for, respectively, a stronger and a weaker species (e.g., a superior vs. an inferior competitor, or an invasive predator vs. an unadapted native prey), and calculates the level of threat that the former may potentially pose to the latter in each locality, building on theory by Richerson & Lum (1980), later developed by Acevedo et al. (2010, 2012) and further studies (e.g., Romero et al. 2014, Muñoz et al. 2015, Chamorro et al. 2019). It depends on the FavClass function (also shown below) for classifying favourability values into low, intermediate or high (Muñoz & Real 2016). For more details, see the help files of these functions in package fuzzySim >= 2.2.

favClass = 0,
fav = 0,
breaks <= 1
fclass <- rep(NA, length(fav))
b1 <- breaks[1]
b2 <- breaks[2]
for (f in 1:length(fav)) {
if ([f])) next
if (fav[f] < b1) fclass[f] = b1 & fav[f] < b2) fclass[f] = b2) fclass[f] <- 3
if (character) {
fclass[! & fclass == 1] <- "low"
fclass[! & fclass == 2] <- "intermediate"
fclass[! & fclass == 3] <- "high"
else fclass <- as.integer(fclass)

bioThreat <- function(strong_F, weak_F, character = FALSE, ...) {
stopifnot(length(strong_F) == length(weak_F))
strong_F <- favClass(strong_F, ...)
weak_F <- favClass(weak_F, ...)
threat <- rep(NA, length(strong_F))
threat[strong_F == 3 & weak_F <= 2] <- 4
threat[strong_F == 2 & weak_F == 2] <- 3
threat[strong_F <= 2 & weak_F == 3] <- 2
threat[strong_F == 3 & weak_F == 3] <- 1
threat[strong_F == 1 | weak_F == 1] <- 0
if (character) {
threat[! & threat == 4] <- "red"
threat[! & threat == 3] <- "orange"
threat[! & threat == 2] <- "yellow"
threat[! & threat == 1] <- "green"
threat[! & threat == 0] <- "white"
else threat <- as.integer(threat)

To try it out, just copy, paste and run the code of both the above functions in your R session (or install and load package fuzzySim >= 2.2, currently available on R-Forge), and run it it for two favourability vectors for, respectively, a stronger and a weaker competitor species:

threat <- bioThreat(Fav_strong, Fav_weak)
threat_colour <- bioThreat(Fav_strong, Fav_weak, character = TRUE)

You can then map the threat values to see where your stronger competitor is most and least threatening for your weaker one.


Acevedo P., Ward A.I., Real R. & Smith G.C. (2010) Assessing biogeographical relationships of ecologically related species using favourability functions: a case study on British deer. Diversity and Distributions, 16: 515-528

Acevedo P., Jiménez-Valverde A., Melo-Ferreira J., Real R. & Alves, P.C. (2012) Parapatric species and the implications for climate change studies: a case study on hares in Europe. Global Change Biology, 18: 1509-1519

Chamorro D., Muñoz A.R., Martínez-Freiría F. & Real R. (2019) Using the fuzzy logic in the distribution modelling of competitive interactions. Poster, IBS Malaga 2019 – 9th Biennial Conference of the International Biogeography Society:

Muñoz A.R. & Real R. (2006) Assessing the potential range expansion of the exotic monk parakeet in Spain. Diversity and Distributions, 12: 656-665

Muñoz A.R., Real R. & Márquez A.L. (2015) Interacciones a escala nacional entre rapaces rupicolas en base a modelos de distribucion espacial. Los casos del buitre leonado, alimoche y aguila perdicera. Informe técnico, Universidad de Málaga & Fundación EDP:

Richerson P.J. & Lum K. (1980) Patterns of plant species diversity in California: relation to weather and topography. American Naturalist, 116:504-536

Romero D., Báez J.C., Ferri-Yáñez F., Bellido J. & Real R. (2014) Modelling favourability for invasive species encroachment to identify areas of native species vulnerability. ScientificWorldJournal, 2014: 519710

Improvements to fuzzySim functions

A few improvements were recently made to several functions in the fuzzySim package. Mainly, function modelTrim is now more independent (it used to require “attach” sometimes); multTSA allows either “AIC” or “significance” as a backward stepwise selection criterion, and provides more descriptive polynomial term names (visible when save.models = TRUE); and multGLM now has the option to include a TSA (trend surface analysis, see ?multTSA) as a spatial variable calculated individually for each species. Some other minor changes are also implemented in version 1.9 of fuzzySim, which was just released on R-Forge. The next version, 2.0, will be submitted to CRAN, so all feedback is welcome to further improve the package now!


Change the background of your plotting window

Changing the background of your screen to dark instead of light colours can save energy and also spare your eyes. RStudio offers several dark-background themes, but the plotting window is still white by default. The function below allows you to define a different colour for the background and for the remaining elements of the plots:

plot_bg <- function(back = "black", front = "white") {
par(bg = back)
par(col = front, col.axis = front, col.lab = front, col.main = front, col.sub = front, fg = front)

To try it out, just copy the function above to R, press “enter”, and then run the following commands while watching the plotting window:

plot_bg(back = "darkblue", front = "turquoise")
plot_bg(back = "black", front = "white")

The default of this function is white plots on a black background so, if that’s what you want, you just need to run the function and this simpler command:


All plots produced afterwards (at least those made with base R, or with no conflicting graphical parameters) should appear on a black window — for example:



Crop or mask rasters from a zip file

Imagine you have a zipped folder with global raster maps representing different variables (e.g. the ones you can download from WorldClim, CliMond, ENVIREM, CHELSA or ecoClimate), and you need to extract just a particular region or country that you’ll be working with. The cropZippedRasters function, provided below, automates this process without the need to previously unzip the entire folder (which can require substantial disk space). The function extracts one raster at a time from the given .zip file, imports the raster to R, crops or masks it with a user-provided raster or polygon vector map, and then (by default) deletes the unzipped raster file before unzipping the next one, therefore requiring minimal disk space.


cropZippedRasters <- function(zip.file,, msk = FALSE, rast.file.ext = ".tif", aux.file.ext = NULL, delete.unzipped = TRUE)
  # zip.file: path to the zip containing the raster maps to crop or mask
  # raster or SpatialPolygons map (in your R workspace) defining the region to crop or mask
  # msk: logical, whether to actually mask with the borders (slower) or to just crop to the extent of (the default, which can be considerably faster)
  # rast.file.ext: file extension for the raster maps on disk (e.g. ".tif", ".bil", ".rst")
  # aux.file.ext: file extension for the auxiliary files when they exist (e.g. ".hdr" for .bil raster files, or ".rdc" for Idrisi .rst files)

  rast.files <- unzip(zip.file, list = TRUE) $ Name
  var.names <- unique(tools::file_path_sans_ext(rast.files))
  n.var <- length(var.names)
  cut.rasts <- vector("list", length(var.names))
  names(cut.rasts) <- var.names

  for (i in 1:n.var) {
    message("Importing map ", i, " of ", n.var)
    message("  - unzipping file...")
    unzip(zip.file, files = paste0(var.names[i], rast.file.ext))
    if (!is.null(aux.file.ext)) unzip(zip.file, files = paste0(var.names[i], aux.file.ext))
    var.rast <- raster(paste0(var.names[i], rast.file.ext))

    message("  - cropping/masking map...")
    var.cut <- crop(var.rast,, snap = "out")
    if (msk)  var.cut <- mask(var.cut,
    cut.rasts[[i]] <- var.cut

    if (delete.unzipped) {
      message("  - deleting unzipped file...")
      unlink(list.files()[grep(pattern = paste0(var.names[i], rast.file.ext), x = list.files())])
      if (!is.null(aux.file.ext)) unlink(list.files()[grep(pattern = paste0(var.names[i], aux.file.ext), x = list.files())])
    }  # end if delete

  }  # end for i

Mind that you need to have the raster package installed for this, and a vector or raster map of the region with which you want to crop or mask the raster variables. This function is closely related to the zonalFromZip function posted a while ago; I may eventually get round to combining them into a single function.

Usage example:

wclim_we <- cropZippedRasters(zip.file = "", = westeurope, rast.file.ext = ".tif", aux.file.ext = NULL)

New features in fuzzySim functions

Just two quick notes on recent improvements to functions in the fuzzySim package:

– Thanks to a bug report by José Carlos Guerrero, the multTSA function now also works when there are not multiple but only a single species to analyse;

– Thanks to a feature request and a bit of assistance by Alba Estrada, the getPreds function now also works if the data argument is a RasterStack rather than a data frame.

Check out the new versions of these functions in the latest version of fuzzySim!