# Clustering Future Scenarios Based on Predicted Range Maps

Davidow, Matthew  Merow, Cory  Che-Castaldo, Judy  Schafer, Toryn L. J  Düker,
Marie-Christine  Corcoran, Derek  Matteson, David S

(November 2020)

###### Summary

  1. 1.

Predictions of biodiversity trajectories under climate change are crucial in
order to act effectively in maintaining the diversity of species. In many
ecological applications, future predictions are made under various global
warming scenarios as described by a range of different climate models. The
outputs of these various predictions call for a reliable interpretation.

  2. 2.

We propose a interpretable and flexible two step methodology to measure the
similarity between predicted species range maps and cluster the future
scenario predictions utilizing a spectral clustering technique.

  3. 3.

We find that clustering based on ecological impact (predicted species range
maps) is mainly driven by the amount of warming. We contrast this with
clustering based only on predicted climate features, which is driven mainly by
climate models.

  4. 4.

The differences between these clusterings illustrate that it is crucial to
incorporate ecological information to understand the relevant differences
between climate models. The findings of this work can be used to better
synthesize forecasts of biodiversity loss under the wide spectrum of results
that emerge when considering potential future biodiversity loss.

Key-words: biodiversity; clustering; similarity measures; future scenarios;
animal species; climate change.

##  1 Introduction

Predicting ecological responses to a rapidly changing climate is essential to
enact effective conservation policies (Parry et al., 2007; Hannah et al.,
2013, 2020). Future range maps of various species can be predicted based on
predicted future patterns of climate (Burrows et al., 2014; Jones and Cheung,
2015; Molinos et al., 2016). However, predicting future patterns of climate is
a challenging task due to the many sources of uncertainty, and a plethora of
climate predictions are possible that are consistent with various unknown
factors (e.g. differing human policy responses, climate models). We are
interested in the analysis of comparing the outputs of these various
predictions and deducing common patterns among predictions.

We make use of climate predictions from the Coupled Model Intercomparison
Project 6 (CMIP6). This project provides multiple climate predictions which
vary by the underlying global climate model (GCM) and what representative
concentration pathway (RCP) is used. An RCP is a greenhouse gas concentration
trajectory (Eyring et al., 2016). Four such trajectories are included in
CMIP6, varying in the quantity of greenhouse gas emissions to capture the
uncertainty of future emissions. We refer to the RCP trajectory of least
emission as “optimistic”, and refer to the most pessimistic trajectory as the
“extreme” scenario. In this work we refer to a _scenario_ as a (GCM, RCP)
pair, these scenarios represent uncertainty both in the evolution of climate,
and in future greenhouse gas emissions.

In order to interpret the differences among climate predictions, we propose a
methodology to cluster the scenarios. We create such a clustering both from
the climate features, and from predicted range maps for 1101 mammalian
species. These predicted range maps are based on the predicted climate
features, the details of these range maps are discussed in Section 2. These
clusterings reveal the important differences between climate models, such as
whether the global climate model or the RCP differences are the most salient
features. The clustering will also group the scenarios into interpretable
collections, such as an “optimistic” collection of scenarios of lesser
ecological impact, and an “extreme” collection of scenarios of greater
ecological impact.

Clustering the scenarios based on these predicted range maps is a difficult
task due to the discrete and high dimensional nature of the range maps.
Although Principal Component Analysis (PCA) is a common tool to analyze data,
it has two significant drawbacks in this setting. One is that PCA implicitly
measures distance between range maps in Euclidean space, whereas we present a
flexible alternative, allowing for any similarity or distance between range
maps. Secondly, it is not clear how to incorporate information from multiple
species into a PCA-based approach.

Currently, a popular metric of change in species richness due to climate
change is climate velocity (Burrows et al., 2014; Jones and Cheung, 2015;
Molinos et al., 2016). However, these climate velocities rely heavily on
climate information while ignoring ecological data. Such ecological
information, such as species’ exposure to climate conditions not found in
their current niches, are important factors to predict the species future
ranges (Trisos et al., 2020). We incorporate such ecological data by
clustering based on predicted range maps based on models of historical ranges.

We will contrast this ecologically driven clustering with a climate driven
clustering to emphasise the importance of incorporating ecological
information. This climate driven clustering will use only predicted climate
features taken from CMIP6, such as annual mean temperature and annual mean
precipitation. By comparing the climate driven clustering and the ecological
one obtained from the range map similarities, we demonstrate the need for
incorporating ecological information. Clustering based only on climate
features is a similar but distinct task to the ecological based clustering;
the climate features are continuous whereas the ecological range maps are
binary. Previously these scenarios have been clustered using the climate
features by averaging these climate features over global regions Giorgi and
Francisco (2000); Cannon (2015). However, spatially averaging the climate
features this way loses significant information about the spatial variability
of the features.

