Ensemble ecological niche modelling
Presence-pseudoabsence ensemble ecological niche models were generated
for each species based on the current climatic data using the BIOMOD
framework (Thuiller et al., 2009), implemented in the biomod2package (Thuiller et al., 2021) in R. The ensemble included four
algorithms known to be robust and perform well across a range of
distribution scales (Meller et al., 2014): multivariate adaptive
regression splines (MARS), artificial neural networks (ANN), random
forests (RF), and maximum entropy (MAXENT, implemented usingmaxnet ; Phillips et al., 2017). This diversity in computational
models provides low inter-correlation in model components, which is
ideal for higher predictive performance in ensemble ENMs (Elith, 2019;
Valavi et al., 2022). The models were calibrated according to the
default parameters and parameterization processes in biomod2 .
To validate model results we used five-fold cross-validation, dividing
the occurrence data into five subsets, and testing each subset against a
model calibrated on the remaining four subsets. Hold-out validation,
where a small subset of the data is used only for evaluation, while the
remaining data are used for calibration-validation, was unsuitable due
to the small sample size of occurrences in several species. Model
performance was evaluated using the area under the receiver operating
characteristic (ROC) curve (AUC) and the true skill statistic (TSS).
Variable contributions were calculated by performing five permutations
for each algorithm, each involving removal of the focal variable from
the model to calculate the difference in performance. We report the
average permutation importance.
The individual models with a final TSS value > 0.7 were
then combined into ensemble consensus models for each species which we
used with future climate data to predict the future climatically
suitable areas. Predicted suitability was reclassified into binary
suitable-unsuitable predictions using the threshold at which TSS was
maximised (values above this threshold represent climatic suitability).
These binary maps were used to characterise changes in climatically
suitable areas from current to future climatic scenarios using four
metrics: the relative change in climatically suitable areas (% change),
the percentage of the current area retained into the future, and the
distance and azimuth (angle of direction) between centroids of current
and future suitable areas. Change and retention were calculated
respectively as:
\begin{equation}
Change\ =\ \frac{Current\ suitable\ area\ -\ Future\ suitable\ area}{\text{Current\ suitable\ area}}\ *\ 100\nonumber \\
\end{equation}\begin{equation}
Retention\ =\ \frac{Current\ suitable\ area\ \land\ Future\ suitable\ area}{\text{Current\ suitable\ area}}\ *\ 100\nonumber \\
\end{equation}To estimate distance and azimuth we first defined distinct fragments
within suitable areas. Fragments represent connected suitable cells from
the binary maps, estimated by first applying morphological dilation
using a buffer of radius 2.5 km, followed by a morphological erosion
using a negative buffer of the same radius. As a result of this process,
suitable areas within 2.5 km of each other were connected to reflect the
assumption that bats can move easily within this distance. The resulting
number of fragments varied among species depending on the configuration
of the suitable area. For all fragments we then identified centroids and
the shortest straight paths between each centroid in the current
scenario and its nearest neighbour in the future scenario, assuming zero
movement costs. The length and azimuth of that path was used to estimate
distances and angles of direction. For each species we calculated a mean
distance and mean azimuth. Azimuths were classified into cardinal and
intercardinal directions based on angle, assuming North to be 0 and 360.
All centroid-based calculations were performed using thegeosphere package (Hijmans, 2022) in R.
Summary models of climatically suitable areas were generated for all
species together by creating a map of species richness for each time
period. All binary models for each time period were summed together into
a model that ranged from 0 to 110, representing the number of species
for which a cell is climatically suitable according to their binary
models. This model was then rescaled to a 0 to 1 scale for consistent
comparison and converted into a binary map using a threshold of 0.3,
such that in the final binary models for current and future, the
positive cells represented suitability hotspots - areas that were
climatically suitable for ≥ 30% (more than 33) of the study species.
Certainty in projected models was calculated as a continuous model by
first averaging individual replicates within a species to assess model
agreement, and then averaging these species models into a final model
ranging from 0 to 1, which can be interpreted as high certainty of
climatic unsuitability through uncertainty to high certainty of climatic
suitability.