index.Rmd

---
title: "Flexible species distribution modelling methods perform well on spatially separated testing data. Global Ecology and Biogeography"
description: |
  Appendix S1
author:
  - name: Roozbeh Valavi 
    url: https://github.com/rvalavi
    affiliation: The University of Melbourne, Australia
    affiliation_url: https://ecosystemforest.unimelb.edu.au/
    orcid_id: 0000-0003-2495-5277
  - name: Jane Elith 
    affiliation: The University of Melbourne, Australia
    affiliation_url: https://ecosystemforest.unimelb.edu.au/
    orcid_id: 0000-0002-8706-0326
  - name: José J. Lahoz-Monfort 
    affiliation: Pyrenean Institute of Ecology, Spanish National Research Council (CSIC), Spain
    affiliation_url: http://www.ipe.csic.es/conservacion-bio/
    orcid_id: 0000-0002-0845-7035
  - name: Gurutzeta Guillera-Arroita 
    affiliation: Pyrenean Institute of Ecology, Spanish National Research Council (CSIC), Spain
    affiliation_url: http://www.ipe.csic.es/conservacion-bio/
    orcid_id: 0000-0002-8387-5739
twitter: 
  site: "@ValaviRoozbeh"
  creator: "@ValaviRoozbeh"
journal: 
  title: "Global Ecology and Biogeography"
  # pasge: "1-27"
  # issn: 2490-1752
  # publisher: ESA
# volume: 44
# issue: 4
doi: "10.1111/GEB.13639"
date: "2023-01-06"
# bibliography: references.bib
# bib-humanities: true
output: 
  distill::distill_article:
    toc: true
    toc_depth: '2'
    toc_float: true
    theme: theme.css
creative_commons: CC BY
header-includes: 
  - \usepackage{caption}
  - \captionsetup[figure]{labelformat=empty}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = FALSE,
	message = FALSE,
	warning = FALSE,
	tidy = TRUE,
	# tidy.opts = list(width.cutoff = 80),
	fig.align = "center"
)
```

Appendix S1 for *Valavi, R., Elith, J., Lahoz-Monfort, J. J., & Guillera-Arroita, G. (2023). Flexible species distribution modelling methods perform well on spatially separated testing data. Global Ecology and Biogeography, 32, 369–383.* [DOI: 10.1111/geb.13639](https://doi.org/10.1111/geb.13639)



```{r}
# for automatic numbering; 
# the section number if not working for distill pages/html files
sn <- 0
```


## `r sn <- floor(sn) + 1; sn`- Covariates and TGB samples

To model the presence-only data with background samples and deal with biases in these data, we used Target-Group-Background (TGB) samples  introduced by Phillips *et al.* (2009). TGB provide a background sample with similar biases to that of the presence records. For each species, a background sample is generated by collating all presence records from the same biological group and region (including the target species), with minor adjustments to avoid multiple records per grid cell. For details see Phillips *et al.* (2009) and for code see Elith *et al.* (2020). 
    
Table 1 to 6 below show the environmental covariates used in each region. These are the variables that do not have a pairwise correlation of more than 0.8. You can find a complete list of variables with more details in Elith *et al.* (2020), and in the help of the [`disdat`](https://CRAN.R-project.org/package=disdat) R package.

```{r}
library(kableExtra)
library(magrittr)

knitr::kable(read.csv("results/AWT_environment.csv"), caption = "Environmental vartiables for AWT region (~80m spatial resolution).")

```

```{r}
knitr::kable(read.csv("results/CAN_environment.csv"), caption = "Environmental vartiables for CAN region (~1000m spatial resolution).")

```

```{r}
knitr::kable(read.csv("results/NSW_environment.csv"), caption = "Environmental vartiables for NSW region (~100m spatial resolution).")

```

```{r}
knitr::kable(read.csv("results/NZ_environment.csv"), caption = "Environmental vartiables for NZ region (~100m spatial resolution).")

```

```{r}
knitr::kable(read.csv("results/SA_environment.csv"), caption = "Environmental vartiables for SA region (~1000m spatial resolution).")

```

```{r}
knitr::kable(read.csv("results/SWI_environment.csv"), caption = "Environmental vartiables for SWI region (~100m spatial resolution).")

```



## `r sn <- sn + 1; sn`- Code and data availability

You can find species data in `disdat` R package, with the rasters available on OSF (Elith *et al.,* 2020). To run the models, you can use the codes and data provide in [OSF repository](https://osf.io/g6dc3/?view_only=2f91390cb8f14ae3ba75c35427e017e4). 
    
To reproduce our results, you can use any of the following (details below): 

* Using [`renv`](https://rstudio.github.io/renv/articles/renv.html) package recovery in local R environment
    + Use the `revn.lock` file in [OSF repository](https://osf.io/g6dc3/?view_only=2f91390cb8f14ae3ba75c35427e017e4) and retrieve the R package versions in your local system with the `renv` package

* Using [Docker](https://www.docker.com/) containers
    + Use the pre-built Docker image (`rvalavi_image.tar` stored in OSF) and create a virtual environment with the same R packages (recommended)
    + Build a Docker image based on the `Dockerfile` provided in the [OSF repository](https://osf.io/g6dc3/?view_only=2f91390cb8f14ae3ba75c35427e017e4).


### `r sn <- sn + 0.1; sn`- Using `renv` package recovery

First, create a new RStudio project and place the `renv.lock` file in it. Use the commands below to restore the R package for modelling.

```{r eval=FALSE, echo=TRUE}
install.packages("renve")

renv::init()

renv::restore()

```

You need to install a few other packages that are not listed in `renv` file.

```{r eval=FALSE, echo=TRUE}
install.packages('remotes')
install.packages('rJava')

remotes::install_github('meeliskull/prg/R_package/prg')
remotes::install_github('b0rxa/scmamp')
remotes::install_github('rvalavi/myspatial')

remotes::install_version('gam', version = '1.20', repos = 'http://cran.us.r-project.org')
remotes::install_version('gbm', version = '2.1.5', repos = 'http://cran.us.r-project.org')

```


### `r sn <- sn + 0.1; sn`- Using Docker containers

Here we explain how to use Dockers for creating a virtual system to reproduce our results. Docker can be installed on different platforms. You can find the instructions for installing Docker on different operating systems on their [website](https://docs.docker.com/get-docker/). Docker provides an RStudio installed and all the R and system packages required for running our analysis.
    
To use Docker, you can either A) load the pre-built image (recommended), or B) build a new image from `Dockerfile`. 

`A)` To load the Docker image, first download the files from OSF, then use:

```{bash eval=FALSE, echo=TRUE}
docker load --input rvalavi_image.tar