We propose a flexible two step approach to cluster both predicted species
range maps and predicted climate features. The first step is to measure the
pairwise similarity between the prediction maps. The second step is to use the
pairwise similarities for spectral clustering, whose implementation is
discussed in Section 2. This two-step procedure is highly flexible, any
similarity measure can be used between range maps.

We propose to measure the similarity between range maps as the cosine
similarity of the range maps. This choice allows for considerable flexibility;
the modeller can weight absences and presences separately, and give certain
sets of cells higher importance. In addition the cosine similarity is
interpretable and comparable across species; it is always in the range
$[-1,1]$. This interpretability allows for a simple method to combine
information across species, as will be discussed in Section 2.2.

##  2 Materials and Methods

We first describe prior modelling work to predict the future range maps given
climate information. Then we discuss our methodology for clustering the
scenarios; first based on on predicted range maps, and secondly based on
predicted climate features. These clusterings reveal interpretable
relationships between the scenarios, and the contrast of these two clusterings
demonstrate the importance of incorporating ecological information.

###  2.1 Data: Predicted Range Map Modelling

A set of 9 GCMs are used from CMIP6 (Eyring et al., 2016; Stouffer et al.,
2017), which span four different levels of RCPS. Five climate features are
chosen from the set of 19 commonly used in WordCLIM (Fick and Hijmans, 2017),
which were chosen to minimize correlations. These five climate features are
annual mean temperature, temperature seasonality (standard deviation of
temperature), annual precipitation, precipitation seasonality (standard
deviation), and precipitation in the driest quarter of a year.

The five chosen climate features are used to fit a Poisson point process model
to explain the present day range maps (Merow et al., 2013; Elith et al.,
2011). The present day occurrences of 1101 mammals are obtained from Miller
(2020), whose range maps contain at 10 unique presence cells on a 10km grid.
This Poisson point process model was used to predict future spatial occurrence
based on these five predicted climate values. Binary maps are obtained from
these abundance maps by thresholding based on the $5$th percentile of
predicted values at training presences. This approach is used to make
predictions on all 1101 mammals using the 34 different sets of predicted
climate features. We make use of these predicted range maps to measure the
similarities between the 34 different climate scenarios.

###  2.2 Range Map Clustering Methodology

We present a novel methodology to cluster scenarios based on the binary
presence maps such as those shown in Figure 1(b). We now discuss the procedure
for achieving this scenario clustering. This clustering illuminates the
important similarities between scenarios, for instance if GCM or RCP is the
main differentiating feature among scenarios, and additionally gives insight
into the variation among scenarios.

####  2.2.1 Range Map Notation

For notational simplicity we focus the presentation of the methodology on a
single species. Furthermore, we suppose our maps are on a $n_{r}$ by $n_{c}$
grid represented as $(r,c):r=1,\ldots,n_{r},c=1,\ldots,n_{c}$. We denote
$B_{s}$ as the binary presence map of this species according to scenario $s$,
to be more precise $B_{s}(r,c)=1$ if the cell at $(r,c)$ represents a
presence, and $B_{s}(r,c)=0$ otherwise. We recognize it is important to
consider how these range maps differ from the present day as will be discussed
further. For this reason we denote by $P$ the present day map. Similarly as
for $B(r,c)$, we write $P(r,c)=1$ if the cell at $(r,c)$ is presently
occupied, and $0$ otherwise. We let the binary map $A$ (same size as $P$ and
each $B_{s}$) denote the background set, where valid absences may occur. For
instance, the set $A$ can take the value 1 only when there is land, or
alternatively only on the same continent or regional areas as $P$. We have
chosen to take $A$ as the union of the original map and all scenarios, thus
$A(r,c)=1$ if $P(r,c)=1$, or there is a scenario $s$ such that $B_{s}(r,c)=1$.

####  2.2.2 Cosine Similarity

In order to interpret the differences between range map predictions, a
similarity or distance measure is called for between pairs of range maps. The
quantification of the similarity or distance between range maps allows us to
cluster the scenarios, which is a desirably interpretable result.

There exists a plethora of alternative measures to quantify the similarity or
distance between these presence maps (Visser and De Nijs, 2006; Wilson, 2011;
Hagen, 2002; Hill et al., 2013; Gritti et al., 2013). The Hellinger distance
and Kullback-Leibler divergence rely on a probabilistic interpretation, and
thus it is difficult to incorporate absences into these measures (Wilson,
2011). The Kappa statistic (Hagen, 2002) can be used to measure similarities
between categorical maps, however it is not as clear how to weight important
cells, such as novel absences. In addition, the Kappa statistic does not
directly incorporate the relative frequencies of absences and presences. For
instance if two maps both predict all presences except a single absence cell,
these two maps will have a negative Kappa statistic if their one absence cell
differs. However, for our purposes we would like to consider such a pair of
maps very similar. The Wasserstein distance (or Earth Mover’s Distance) (Peyré
et al., 2019; Kranstauber et al., 2017) is an attractive alternative which
captures the idea of movement of a species. However, the Wasserstein distance
does not model the disappearance of regions (as opposed to the shift/movement
of regions), and computing the distance via optimal discrete transport (a
linear program of quadratic size in the number of present cells) proved too
slow or even infeasible (no solution found to the linear program) in all but
the smallest maps (the maps with the smallest number of presence cells).

We chose the cosine similarity function, which has several attractive
features. It is flexible; the modeller can weight absences and presences
separately, and the modeller has the flexibility to give certain sets of cells
higher importance. The cosine similarity is interpretable and comparable
across species as we describe below.

We compute the pairwise similarity between range maps as the cosine similarity
of their weighted range maps. The motivation and details for weighting the
range maps is discussed in Section 2.2.3. In short, we construct a weighted
map $W_{s}$ based on the binary range map $B_{s}$ and present day map, $P$.
For a pair of scenarios $s,s^{\prime}$ with corresponding binary range maps
$B_{s},B_{s^{\prime}}$ and present day map $P$, we measure their similarity as
the cosine similarity of their weighted maps,
$\text{CS}(W_{s},W_{s^{\prime}})$. We suppose that a map $G$ on an $n_{r}$ by
$n_{c}$ grid can be interpreted as a vector with length $n_{r}\times n_{c}$.
Then the cosine similarity between two vectors $W_{s},W_{s^{\prime}}$ can be
written as

$\text{CS}(W_{s},W_{s^{\prime}})=\frac{W_{s}\cdot
W_{s^{\prime}}}{\|W_{s}\|_{2}\|W_{s^{\prime}}\|_{2}}=\frac{\sum\limits_{r=1}^{n_{r}}\sum\limits_{c=1}^{n_{c}}W_{s}(r,c)W_{s^{\prime}}(r,c)}{\|W_{s}\|_{2}\|W_{s^{\prime}}\|_{2}},\hskip
5.69046pt\text{ with }\hskip
5.69046pt\|W_{s}\|_{2}=\Big{(}\sum\limits_{r=1}^{n_{r}}\sum\limits_{c=1}^{n_{c}}W_{s}(r,c)^{2}\Big{)}^{1/2}.$
(1)

####  2.2.3 Weighted Map

Weighting the cells of the binary maps is required to quantitatively measure
the similarity between range maps. For each scenario $s$, we denote the
corresponding weighted map $W_{s}$. The value of $W_{s}(r,c)$ is 0 if
$A(r,c)=0$, but when $A(r,c)=1$ we choose one of four values for $W_{s}(r,c)$
depending on whether the cell represents a presence/absence for the present
day, and a presence/absence for the scenario. These four choices are shown in
Table 1.

| $P(r,c)=1$ | $P(r,c)=0$
---|---|---
$B(r,c)=1$ | $p_{\text{keep}}{}$ | $p_{\text{new}}{}$
$B(r,c)=0$ | $a_{\text{new}}{}$ | $a_{\text{keep}}{}$
Table 1:  Cell Weighting Values for $W(r,c)$ given $A(r,c)=1$.

The choice of weightings in Table 1 is dependent on the use of cosine
similarity. For instance using cosine similarity but representing presences
with the value 1 and absences with 0 is inappropriate because any two maps
will have positive similarity, and this choice will also have the property
that absences and presences are counted differently in the denominator of the
cosine similarity.

These first four cases (reading row by row) represent unchanged (kept)
presences, new presences, new absences, and unchanged absences respectively,
where “keep” and “new” are with respect to the present day distribution, $P$.
For example, a cell corresponding to $p_{\text{keep}}{}$ is a cell that is
present in $P$, and is kept present according to scenario $B_{s}$. Presences
are given positive weights, absences negative weights, and we choose
$|a_{\text{new}}{}|>|a_{\text{keep}}{}|$, to emphasize those cells whose
ecological suitability for this species is vanishing. These “lost” cells
corresponding to novel absences are particularly important; they represent
cells where the climate is changing so drastically that the species cannot
continue to live there. In addition these cells should be weighted higher
because we are more confident about their prediction; by definition they exist
within the training data’s presences. By contrast the cells corresponding to
$p_{\text{new}}{}$ represent regions where the species are predicted to move
to, however such predictions are more uncertain as the movement of a species
is complex and not directly taken into account by the Poisson point process
model. Thus we also make the choice $|p_{\text{keep}}{}|>|p_{\text{new}}{}|$.
We make the choice
$|a_{\text{new}}{}|=|p_{\text{keep}}{}|=1,|a_{\text{keep}}{}|=|p_{\text{new}}{}|=0.5$
to represent the fact we are more confident about all cells in the region of
present day occurrences (those cells corresponding to $P(r,c)=1$). However, we
emphasize the flexibility of our model, alternative choices can be made
depending on the modeler’s goals. A visualization of this weighting scheme is
shown in Figure 1.

(a) Present Day Range Map

(b) Scenario Range Map

(c) Overlap of Present Day and Scenario