```

`B)` If you want to build the image in your local system, you need to download all the files in OSF (except "docker" folder). Then run the following terminal commands. In Linux systems you might need to use `sudo` before `docker` commands. 
    
Go to the directory of downloaded files from the OSF within terminal and run:
    
```{bash eval=FALSE, echo=TRUE}
docker build -t rvalavi:4.0 .

```

Wait until the build is complete. Then check to see the images is created.

```{bash eval=FALSE, echo=TRUE}
docker images

```

You should see `rvalavi` with TAG 4.0 listed as a Docker image.
    
After the image is loaded (A) or created (B), you need to run a container to get access to RStudio and the R packages. The Docker container is a live instance of the image. Use the following command to run a container to access RStudio.

```{bash eval=FALSE, echo=TRUE}
docker run --name rstudio -p 8787:8787 -e PASSWORD=123 -d rvalavi:4.0

```

This code has several components:  
`--name`: name of the container  
`-p`: port on which container is running. We use this to connect to rstudio  
`-e PASSWORD`: password for the rstudio server (you can choose any password)  
`-v`: mapping a directory in the local system to a directory in the container (this was not used in the code above). This will allow you to save the generated files and code in local drive and also access to code/data inside your system.  
`-d`: run the container in the background   
`rvalavi:4.0`: name and tag of the Docker image  


Now, run RStudio server from container. Open an Internet browser and go to `localhost:8787` to open RStudio. Use "**rstudio**" as username and the password you specified in the previous step (here password is "123") to open RStudio.

![](rstudio.png)

You can run the models by running the `R/nceas_modelling.R` script. The model predictions will be stored in `output/nceas_model_output` folder. To calculate evaluations, run `R/nceas_evaluation.R` script.



## `r sn <- floor(sn) + 1; sn`- Modelling methods parameters

A summary of model implementation settings are presented in Table 7, below. The “parameters” column shows model arguments in R programming that are selected in the modelling process. The “values” column shows the value or ranges of values selected for model fitting and tuning for each R function. All models are fitted in R `v4.0.0`.
     
Some methods accept weights. GAM, GLMs, and BRT use case weights (i.e., there is a weight for each training sample), and SVM utilizes class weights (i.e., there is a weight for each class; here presence and background samples are the two classes so there are only two weights). The weights are generated by giving a weight of 1 to every presence point and giving the weights to the background in a way that the sum of the weights for the presence and background are equal. For class weights in SVM, an inverse proportional weight was used.


```{r}
# knitr::kable(read.csv("model_parameters.csv"), caption = "Parameters used for implementing different modes.")
knitr::kable(read.csv("model_parameters.csv"), caption = "Parameters used for implementing different modes.")


```
\* GLM-step and GLM-lasso were fitted allowing linear and quadratic terms only, with no interactions.
     
Parameters of MaxEnt variants are presented in the following table.


```{r}
knitr::kable(read.csv("maxent_parameters.csv"), caption = "Parameters of MaxEnt model in different MaxEnt variants.")

```


## `r sn <- sn + 1; sn`- Evaluation metrics

We used $AUC_{ROC}$, $AUC_{PRG}$ and COR for evaluating the models. $AUC_{ROC}$ measures how well a model discriminates between presence and absence records in the test dataset. It can range from 0 to 1, with 1 indicating a model has perfect discrimination abilities and 0.5 showing discrimination is equivalent to that from random predictions (Pearce & Ferrier, 2000; Elith *et al.*, 2006). $AUC_{PRG}$ is similar to $AUC_{ROC}$, but less commonly used in ecology. It puts more focus on correctly predicted presences (Flach & Kull, 2015). An $AUC_{PRG}$ value of 1 shows perfect discrimination, 0 indicates random discrimination and negative denotes worse than random. Since there is no lower limit for negative values in $AUC_{PRG}$, we only estimated ranks, not mean performance, for this metric. COR is the correlation between model predictions and the presence-absence testing data (Elith *et al.*, 2006). 



## `r sn <- sn + 1; sn`- Dispersion of $AUC_{ROC}$

To further highlight the difference between the performance of models (dispersion of validation metrics), we calculated, for each species, the difference between the $AUC_{ROC}$ of each method and the average $AUC_{ROC}$ of all modelling methods for that species (Figure 1). Values higher than zero indicate $AUC_{ROC}$ higher than average.

```{r}
library(tidyverse)

cv_result <- read.csv("results/nceas_blockcv_total.csv")

cv_result <- filter(cv_result, !model %in% c("GAM-unweighted", "GLM-unweighted"))

# changing the names
cv_result$model <- ifelse(cv_result$model == "RF", "RF down-sample", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Ranger", "RF-shallow", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "GLM", "GLM-step", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Lasso", "GLM-lasso", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "MaxEnt", "MaxEnt (default)", cv_result$model)
cv_result <- filter(cv_result, !model %in% c("MaxEnt-tuned", "MaxEnt-spatial-tuned", 
                                           "MaxEnt-noclamp", "MaxEnt-H2",
                                           "MaxEnt-LQ", "MaxEnt-LQP"))

# ggplot(data = cv_result, aes(x = model, y = ROC)) +
#   geom_violin() + 
#   facet_wrap(~ cv)
morder <- c("SVM",
            "MARS",
            "RF-shallow",
            "GAM",
            "GLM-step",
            "GLM-lasso",
            "BRT",
            "RF down-sample",
            "MaxEnt (default)",
            "Ensemble")

```


```{r, fig.width=8.2, fig.height=4.2, fig.cap="Difference from average $AUC_{ROC}$."}
stat_fun <- function(x){
  mn <- mean(x, na.rm = TRUE)
  out <- x - mn
  return(out)
}

diff_mean_rcv <- cv_result %>% 
  filter(cv == "random") %>%
  dplyr::select(species, cv, model, ROC) %>% 
  pivot_wider(names_from = model, values_from = ROC) %>% 
  dplyr::select(-species, -cv) %>% 
  apply(MARGIN = 1, FUN = stat_fun) %>%
  as.data.frame() %>% 
  mutate(model = rownames(.)) %>% 
  pivot_longer(cols = 1:(ncol(.) - 1)) %>% 
  mutate(cv = "Random partitioning")