Figure 1:  Illustration for weighted presences and absences. Presences are
shown in blue, absences in red/pink. The overlap with the present day in shown
in darker regions which represent more significant cells, which are dark blue
kept presences, and dark red new absences, which represent “lost” cells. A
strength of our proposed methodology is the flexibility to weight these cells
differently

The computation and resultant cosine similarity using these weighted maps are
highly interpretable; when two scenarios agree on the presence or absence of a
cell, this cell has a positive contribution to the cosine similarity, whereas
the cell has a negative contribution when the two scenarios disagree. The
resulting similarity is always in the range $[-1,1]$ (for any choice of
weightings), and takes the value $1$ when the two maps are identical and $-1$
if the two maps completely disagree on presences, as long as
$|a_{\text{new}}{}|=|p_{\text{keep}}{}|,|a_{\text{keep}}{}|=|p_{\text{new}}{}|$,
which was argued for previously.

####  2.2.4 Spectral Clustering

For all species we compute the pairwise scenario similarity matrix by
computing the cosine similarity, eq. (1) on each pair of scenarios. That is
for each species $m$, we construct the $n_{s}$-by-$n_{s}$ matrix
$S^{m}(s,s^{\prime})=CS(W^{m}_{s},W^{m}_{s^{\prime}})$, where $n_{s}$ is the
number of scenarios and $W^{m}_{s}$ is the weighting map for species $m$ on
scenario $s$. Spectral clustering is well suited to cluster the scenarios
based on this similarity measure (Von Luxburg, 2007).

The properties of spectral clustering are understood from a graph theory
perspective. The similarity matrix $S^{m}$ can be thought of as an undirected
graph whose nodes are the scenarios and edge weights between a pair of
scenarios $s$ and $s^{\prime}$ given by $S^{m}(s,s^{\prime})$. Spectral
clustering has best performance on sparse graphs, thus for the dense
similarity matrix $S^{m}$, the first step is to sparsify it by means of taking
the $k$-nearest neighbor graph, that is retaining an edge from $s$ to
$s^{\prime}$ only if $s^{\prime}$ is within the top $k$ neighbors of $s$ (i.e.
it is within the top $k$ nodes of maximal similarity to $s$). However this
would lead to a directed graph as this definition of nearest neighbor is not
symmetric. Thus we retain the undirected edge from $s$ to $s^{\prime}$ if
either $s^{\prime}$ is within the top $k$ neighbors of $s$, or vice-versa. The
retained edges are still weighted by the similarity of their endpoints. We
denote by $E$ this matrix of retained weights.

The main tool of spectral clustering is the graph Laplacian $L=D-E$, where $D$
is a diagonal matrix of node degree, $D_{ii}=\sum_{j=1}^{n_{s}}E(i,j)$. We use
the random-walk normalized graph Laplacian, $L_{rw}=D^{-1}L$ as suggested in
Von Luxburg (2007). Both $L$ and $L_{rw}$ have several nice properties, they
are positive semi-definite and the multiplicity of their zero eigenvalue is
the number of connected components of the graph. For most real-world graphs
including the sparsified cosine similarity matrix $E$, the graph is fully
connected and thus the number of connected components of this graph is one.
When this is the case the eigenvectors corresponding to the smallest non-zero
eigenvalues can be used as an embedding. This can be thought of from a
perturbation perspective, if the graph had true clusters of connected
components than these eigenvectors would have the same span as vectors
representing clustering membership indicators. However in real world graphs
there are a few “noisy” edges between clusters, thus these first few
eigenvectors instead are near piecewise on those indicators.

Drawing from this insight, the space associated with the eigenvectors of
$L_{rw}$, which we denote as columns of a matrix $U$, is used as “spectral
embedding”. We use the second and third eigenvectors of $U$ (corresponding to
the second and third smallest eigenvalue) as the spectral embedding , that is
scenario $s$ is represented by $(U_{s,2},U_{s,3})$. We choose a simple
clustering algorithm, single-linkage clustering, to cluster in this embedded
space as it performs well according to the Davies-Bouldin criterion, a common
clustering criterion for how well separated clusters are (Davies and Bouldin,
1979).

This clustering can be performed for a single species using only $S^{m}$, for
example this is performed for four different species, shown in Figure 6. These
results will be discussed in the next section. One way to combine information
across animals is to average the similarity matrices across all $n_{m}=1101$
animals, as in $S=(1/n_{m})\sum_{a=1}^{n_{m}}S^{m}$. Performing spectral
clustering using this averaged similarity matrix, $S$, is shown in Figure
5(b).

###  2.3 Global Spatial Diversity Loss: Data Summary Visualization

One can get an overall sense of the changes predicted in the range maps from
Figure 2, which shows how the diversity of mammals is spread spatially
throughout the globe. We see from both Figures 2(c) and 2(d) that significant
diversity losses are predicted around the equator in South America and Africa,
whereas there is some diversity increases further up north, consistent with
previous findings (Chen et al., 2011).