diff_mean_scv <- cv_result %>% 
  filter(cv == "spatial") %>%
  dplyr::select(species, cv, model, ROC) %>% 
  pivot_wider(names_from = model, values_from = ROC) %>% 
  dplyr::select(-species, -cv) %>% 
  apply(MARGIN = 1, FUN = stat_fun) %>%
  as.data.frame() %>% 
  mutate(model = rownames(.)) %>% 
  pivot_longer(cols = 1:(ncol(.) - 1)) %>% 
  mutate(cv = "Spatial partitioning")

diff_mean <- bind_rows(diff_mean_rcv, diff_mean_scv)

ggplot(data = diff_mean, aes(x = model, y = value, col = value)) +
  geom_jitter(alpha = 0.4, width = 0.3) +
  geom_violin(fill = NA) +
  # viridis::scale_colour_viridis(option = "A", direction = -1) +
  scale_colour_gradientn(colours = RColorBrewer::brewer.pal(11, "Spectral"),
                         limits = c(-0.2, 0.2)) +
  facet_wrap(~ cv) +
  # theme_minimal() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(limits = morder) +
  labs(y = "Difference from average", x = "", col = "Difference")

```


## `r sn <- sn + 1; sn`- Top rank methods in different evaluations

In the main text, the mean performance and the average rank of performance metrics (i.e., $AUC_{ROC}$, $AUC_{PRG}$, and COR) were used to assess and compare models. However, this approach does not show how frequently a method was the top method or whether it was among the top 2 or 3 methods. Here we calculated the percentage of species (171 species in our study) for which a method was in the top 1, top 2, or top 3 methods (Figure 2 and 3). For example, for $AUC_{ROC}$ and random partitioning (Figure 2), the Ensemble was within the top 3 methods for 63.7% of the species.

Notice that all methods performed the best (top 1) for at least a few species in both random and spatial partitioning. A noticeable result here is that while some methods like GLM-step, MARS, or RF-shallow are average performers overall, they may be top methods more frequently than better average performers (for example, GLM-step vs BRT in the top 1 methods for random partitioning; or RF-shallow vs RF down-sample or MaxEnt in top 1 method for spatial partitioning).

Another highlight is that although the Ensemble was not the top performer for many species, it was among the top 3 methods for more than half of them in both partitioning strategies (Figure 2 and 3).


```{r fig.height=4, fig.width=6, fig.cap="Top rank methods in random partitioning."}
library(tidyverse)

cv_result <- read.csv("results/nceas_blockcv_total.csv")

# changing the names
cv_result$model <- ifelse(cv_result$model == "RF", "RF down-sample", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Ranger", "RF-shallow", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "GLM", "GLM-step", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Lasso", "GLM-lasso", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "MaxEnt", "MaxEnt (default)", cv_result$model)

cv_result <- filter(cv_result, !model %in% c("GAM-unweighted", "GLM-unweighted"))
cv_result <- filter(cv_result, !model %in% c("MaxEnt-tuned", "MaxEnt-spatial-tuned", 
                                             "MaxEnt-noclamp", "MaxEnt-H2",
                                             "MaxEnt-LQ", "MaxEnt-LQP"))

morder <- c("SVM",
            "MARS",
            "RF-shallow",
            "GAM",
            "GLM-step",
            "GLM-lasso",
            "BRT",
            "RF down-sample",
            "MaxEnt (default)",
            "Ensemble")

lastn <- function(x, n = 3, by = "COR"){
  x1 <- unique(x[, by, drop = TRUE])
  idx <- tail(order(x1), n)
  xx <- x1[idx]
  idx2 <- which(x[, by, drop = TRUE] %in% xx)
  return(x[idx2, ])
}

topModFun <- function(evalmetric, nm){
  cv_result %>% 
    filter(cv == "random") %>% # filter by cv 
    group_by(species) %>% 
    nest() %>% 
    mutate(df = map(data, ~lastn(x = ., n = nm, by = evalmetric)),
           models = map(df, pluck("model"))) %>% 
    select(models) %>% 
    map(unlist) %>% 
    pluck("models") %>% 
    table() %>% 
    as.data.frame() %>% 
    setNames(c("Model", "freq")) %>% 
    mutate(prop = (freq / 171) * 100,
           percent = round(prop, 1),
           metric = evalmetric,
           position = paste("Top", nm))
}

topModels <- map2(rep(c("ROC", "PRG", "COR"), each = 3), 
                  rep(c(1,2,3), 3), 
                  topModFun) %>% 
  bind_rows()

# create labes for the facet's labeller
topModels$metlabel <- factor(topModels$metric, labels = c(
  'COR',
  '"AUC"["PRG"]',
  '"AUC"["ROC"]'
))
topModels$metlabel <- fct_relevel(topModels$metlabel, c(c('"AUC"["ROC"]',
                                                          '"AUC"["PRG"]',
                                                          'COR')))



# plot all the metrics
ggplot(data = topModels, aes(x = Model, y = forcats::fct_rev(as.factor(position)), fill = percent)) + 
  geom_tile(color = "gray") +
  facet_wrap(vars(metlabel), nrow = 3, strip.position = "right", labeller = label_parsed) +
  geom_text(aes(label = percent, colour = percent), size = 3.5) +
  viridis::scale_fill_viridis(option = "A", direction = -1) +
  viridis::scale_colour_viridis(option = "E", direction = 1, begin = 0.2, end = 0.8) +
  labs(x = NULL, y = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 8, angle = 60, hjust = 1),
        axis.text.y = element_text(size = 8),
        axis.title.y = element_text(margin = margin(r = 10)),
        text = element_text(size = 11, family = "Helvetica")) +
  guides(fill = "none", colour = "none") +
  scale_x_discrete(limits = morder) +
  ggtitle("Random partitioning")

```


```{r fig.height=4, fig.width=6, fig.cap="Top rank methods in spatial partitioning"}

topModFun <- function(evalmetric, nm){
  cv_result %>% 
    filter(cv == "spatial") %>% # filter by cv 
    group_by(species) %>% 
    nest() %>% 
    mutate(df = map(data, ~lastn(x = ., n = nm, by = evalmetric)),
           models = map(df, pluck("model"))) %>% 
    select(models) %>% 
    map(unlist) %>% 
    pluck("models") %>% 
    table() %>% 
    as.data.frame() %>% 
    setNames(c("Model", "freq")) %>% 
    mutate(prop = (freq / 171) * 100,
           percent = round(prop, 1),
           metric = evalmetric,
           position = paste("Top", nm))
}

topModels <- map2(rep(c("ROC", "PRG", "COR"), each = 3), 
                  rep(c(1,2,3), 3), 
                  topModFun) %>% 
  bind_rows()

# create labes for the facet's labeller
topModels$metlabel <- factor(topModels$metric, labels = c(
  'COR',
  '"AUC"["PRG"]',
  '"AUC"["ROC"]'
))
topModels$metlabel <- fct_relevel(topModels$metlabel, c(c('"AUC"["ROC"]',
                                                          '"AUC"["PRG"]',
                                                          'COR')))

# plot all the metrics
ggplot(data = topModels, aes(x = Model, y = forcats::fct_rev(as.factor(position)), fill = percent)) + 
  geom_tile(color = "gray") +
  facet_wrap(vars(metlabel), nrow = 3, strip.position = "right", labeller = label_parsed) +
  geom_text(aes(label = percent, colour = percent), size = 3.5) +
  viridis::scale_fill_viridis(option = "A", direction = -1) +
  viridis::scale_colour_viridis(option = "E", direction = 1, begin = 0.2, end = 0.8) +
  labs(x = NULL, y = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 8, angle = 60, hjust = 1),
        axis.text.y = element_text(size = 8),
        axis.title.y = element_text(margin = margin(r = 10)),
        text = element_text(size = 11, family = "Helvetica")) +
  guides(fill = "none", colour = "none") +
  scale_x_discrete(limits = morder) +
  ggtitle("Spatial partitioning")

```


An interesting result is that although Ensemble is the most frequent top performer in terms of $AUC_{ROC}$ and $AUC_{PRG}$ when predicting spatially separated testing data, and second best in terms of COR in spatial partitioning. Under random partitioning, Ensemble was the second or third best average performer. We explored further how Ensemble is performing compared to its component models in the following section.

## `r sn <- sn + 1; sn`- Rank of the Ensemble vs its components

Here we calculated the same plots but only for the Ensemble and its component models i.e, GLM-lasso, GAM, MaxEnt (default), BRT, and RF down-sample (Figure 4). The Ensemble appeared in the top 2 and 3 models for more species than its component in both random and spatial partitioning. For top 1 methods, it was only best for $AUC_{ROC}$ (along with RF down-sample) in random partitioning, but the best for both AUCs and the second-best for COR in spatial partitioning. The fact that Ensemble appears better than its component in spatial partitioning may be evidence that ensembling of tuned models can lead to better generalisation.

```{r}
ensmodels <- c(
  "GAM",
  "GLM-lasso",
  "BRT",
  "RF down-sample",
  "MaxEnt (default)",
  "Ensemble"
)
cv_result <- cv_result %>%
  filter(model %in% ensmodels)

topModFun <- function(evalmetric, nm){
  cv_result %>% 
    filter(cv == "random") %>% # filter by cv 
    group_by(species) %>% 
    nest() %>% 
    mutate(df = map(data, ~lastn(x = ., n = nm, by = evalmetric)),
           models = map(df, pluck("model"))) %>% 
    select(models) %>% 
    map(unlist) %>% 
    pluck("models") %>% 
    table() %>% 
    as.data.frame() %>% 
    setNames(c("Model", "freq")) %>% 
    mutate(prop = (freq / 171) * 100,
           percent = round(prop, 1),
           metric = evalmetric,
           position = paste("Top", nm))
}

topModels <- map2(rep(c("ROC", "PRG", "COR"), each = 3), 
                  rep(c(1,2,3), 3), 
                  topModFun) %>% 
  bind_rows()

# create labes for the facet's labeller
topModels$metlabel <- factor(topModels$metric, labels = c(
  'COR',
  '"AUC"["PRG"]',
  '"AUC"["ROC"]'
))
topModels$metlabel <- fct_relevel(topModels$metlabel, c(c('"AUC"["ROC"]',
                                                          '"AUC"["PRG"]',
                                                          'COR')))

# plot all the metrics
p1 <- ggplot(data = topModels, aes(x = Model, y = forcats::fct_rev(as.factor(position)), fill = percent)) + 
  geom_tile(color = "gray") +
  facet_wrap(vars(metlabel), nrow = 3, strip.position = "right", labeller = label_parsed) +
  geom_text(aes(label = percent, colour = percent), size = 3.5) +
  viridis::scale_fill_viridis(option = "A", direction = -1) +
  viridis::scale_colour_viridis(option = "E", direction = 1, begin = 0.2, end = 0.8) +
  labs(x = NULL, y = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 8, angle = 60, hjust = 1),
        axis.text.y = element_text(size = 8),
        axis.title.y = element_text(margin = margin(r = 10)),
        text = element_text(size = 11, family = "Helvetica")) +
  guides(fill = "none", colour = "none") +
  scale_x_discrete(limits = ensmodels) +
  ggtitle("Random partitioning")

```

```{r}
topModFun <- function(evalmetric, nm){
  cv_result %>% 
    filter(cv == "spatial") %>% # filter by cv 
    group_by(species) %>% 
    nest() %>% 
    mutate(df = map(data, ~lastn(x = ., n = nm, by = evalmetric)),
           models = map(df, pluck("model"))) %>% 
    select(models) %>% 
    map(unlist) %>% 
    pluck("models") %>% 
    table() %>% 
    as.data.frame() %>% 
    setNames(c("Model", "freq")) %>% 
    mutate(prop = (freq / 171) * 100,
           percent = round(prop, 1),
           metric = evalmetric,
           position = paste("Top", nm))
}

topModels <- map2(rep(c("ROC", "PRG", "COR"), each = 3), 
                  rep(c(1,2,3), 3), 
                  topModFun) %>% 
  bind_rows()

# create labes for the facet's labeller
topModels$metlabel <- factor(topModels$metric, labels = c(
  'COR',
  '"AUC"["PRG"]',
  '"AUC"["ROC"]'
))
topModels$metlabel <- fct_relevel(topModels$metlabel, c(c('"AUC"["ROC"]',
                                                          '"AUC"["PRG"]',
                                                          'COR')))

# plot all the metrics
p2 <- ggplot(data = topModels, aes(x = Model, y = forcats::fct_rev(as.factor(position)), fill = percent)) + 
  geom_tile(color = "gray") +
  facet_wrap(vars(metlabel), nrow = 3, strip.position = "right", labeller = label_parsed) +
  geom_text(aes(label = percent, colour = percent), size = 3.5) +
  viridis::scale_fill_viridis(option = "A", direction = -1) +
  viridis::scale_colour_viridis(option = "E", direction = 1, begin = 0.2, end = 0.8) +
  labs(x = NULL, y = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 8, angle = 60, hjust = 1),
        axis.text.y = element_text(size = 8),
        axis.title.y = element_text(margin = margin(r = 10)),
        text = element_text(size = 11, family = "Helvetica")) +
  guides(fill = "none", colour = "none") +
  scale_x_discrete(limits = ensmodels) +
  ggtitle("Spatial partitioning")