(a) Present Day Mammal Variety

(b) Predicted Mammal Variety

(c) Net Change

(d) Fraction Change

Figure 2:  Visualization of the mammal richness in the dataset over space.
White cells correspond to locations with no recorded presences. Our findings
are consistent with Chen et al. (2011) which finds that species are moving
poleward and towards higher elevations, there is loss around the equator and
some increase in diversity towards the northern pole.

####  2.3.1 Frequency Analysis

One way to qualitatively measure the performance of the clustering is to look
at the spatial overlap of presences for each scenario. That is we define the
frequency map $F^{m}$ based on the scenario maps:
$F^{m}(r,c):=\sum_{s=1}^{n_{s}}B^{m}_{s}(r,c)$ For example the overlap of all
34 scenarios for the African cheetah is shown in Figure 8(a).

####  2.3.2 Principal Component Analysis

Principal component analysis is a great tool to understand the main directions
of variation in data. We use PCA to visualize the correlations of presence
cells. This can be performed for a single species by considering each scenario
as an observation, with $n_{r}\times n_{c}$ binary features corresponding to
the vectorized (flattened) binary range map.

####  2.3.3 Climate Based Scenario Clustering

We discuss the process of clustering using only the five climate features that
were used to predict the range maps. This clustering when contrasted with the
ecologically driven clustering demonstrates the importance of incorporating
ecological information. Each of these five globally distributed features is
predicted across the 34 scenarios. In order to directly compare to our
ecologically based clustering, we utilize the climate features to perform a
clustering in a similar fashion to the ecologically based clustering. However,
for continuous data, the cosine similarity is not appropriate, as two maps
that are scaled versions of each other would be considered very similar to
each other, which is not desirable. For instance if one map predicted two
degrees warmer everywhere than another, the cosine similarity between these
two maps would very high, which is not desirable as these two maps represent
significantly different predictions. Instead we use the $L_{2}$ distance
between maps, which will effectively use both the difference between the means
of maps, and differences in the spatial variation. To incorporate all five
climate features, each feature is normalized before applying the $L_{2}$
distance. We denote $T_{s}^{f}(r,c)$ the value of feature $f$ according to
scenario $s$ at location $(r,c)$. The feature scaled $L_{2}$ distance between
a pair of scenarios $s$ and $s^{\prime}$ using all five features is given by:

$H(s,s^{\prime})=\sum_{f=1}^{5}\sum_{r=1}^{n_{r}}\sum_{c=1}^{n_{c}}[(T_{s}^{f}(r,c)-T_{s^{\prime}}^{f}(r,c))/\sigma_{f}]^{2}.$

Where $\sigma_{f}$ is the standard deviation of feature $f$ measured across
all locations and scenarios, that is we calculate the mean of each feature
across all locations and scenarios,

$\sigma_{f}^{2}=\sum_{s=1}^{n_{s}}\sum_{r=1}^{n_{r}}\sum_{c=1}^{n_{c}}[T_{s}^{f}(r,c)-\mu_{f}]^{2}\hskip
5.69046pt\text{ with }\hskip 5.69046pt\mu_{f}=(n_{s}\cdot n_{r}\cdot
n_{c})^{-1}\sum_{s=1}^{n_{s}}\sum_{r=1}^{n_{r}}\sum_{c=1}^{n_{c}}T_{s}^{f}(r,c).$
(2)

Spectral clustering requires a similarity matrix (instead of a dissimilarity
or distance matrix), thus we transform the pairwise $L_{2}$ distances into
similarities with any monotonically negative function, such as $x\rightarrow
1/x$. An important step in spectral clustering is to sparsify the similarity
matrix by keeping only the largest similarities (discussed above), and so any
monotonically negative function will produce the same spectral embedding and
clusters Von Luxburg (2007). Thus we use the similarity matrix
$L(s,s^{\prime}):=1/H(s,s^{\prime})$ as the similarity matrix for spectral
clustering, which is shown in Figure 5(b).

(a) Annual Mean Temperature

(b) Temperature Seasonality

(c) Annual Precipitation

(d) Precipitation Seasonality

Figure 3:  Spectral clustering of scenarios using individual climate features
(one of the five is not shown due to space). The annual temperature clustering
is most similar to the ecological clustering of Figure 5(b), clustering mainly
by RCP. The other features mainly cluster mainly by GCM. This is why the
climate driven clustering of Figure 5(a) is driven mainly by GCM, most of the
features are. The ecological clustering is important to discern which of these
features are most ecologically relevant, these plots show that annual
temperature contains the most ecologically relevant differences between
climate models.

###  2.4 Variation of Scenarios