```

```{r fig.height=5, fig.width=7.5, fig.cap="Top rank methods among Ensemble and its components."}
cowplot::plot_grid(p1, p2)

```


To assess whether the Ensemble improves with respect to the best method in the set or not, we further explored (Figure 5) and realised that in all cases of “Top 1”, Ensemble actually somewhat outperformed its component (rather than being just as good as the best of its components). This could be an indication that the ensemble gains by combining “complementary” predictions.

```{r fig.height=5.5, fig.width=4.5, fig.cap="The differnce between $AUC_{ROC}$ of Ensemble and $AUC_{ROC}$ of the best component method."}

cv_result <- read.csv("results/nceas_blockcv_total.csv")

cv_result <- filter(cv_result, !model %in% c("GAM-unweighted", "GLM-unweighted"))

# changing the names
cv_result$model <- ifelse(cv_result$model == "RF", "RF down-sample", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Ranger", "RF-shallow", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "GLM", "GLM-step", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Lasso", "GLM-lasso", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "MaxEnt", "MaxEnt (default)", cv_result$model)
cv_result <- filter(cv_result, !model %in% c("MaxEnt-tuned", "MaxEnt-spatial-tuned", 
                                           "MaxEnt-noclamp", "MaxEnt-H2",
                                           "MaxEnt-LQ", "MaxEnt-LQP"))

topmodels <- cv_result %>% 
  filter(model %in% c("GAM", "GLM-lasso", "MaxEnt (default)",
                      "BRT", "RF down-sample", "Ensemble")) %>%
  # filter(cv == "random") %>%
  dplyr::select(species, model, cv, COR) %>% 
  pivot_wider(names_from = model, values_from = COR) %>% 
  relocate(species, cv, Ensemble)

topmodels$top <- topmodels %>% 
  dplyr::select(-species, -Ensemble, -cv) %>% 
  apply(1, max)

ens_gain <- topmodels %>% 
  mutate(top = as.numeric(top),
         enhigh = Ensemble - top) # %>% as.data.frame() %>% head()
  # group_by(cv) %>%
  # summarise(sum(enhigh) / 171 * 100)
  # pull(enhigh) %>%
  # sum() / 171 * 100

ggplot(data = ens_gain, aes(paste(str_to_title(cv), "partitioning"), enhigh, col = enhigh)) +
  geom_violin(fill = NA) +
  geom_jitter(width = 0.1, size = 2.5, alpha = 0.8) +
  scale_colour_gradientn(colours = RColorBrewer::brewer.pal(11, "Spectral"),
                         limits = c(-0.12, 0.12),
                         breaks = c(-0.1, 0, 0.1)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "", y = "Difference from the top component method", 
       col = "AUC difference") 

```


## `r sn <- sn + 1; sn`- Interactions in BRT and MaxEnt

To assess the impacts of interactions in flexible methods, we compared two of our top methods, BRT and MaxEnt, with and without interactions (Figure 6). We implemented a BRT model with a set tree-complexity of `1` also known as **stump**. The main difference between the BRT-stump and the default BRT (with tree-complexity 1 or 5) is that the default BRT is allowed to fit a higher level of interaction between the covariates if there are more than 50 species records in the training data (57% of cases). The implemented BRT method (Elith *et al.* 2008) utilized internal cross-validation to find the best number of trees for the model. Thus, by limiting tree-complexity to 1, the model adds more trees to find a similar balance in the fitted model as the BRT with tree-complexity 5.   

MaxEnt is also presented as two variants here, the MaxEnt with enforced LQ features and one with enforced LQP features. The main difference between these two is that MaxEnt-LQP accommodates interaction as the product of the linear features.


**The main BRT was modelled with a tree-complexity of 1 (stump) for 43% of the times.**

```{r fig.width=8, fig.height=4, fig.cap="Perfromance of implementation of BRT and MaxEnt with forced limited flexibility."}
# BRT and Maxent second run -----------------------------------------------
cv_result <- read.csv("results/nceas_blockcv_total.csv")
brtmxnt <- read.csv("results/nceas_blockcv_brt_maxent.csv")

# changing the names
cv_result$model <- ifelse(cv_result$model == "RF", "RF down-sample", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Ranger", "RF-shallow", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "GLM", "GLM-step", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Lasso", "GLM-lasso", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "MaxEnt", "MaxEnt (default)", cv_result$model)

cv_result <- filter(cv_result, !model %in% c("GAM-unweighted", "GLM-unweighted"))
cv_result <- filter(cv_result, !model %in% c("MaxEnt-tuned", "MaxEnt-spatial-tuned", 
                                             "MaxEnt-noclamp", "MaxEnt-H2",
                                             "MaxEnt-LQP"))

cv_result <- bind_rows(cv_result, brtmxnt)

mean_cv2 <-  cv_result %>%
  group_by(model, cv) %>% 
  summarise(
    ROC_mean = mean(ROC), ROC_se = 1 * (sd(ROC) / sqrt(n())),
    PRG_mean = mean(PRG), PRG_se = 1 * (sd(PRG) / sqrt(n())),
    COR_mean = mean(COR, na.rm = TRUE), COR_se = 1 * (sd(COR, na.rm = TRUE) / sqrt(n()))
  )


cols2 <- c(
  "SVM" = "gray90",
  "Ensemble" = "gray90",
  
  "MARS" = "gray90",
  "GLM-step" = "gray90",
  "GLM-lasso" = "gray90",
  "GAM" = "gray90",
  "MaxEnt (default)" = "gray90",
  "MaxEnt-LQP" = "#BE2207",
  "MaxEnt-LQ" = "#35274A",
  
  "RF-shallow" = "gray90",
  "BRT" = "#046C9A",
  "BRT-stump" = "#0B775E",
  "RF down-sample" = "gray90"
)


mean_cv3 <- filter(mean_cv2, model %in% c("BRT", 
                                          "BRT-stump",
                                          "MaxEnt-LQ",
                                          "MaxEnt-LQP"))
mean_cv2 <- filter(mean_cv2, !model %in% c("BRT", 
                                           "BRT-stump",
                                           "MaxEnt-LQ",
                                           "MaxEnt-LQP"))