We summarize some overall patterns of the projected range maps by counting
both presences that differ from the present day range map, and absences that
differ from the present day range maps (cells corresponding to
$a_{\text{new}}{}$ and $p_{\text{new}}{}$ respectively). The fraction of each
of these cell types is averaged over species for each scenario. The resulting
average of fraction of $a_{\text{new}}{}$ and $p_{\text{new}}{}$ type cells is
shown in Figure 4, which shows that these novel absences and novel presences
are correlated.

##  3 Results and Discussion

In this section we discuss the scenario clusterings, with emphasis on the
difference between the climate based and ecological based clustering. We then
demonstrate that we have detected meaningful clusters, with range maps within
clusters similar to each other, but different than range maps in other
clusters.

###  3.1 Range Shifts

By comparing the fraction of cells corresponding to $a_{\text{new}}{}$ and
$p_{\text{new}}{}$ we get a sense of the changes of the range maps compared to
the present day range map. This is shown in Figure 4, which demonstrates these
range maps tend to predict range shifts, as opposed to range expansions or
contractions. We conclude this because the number of new presences and new
absences grow together, which would occur as range shifts, instead of say
absences growing as presences shrink, which would be indicative of an overall
range decrease.

On average there are significantly more novel absences (average 38% (28% s.d.)
of current range) than novel presences (average 25% (26% s.d.)), which implies
that the species ranges are shrinking (average 13% net loss, (31% s.d.)) due
to the changing climate.

Figure 4:  Shows that predicted range maps tend to be shifted, that is lost
and new cells grow together, with losses being larger than gains, as the
scenarios are above the 45 degree dashed black line.

###  3.2 Clustering Plots

The spectral clustering results using the species specific similarity matrix
is shown in Figure 6 for four species, illustrating the main types of patterns
observed. We see an interesting mix of RCP and GCM dependence. It appears that
RCP is a major driving factor that separates these clusters; in most of these
clusterings the far left (“optimistic”) cluster contains mainly scenarios with
low RCP, and the far right (“extreme”) cluster only scenarios with high RCP.
However, the clustering using the slender treeshew (_Tupaia gracilis_) is
driven mainly by GCM. This clustering mainly by GCM was found in many species
(30%), although the most common trend is RCP dependence (70% of species). The
variability among animals is further evidence supporting the importance of the
climate ecology relationship: the most important difference between climate
models (GCM or RCP) varies depending on the individual species’ response to
the climate.

The spectral cluster results from using the similarity matrix averaged over
all 1101 mammals is shown in Figure 5(b). We see in the bottom cluster of
Figure 5(b) an interesting mix of varying RCP and GCM. This suggests that RCP
alone does not account for the variation, the GCM is also important. However,
there are still variabilities among the GCMs. For example, the “ca” climate
model appears twice in the far right cluster.

(a) Climate Feature Based

(b) Predicted Range Map Based

Figure 5:  In the climate based clustering shown on the left, the clusters are
mainly determined by GCM. This clustering is significantly different from the
ecological based one shown on the right, which demonstrates the importance of
incorporating ecological information. B: Scenario Embedding and Clustering.
The clusters are mainly driven by RCP. However, there is still some
relationships among the GCM, for instance the red ’ca’ GCM predicts more
extreme outcomes, and at each level of RCP the green black and yellow
(cc,ce,ip) climate models are together. This clustering is significantly
different from the climate based one (Panel A), which demonstrates the
importance of incorporating ecological information.

Instead of averaging over all species, we also performed clustering for only
the species most at risk, defined by those species whose fraction of area lost
is among the highest $10\%$. This loss in area can be used to approximate a
loss in population abundance using the techniques in He (2012); Che-Castaldo
and Neel (2016). Performing spectral clustering by using the average of
similarity matrices of this subset of animals at risk in shown in Figure 7.
This clustering puts a stronger emphasis on RCP, that is, the animals most at
risk are more sensitive to RCP.

(a) _Acinonyx_ jubatus (cheetah)

(b) _Zapus princeps_

(c) _Procyon lotor_

(d) _Tupaia gracilis_

Figure 6:  Spectral clustering of scenarios using individual species. We see
mainly grouping by RCP for the first three, but a starkly differently mainly
GCM driven clustering for the treeshew. This variety was found among animals,
most species driven clusterings are strongly connected with RCP, but some are
more connected with GCM. These differences between species further
demonstrates that ecological information is important to interpret the
differences between scenarios. Figure 7:  Scenario Embedding and Clustering
using only animals in riskiest quantiles. Clustering based only these species
most at risks puts an even higher emphasis on RCP than the clustering in
Figure 5(b). This further demonstrates the importance of accounting for
ecological information, different subsets of ecological populations emphasis
RCP even more strongly.

###  3.3 Cluster Quality Analysis

One way to qualitatively measure the performance of the clustering is to look
at the spatial overlap of presences for each scenario. For example the overlap
of all 34 scenarios is shown in Figure 8(a). We see that the cluster
associated with the smallest RCP scenarios (“optimistic” scenarios) accounts
for most of the presences in the discrepant regions, whereas the cluster
associated with the highest RCP (“extreme” scenarios) accounts for many of the
absences. This demonstrates that we have discovered meaningful clusters; there
is agreement within clusters but disagreements across clusters.

We visualize the correlations of presence cells over scenarios using principal
component analysis. This can be performed for a single species by considering
each scenario as an observation, with $n_{r}\times n_{c}$ binary features
corresponding to the vectorized (flattened) binary range map. Performing PCA
in this way is shown for the cheetah in Figure 8(b). The mix of positive and
negative coordinates illustrate the fact that the most extreme scenarios do
not only “lose” certain cells, but also predict novel presences.

(a) All Scenarios

(b) First PC

(c) Optimistic Scenarios

(d) Exteme Scenarios

Figure 8:  Sum of presences over scenarios. A meaningful clustering should
have strong similarities within cluster, and differences across cluster. The
circled regions showed that indeed we have discovered meaningful clusters;
there is agreement within cluster in these regions, but differences across
clusters.

###  3.4 Individual Climate Feature Clustering

We see that the climate based clustering is significantly different than the
ecological one, the clustering is driven mainly by GCM instead of by RCP in
the ecological based clustering. This is evidence for the importance of
considering the specific climate niche occupied by a species in relation to
how those conditions are projected to change. Although these same five
features are used to predict the ecological range maps, these predicted range
maps paint a different picture of the scenario clustering because of how the
animals are influenced by the climate features. In fact, we can get a sense of
feature importance by clustering based on individual climate features, shown
in Figure 3. We see that clustering using only the annual temperature creates
the most similar clustering to the ecologically driven one, suggesting that
annual temperature is the most important feature of these five.

##  4 Conclusion

We have proposed a novel methodology to cluster scenarios based on ecological
range maps. The presented approach is interpretable, flexible, and fast. We
have demonstrated different patterns of clustering depending on which ]subsets
of species are included. For instance, the animals most at risk group many of
the higher RCP scenarios together. The differences between the climate based
clustering and ecological based clustering highlights the importance of
considering ecological response; the interaction of climate and ecology is
essential to understand the ecologocially most important differences between
future scenario predictions. An interesting direction to explore further is to
uncover subsets of animals that respond differently than others. For instance
it may be the case that rodents tend to fare worse under the “bc” climate
model, compared to other mammals. A similar area of future research is to
determine why some species like the slender treeshew cluster more by GCM
instead of the more pattern of RCP. Are these species less sensitive to
temperature changes? Another extension is to consider how to combine
information across animals in a more holistic manner. For instance, Dong et
al. (2013) presents a methodology to cluster according to many graphs, which
could be applied in our use case to the scenario graphs from each species.