ggplot(data = mean_cv2, aes(x = ROC_mean, y = COR_mean, color = model)) +
  scale_color_manual(values = cols2) +
  geom_segment(aes(x = ROC_mean, 
                   y = COR_mean - COR_se, 
                   xend = ROC_mean, 
                   yend = COR_mean + COR_se,
                   colour = model),
               alpha = 0.8, data = mean_cv2) +
  geom_segment(aes(x = ROC_mean - ROC_se, 
                   y = COR_mean, 
                   xend = ROC_mean + ROC_se, 
                   yend = COR_mean, 
                   colour = model), 
               alpha = 0.8, data = mean_cv2) +
  geom_point(size = 2) +
  ggrepel::geom_text_repel(aes(x = ROC_mean, 
                               y = COR_mean,
                               colour = model,
                               label = model),
                           force = 5,
                           data = mean_cv2) +
  geom_segment(aes(x = ROC_mean, 
                   y = COR_mean - COR_se, 
                   xend = ROC_mean, 
                   yend = COR_mean + COR_se),
               colour = "gray50",
               alpha = 0.8, data = mean_cv3) +
  geom_segment(aes(x = ROC_mean - ROC_se, 
                   y = COR_mean, 
                   xend = ROC_mean + ROC_se, 
                   yend = COR_mean), 
               colour = "gray50", 
               alpha = 0.8, data = mean_cv3) +
  geom_point(size = 2, data = mean_cv3, aes(x = ROC_mean, y = COR_mean, color = model)) +
  ggrepel::geom_text_repel(aes(x = ROC_mean, 
                               y = COR_mean,
                               colour = model,
                               label = model),
                           force = 5,
                           data = mean_cv3) +
  facet_wrap(~ paste(str_to_title(cv), "partitioning")) +
  labs(x = expression("AUC"["ROC"]), y = "COR") +
  theme_bw(base_line_size = 0.2) +
  theme(text = element_text(size = 12, family = "Helvetica"),
        legend.position = "none") +
  scale_x_continuous(breaks = seq(from = 0.67, to = 0.75, by = 0.02))


```

There is a small and non significant difference between the average performance of the BRT vs BRT-stump and between MaxEnt-LQP vs MaxEnt-LQ. The MaxEnt-LQP had a lower mean $AUC_{ROC}$ and COR in both partitioning methods which could be because enforced interaction in the model is too complex for some species with a very low number of presence records. On the other hand, BRT performed better than BRT-stump, implying that the additional flexibility of interactions is beneficial.  


## `r sn <- sn + 1; sn`- Extrapolation in testing blocks

It is useful to know whether extrapolation occurs when models are used to predict to spatially separated points. Extrapolation occurs when the testing/predicting sites have environmental values outside of the range of environmental conditions used in the training samples. To measure the amount of extrapolation in testing sites we used Multivariate Environmental Similarity Surface (**MESS**) introduced by Elith *et al.* (2010). We modified the `mess` function in the **dismo** R package to compute MESS values for points, not rasters. We estimated MESS values for the records in the presence-absence evaluation dataset, using the training presence-TGB as the reference sites. We used only continuous covariates for this. For more explanation on MESS, see Elith *et al.* (2010).    

In Figure 7, we summed the number of testing points with extrapolation (negative MESS values) for each species. This gives a good sense of how frequently species from a region experience extrapolation when predicting.


```{r }
library(tidyverse)

messall <- read_csv("results/nceas_blockcv_mess.csv")

```


```{r, fig.width=8, fig.height=4, fig.cap="The sum of extrapolated points of each species in each region."}
messall %>% 
  mutate(cv = paste(str_to_title(cv), "partitioning")) %>% 
ggplot(data = ., aes(x = region, y = extrapolate)) +
  geom_boxplot() +
  geom_jitter(alpha = 0.3) +
  facet_wrap(cv ~ .) +
  # theme_bw() +
  labs(x = "Regions", y = "Number of extraploated points")

```


Figure 8 shows the number of extrapolated sites when using spatial partitioning compared to random partitioning.

```{r, fig.width=5, fig.height=4, fig.cap="The sum of the number of extrapolated points in each species in random vs spatial partitioning."}
messall %>% 
  pivot_wider(names_from = cv, values_from = extrapolate) %>% 
ggplot(data = ., aes(y = spatial, x = random, col = region)) +
  geom_point(size = 2.5, alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, color = "gray") +
  scale_x_continuous(limits = c(0, 43)) +
  scale_y_continuous(limits = c(0, 43)) +
  coord_equal() + 
  theme_classic() +
  labs(x = "Random partitioning", y = "Spatial partitioning", col = "Regions")

```


## `r sn <- sn + 1; sn`- List of species used in modelling

For creating spatial blocks some species did not have enough data to fit and evaluate models. Here, we provide the list of species we used in our study (Table 9). You can see the location of each region on the world map in Figure 9. Read detailed explanation of this dataset in Elith *et al* (2020).

```{r fig.width=14, fig.height=8, fig.cap="Location of each region in the world."}
library(tidyverse)
library(ggrepel)
library(sf)
library(disdat)

basemap <- geodata::world(resolution = 5, path = "temp/") %>% 
  st_as_sf()

robinson <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"


regions <- c('AWT', 'NSW', 'CAN', 'SWI', 'NZ', 'SA')
reg <- regions[2]

xycent <- data.frame("x" = 1, "y" = 1, 'id' = 'aa')
for(i in seq_along(regions)){
  reg <- regions[i]
  map <- disdat::disBorder(reg) %>% 
    st_union() %>% 
    st_transform(crs = robinson)
  pnt <- st_centroid(map)
  xycent[i, 'x'] <- st_coordinates(pnt)[1]
  xycent[i, 'y'] <- st_coordinates(pnt)[2]
  xycent[i, 'id'] <- reg
}

xycent$names <- c('Australian Wet Tropic (AWT)',
                  'New South Wales (NSW)',
                  'Ontario, Canada (CAN)',
                  'Switzerland (SWI)',
                  'New Zealand (NZ)',
                  'South American countries (SA)')

# Swiss location is not correct in the original data
swiss_xy <- data.frame(x=8, y=47) %>% 
  st_as_sf(coords=1:2, crs=4326) %>% 
  st_transform(crs = robinson) %>% 
  st_coordinates()

xycent$x[4] <- swiss_xy[1]
xycent$y[4] <- swiss_xy[2]

ggplot() +
  geom_sf(data = basemap, col = 'white', size = 0.1, fill = 'gray') +
  geom_point(data = xycent, aes(x=x, y=y), color = 'red', size = 4) +
  geom_text_repel(data = xycent, aes(x=x, y=y, label = names), size = 7,
                  box.padding = 0.8, max.overlaps = Inf) +
  coord_sf(crs = robinson) +
  theme_minimal() +
  theme(axis.title =  element_blank())

```



```{r}
sp_list <- read.csv("species_list.csv")
names(sp_list)[1] <- ""

knitr::kable(sp_list, caption = "List of species used for modelling. PO is the number of presence-only recods in the tarining dataset, TGBs is the number of Target-Group-Background samples, and Presence/Absence are the number of records in the evaluation dataset.") 


```


## `r sn <- sn + 1; sn`- Statistical tests
### `r sn <- sn + 0.1; sn`- Statistical test for random partitioning

Here the statistical test on the differences between methods in random partitioning are presented (Figure 10). The plots are $AUC_{ROC}$, $AUC_{PRG}$, and COR from top to bottom, respectively. The number on the top of the x-axis shows the range of the ranks of the models. The average rank of each model is indicated by the thin line connected to the axis. The lines (methods) that are connected by the horizontal thick line are not statistically different at 0.05 significance level.

```{r fig.width=9, fig.cap="Average rank and statistical difference of the methods in random partitioning."}
library(tidyverse)
library(scmamp)

cv_result <- read.csv("results/nceas_blockcv_total.csv")

# changing the names
cv_result$model <- ifelse(cv_result$model == "RF", "RF down-sample", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Ranger", "RF-shallow", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "GLM", "GLM-step", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "Lasso", "GLM-lasso", cv_result$model)
cv_result$model <- ifelse(cv_result$model == "MaxEnt", "MaxEnt-default", cv_result$model)

cv_result <- filter(cv_result, !model %in% c("GAM-unweighted", "GLM-unweighted"))

cv_result2 <- filter(cv_result, !model %in% c("MaxEnt-tuned", "MaxEnt-spatial-tuned", 
                                              "MaxEnt-noclamp", "MaxEnt-H2",
                                              "MaxEnt-LQ", "MaxEnt-LQP"))

################################
auc_ranks1 <- cv_result2 %>% 
  filter(cv == "random") %>% 
  pivot_wider(id_cols = species, names_from = model, values_from = ROC) %>% 
  dplyr::select(- species) %>% 
  as.matrix()

post_results1 <- postHocTest(
  data = auc_ranks1,
  test = "aligned ranks",
  correct = "shaffer",
  control = NULL,
  use.rank = TRUE
)

colnames(post_results1$summary) <- gsub("F.", "F ", colnames(post_results1$summary))
colnames(post_results1$summary) <- gsub("down.sample", "down-sample", colnames(post_results1$summary))
colnames(post_results1$corrected.pval) <- gsub("F.", "F ", colnames(post_results1$summary))
rownames(post_results1$corrected.pval) <- gsub("down.sample", "down-sample", colnames(post_results1$summary))
colnames(post_results1$corrected.pval) <- gsub("F.", "F ", colnames(post_results1$summary))
rownames(post_results1$corrected.pval) <- gsub("down.sample", "down-sample", colnames(post_results1$summary))

################################
auc_ranks2 <- cv_result2 %>% 
  filter(cv == "random") %>% 
  pivot_wider(id_cols = species, names_from = model, values_from = PRG) %>% 
  dplyr::select(- species) %>% 
  as.matrix()

post_results2 <- postHocTest(
  data = auc_ranks2,
  test = "aligned ranks",
  correct = "shaffer",
  control = NULL,
  use.rank = TRUE
)

colnames(post_results2$summary) <- gsub("F.", "F ", colnames(post_results1$summary))
colnames(post_results2$summary) <- gsub("down.sample", "down-sample", colnames(post_results1$summary))
colnames(post_results2$corrected.pval) <- gsub("F.", "F ", colnames(post_results1$summary))
rownames(post_results2$corrected.pval) <- gsub("down.sample", "down-sample", colnames(post_results1$summary))
colnames(post_results2$corrected.pval) <- gsub("F.", "F ", colnames(post_results1$summary))
rownames(post_results2$corrected.pval) <- gsub("down.sample", "down-sample", colnames(post_results1$summary))

################################
auc_ranks3 <- cv_result2 %>% 
  filter(cv == "random") %>% 
  pivot_wider(id_cols = species, names_from = model, values_from = COR) %>% 
  dplyr::select(- species) %>% 
  as.matrix()

post_results3 <- postHocTest(
  data = auc_ranks3,
  test = "aligned ranks",
  correct = "shaffer",
  control = NULL,
  use.rank = TRUE
)

colnames(post_results3$summary) <- gsub("F.", "F ", colnames(post_results3$summary))
colnames(post_results3$summary) <- gsub("down.sample", "down-sample", colnames(post_results3$summary))
colnames(post_results3$corrected.pval) <- gsub("F.", "F ", colnames(post_results3$summary))
rownames(post_results3$corrected.pval) <- gsub("down.sample", "down-sample", colnames(post_results3$summary))
colnames(post_results3$corrected.pval) <- gsub("F.", "F ", colnames(post_results3$summary))
rownames(post_results3$corrected.pval) <- gsub("down.sample", "down-sample", colnames(post_results3$summary))

################################
# c(bottom, left, top, right)
par(mfrow=c(3,1), mai = c(0, 0, 0, 0))
plotRanking(post_results1$corrected.pval, 
            post_results1$summary,
            alpha = 0.05, 
            cex = 1.3,
            decreasing = FALSE)

plotRanking(post_results2$corrected.pval, 
            post_results2$summary,
            alpha = 0.05, 
            cex = 1.3,
            decreasing = FALSE)

plotRanking(post_results3$corrected.pval, 
            post_results3$summary,
            alpha = 0.05, 
            cex = 1.3,
            decreasing = FALSE)

```

<!-- \newpage -->

### `r sn <- sn + 0.1; sn`- Statistical test for MaxEnt variants in spatial partitioning

The Friedman’s Aligned Rank test indicated no significant difference among MaxEnt variants for $AUC_{PRG}$. Thus here we only plot the result of pairwise test for $AUC_{ROC}$ and COR (Figure 11). In general the differences between models are also insignificant for these statistics.


```{r}
library(tidyverse)
library(scmamp)

cv_result <- read.csv("results/nceas_blockcv_total.csv")

cv_result2 <- filter(cv_result, model %in% c("MaxEnt", 
                                             "MaxEnt-tuned", 
                                             "MaxEnt-spatial-tuned", 
                                             "MaxEnt-noclamp", 
                                             "MaxEnt-LQ"))

cv_result2$model <- ifelse(cv_result2$model == "MaxEnt", "MaxEnt (default)", cv_result2$model)

################################
roc_maxents <- cv_result2 %>% 
  filter(cv == "spatial") %>% 
  pivot_wider(id_cols = species, names_from = model, values_from = ROC) %>% 
  dplyr::select(- species) %>% 
  as.matrix()

friedmanAlignedRanksTest(roc_maxents)

post_results1 <- postHocTest(
  data = roc_maxents,
  test = "aligned ranks",
  correct = "shaffer",
  control = NULL,
  use.rank = TRUE
)


################################
prg_maxents <- cv_result2 %>% 
  filter(cv == "spatial") %>% 
  pivot_wider(id_cols = species, names_from = model, values_from = PRG) %>% 
  dplyr::select(- species) %>% 
  as.matrix()

friedmanAlignedRanksTest(prg_maxents)

post_results2 <- postHocTest(
  data = prg_maxents,
  test = "aligned ranks",
  correct = "shaffer",
  control = NULL,
  use.rank = TRUE
)


################################
cor_maxents <- cv_result2 %>% 
  filter(cv == "spatial") %>% 
  pivot_wider(id_cols = species, names_from = model, values_from = COR) %>% 
  dplyr::select(- species) %>% 
  as.matrix()

friedmanAlignedRanksTest(cor_maxents)

post_results3 <- postHocTest(
  data = cor_maxents,
  test = "aligned ranks",
  correct = "shaffer",
  control = NULL,
  use.rank = TRUE
)

```
    

```{r fig.cap="Average rank and statistical difference of the MaxEnt variants in spatial partitioning."}
# c(bottom, left, top, right)
par(mfrow=c(3,1), mai = c(0, 0, 0, 0))
plotRanking(post_results1$corrected.pval, 
            post_results1$summary,
            alpha = 0.05, 
            cex = 1.3,
            decreasing = FALSE)
# 
# plotRanking(post_results2$corrected.pval, 
#             post_results2$summary,
#             alpha = 0.05, 
#             cex = 1.3,
#             decreasing = FALSE)

plotRanking(post_results3$corrected.pval, 
            post_results3$summary,
            alpha = 0.05, 
            cex = 1.3,
            decreasing = FALSE)
```



## `r sn <- floor(sn) + 1; sn`- Explanation for poor performance of some flexible methods

One result of interest is that some flexible methods (SVM, MARS and RF-shallow) did not perform particularly well (Fig. 3-5 in the main text). So why do these perform less well? It could be that they are not flexible enough: they have no or very limited interactions (Table 1 in the main text). To explore this idea, we did explore how much interactions affect results for BRT and MaxEnt (Section 8). Impacts are small and variable – for BRT, allowing interactions for species with `> 50` data points improves performance slightly, and for MaxEnt, enforcing product features decreases it, probably because enforced feature classes perform less well on species with fewer records. Other issues may be at play for SVM, MARS and RF-shallow. It is possible that our implementation of them was not optimal, since we have more experience with the methods that did well than with some of these. SVM was fitted with defaults. Tuning options are available, so it may be worth testing this further. SVM did not perform as well in the “top 3” analysis (see Section 6) as it did in Valavi *et al.* (2022), where it was fitted with the same settings but on presence-random background data (rather than presence-TGB), and using more species and more data per species (because the current test-train experiment requires subsampling). These nuances could be further explored. MARS was fitted with the [`earth`](https://cran.r-project.org/package=earth) package in R. In Elith *et al.* (2006) MARS, fitted there using the [`mda`](https://cran.r-project.org/package=mda) R package, performed relatively better than in Valavi *et al.* (2022) and this study. We used `earth` because `mda` failed to fit MARS with many background samples for many species.  Perhaps different implementations of this model would perform better. Finally, RF-shallow performed better on spatial partitioning than random partitioning (Figs 3-5 in the main text, and Section 6 here), and may benefit from some tuning of tree depth parameter (e.g., Valavi *et al.* 2021). This is an experimental approach for RF, and still needs to be explored further. 



## `r sn <- floor(sn) + 1; sn`- References

* Elith, J., Graham, C.H., Anderson, R.P., Dudík, M., Ferrier, S., Guisan, A., J Hijmans, R., Huettmann, F., R Leathwick, J., Lehmann, A., Li, J., G Lohmann, L., Overton, J.McC.M., Peterson, A.T., Phillips, S.J., Richardson, K., Scachetti-Pereira, R., Schapire, R.E., Soberón, J., Williams, S., Wisz, M.S. & Zimmermann, N.E. (2006) Novel methods improve prediction of species’ distributions from occurrence data. *Ecography*, 29, 129–151.

* Elith, J., Graham, C., Valavi, R., Abegg, M., Bruce, C., Ferrier, S., Ford, A., Guisan, A., Hijmans, R.J., Huettmann, F., Lohmann, L., Loiselle, B., Moritz, C., Overton, J., Peterson, A.T., Phillips, S., Richardson, K., Williams, S., Wiser, S.K., Wohlgemuth, T. & Zimmermann, N.E. (2020) Presence-only and Presence-absence Data for Comparing Species Distribution Modeling Methods. *Biodiversity Informatics*, 15, 69–80.

* Elith, J., Kearney, M. & Phillips, S. (2010) The art of modelling range-shifting species. *Methods in ecology and evolution*, 1, 330–342.

* Elith, J., Leathwick, J.R. & Hastie, T. (2008) A Working Guide to Boosted Regression Trees. *Journal of Animal Ecology*, 77, 802–813.

* Flach, P. & Kull, M. (2015) Precision-Recall-Gain Curves: PR Analysis Done Right. *Advances in Neural Information Processing Systems*, pp. 838–846. Curran Associates, Inc.

* Pearce, J. & Ferrier, S. (2000) Evaluating the predictive performance of habitat models developed using logistic regression. *Ecological modelling*, 133, 225–245.

* Phillips, S.J., Dudík, M., Elith, J., Graham, C.H., Lehmann, A., Leathwick, J. & Ferrier, S. (2009) Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. *Ecological applications*, 19, 181–197.

* Valavi, R., Elith, J., Lahoz‐Monfort, J.J. & Guillera‐Arroita, G. (2021) Modelling species presence‐only data with random forests. *Ecography*, 44, 1731–1742.

* Valavi, R., Guillera‐Arroita, G., Lahoz‐Monfort, J.J. & Elith, J. (2022) Predictive performance of presence‐only species distribution models: a benchmark study with reproducible code. *Ecological Monographs*, 92, 1–27.