## References

  * Parry et al. (2007) Martin L Parry, Osvaldo Canziani, Jean Palutikof, Paul Van der Linden, and Clair Hanson.  _Climate change 2007-impacts, adaptation and vulnerability: Working group II contribution to the fourth assessment report of the IPCC_ , volume 4.  Cambridge University Press, 2007.
  * Hannah et al. (2013) Lee Hannah, Makihiko Ikegami, David G Hole, Changwan Seo, Stuart HM Butchart, A Townsend Peterson, and Patrick R Roehrdanz.  Global climate change adaptation priorities for biodiversity and food security.  _PLoS one_ , 8(8):e72590, 2013.
  * Hannah et al. (2020) Lee Hannah, Patrick R Roehrdanz, Pablo A Marquet, Brian J Enquist, Guy Midgley, Wendy Foden, Jon C Lovett, Richard T Corlett, Derek Corcoran, Stuart HM Butchart, et al.  30% land conservation and climate action reduces tropical extinction risk by more than 50%.  _Ecography_ , 2020.
  * Burrows et al. (2014) Michael T Burrows, David S Schoeman, Anthony J Richardson, Jorge Garcia Molinos, Ary Hoffmann, Lauren B Buckley, Pippa J Moore, Christopher J Brown, John F Bruno, Carlos M Duarte, et al.  Geographical limits to species-range shifts are suggested by climate velocity.  _Nature_ , 507(7493):492–495, 2014.
  * Jones and Cheung (2015) Miranda C Jones and William WL Cheung.  Multi-model ensemble projections of climate change effects on global marine biodiversity.  _ICES Journal of Marine Science_ , 72(3):741–752, 2015.
  * Molinos et al. (2016) Jorge García Molinos, Benjamin S Halpern, David S Schoeman, Christopher J Brown, Wolfgang Kiessling, Pippa J Moore, John M Pandolfi, Elvira S Poloczanska, Anthony J Richardson, and Michael T Burrows.  Climate velocity and the future global redistribution of marine biodiversity.  _Nature Climate Change_ , 6(1):83–88, 2016.
  * Eyring et al. (2016) Veronika Eyring, Sandrine Bony, Gerald A Meehl, Catherine A Senior, Bjorn Stevens, Ronald J Stouffer, and Karl E Taylor.  Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization.  _Geoscientific Model Development_ , 9(5):1937–1958, 2016.
  * Trisos et al. (2020) Christopher H Trisos, Cory Merow, and Alex L Pigot.  The projected timing of abrupt ecological disruption from climate change.  _Nature_ , 580(7804):496–501, 2020.
  * Giorgi and Francisco (2000) Filippo Giorgi and Raquel Francisco.  Uncertainties in regional climate change prediction: a regional analysis of ensemble simulations with the hadcm2 coupled aogcm.  _Climate Dynamics_ , 16(2-3):169–182, 2000.
  * Cannon (2015) Alex J Cannon.  Selecting gcm scenarios that span the range of changes in a multimodel ensemble: application to cmip5 climate extremes indices.  _Journal of Climate_ , 28(3):1260–1267, 2015\.
  * Stouffer et al. (2017) Ronald J Stouffer, Veronika Eyring, Gerald A Meehl, Sandrine Bony, Cath Senior, Bjorn Stevens, and KE Taylor.  CMIP5 scientific gaps and recommendations for CMIP6.  _Bulletin of the American Meteorological Society_ , 98(1):95–105, 2017.
  * Fick and Hijmans (2017) Stephen E Fick and Robert J Hijmans.  Worldclim 2: new 1-km spatial resolution climate surfaces for global land areas.  _International journal of climatology_ , 37(12):4302–4315, 2017.
  * Merow et al. (2013) Cory Merow, Matthew J Smith, and John A Silander Jr.  A practical guide to maxent for modeling species’ distributions: what it does, and why inputs and settings matter.  _Ecography_ , 36(10):1058–1069, 2013.
  * Elith et al. (2011) Jane Elith, Steven J Phillips, Trevor Hastie, Miroslav Dudík, Yung En Chee, and Colin J Yates.  A statistical explanation of maxent for ecologists.  _Diversity and distributions_ , 17(1):43–57, 2011\.
  * Miller (2020) Joe Miller.  Gbif home page, 2020.  URL https://www.gbif.org.
  * Visser and De Nijs (2006) Hans Visser and T De Nijs.  The map comparison kit.  _Environmental Modelling & Software_, 21(3):346–358, 2006.
  * Wilson (2011) Peter D Wilson.  Distance-based methods for the analysis of maps produced by species distribution models.  _Methods in Ecology and Evolution_ , 2(6):623–633, 2011.
  * Hagen (2002) Alex Hagen.  Multi-method assessment of map similarity.  In _Proceedings of the 5th AGILE Conference on Geographic Information Science_ , pages 171–182. Universitat de les Illes Balears Palma, Spain, 2002.
  * Hill et al. (2013) Mark O Hill, Colin A Harrower, and Christopher D Preston.  Spherical k-means clustering is good for interpreting multivariate species occurrence data.  _Methods in Ecology and Evolution_ , 4(6):542–551, 2013.
  * Gritti et al. (2013) Emmanuel S Gritti, Anne Duputie, Francois Massol, and Isabelle Chuine.  Estimating consensus and associated uncertainty between inherently different species distribution models.  _Methods in Ecology and Evolution_ , 4(5):442–452, 2013.
  * Peyré et al. (2019) Gabriel Peyré, Marco Cuturi, et al.  Computational optimal transport: With applications to data science.  _Foundations and Trends® in Machine Learning_ , 11(5-6):355–607, 2019.
  * Kranstauber et al. (2017) Bart Kranstauber, Marco Smolla, and Kamran Safi.  Similarity in spatial utilization distributions measured by the earth mover’s distance.  _Methods in Ecology and Evolution_ , 8(2):155–160, 2017.
  * Von Luxburg (2007) Ulrike Von Luxburg.  A tutorial on spectral clustering.  _Statistics and computing_ , 17(4):395–416, 2007\.
  * Davies and Bouldin (1979) David L Davies and Donald W Bouldin.  A cluster separation measure.  _IEEE transactions on pattern analysis and machine intelligence_ , 2(2):224–227, 1979.
  * Chen et al. (2011) I-Ching Chen, Jane K Hill, Ralf Ohlemüller, David B Roy, and Chris D Thomas.  Rapid range shifts of species associated with high levels of climate warming.  _Science_ , 333(6045):1024–1026, 2011.
  * He (2012) Fangliang He.  Area-based assessment of extinction risk.  _Ecology_ , 93(5):974–980, 2012.
  * Che-Castaldo and Neel (2016) Judy P Che-Castaldo and Maile C Neel.  Species-level persistence probabilities for recovery and conservation status assessment.  _Conservation Biology_ , 30(6):1297–1306, 2016\.
  * Dong et al. (2013) Xiaowen Dong, Pascal Frossard, Pierre Vandergheynst, and Nikolai Nefedov.  Clustering on multi-layer graphs via subspace analysis on Grassmann manifolds.  _IEEE Transactions on signal processing_ , 62(4):905–918, 2013.