See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/310702405

ecospat: An R package to support spatial analyses and modeling of species

niches and distributions

Article  in  Ecography · February 2017

DOI: 10.1111/ecog.02671

CITATIONS READS

646 5,924

15 authors, including:

Valeria Di Cola Olivier Broennimann

Swiss Institute of Bioinformatics University of Lausanne
63 PUBLICATIONS   979 CITATIONS    151 PUBLICATIONS   11,687 CITATIONS   

SEE PROFILE SEE PROFILE

Blaise Petitpierre Robin Engler

44 PUBLICATIONS   4,615 CITATIONS    37 PUBLICATIONS   8,108 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Quantifying the temporal evolution of the protection service of forests against rockfall in the face of climate change View project

Towards a national restoration plan considering connectivity and vulnerability to climate change (connect2restore) View project

All content following this page was uploaded by Rubén G. Mateo on 24 January 2018.

The user has requested enhancement of the downloaded file.

 Ecography 40: 001–014, 2017
doi: 10.1111/ecog.02671
© 2016 The Authors. Ecography © 2016 Nordic Society Oikos
Subject Editor: Thiago Rangel. Editor-in-Chief: Miguel Araújo. Accepted 17 November 2016

ecospat: an R package to support spatial analyses and modeling of

species niches and distributions

Valeria Di Cola, Olivier Broennimann, Blaise Petitpierre, Frank T. Breiner, Manuela D’Amen,
Christophe Randin, Robin Engler, Julien Pottier, Dorothea Pio, Anne Dubuis, Loic Pellissier,
Rubén G. Mateo, Wim Hordijk, Nicolas Salamin and Antoine Guisan
V. Di Cola ([email protected]), O. Broennimann, B. Petitpierre, F. T. Breiner (http://orcid.org/0000-0003-4465-1684), M. D’Amen, C.
Randin, D. Pio, A. Dubuis, R. G. Mateo, W. Hordijk, N. Salamin and A. Guisan, Dept of Ecology and Evolution (DEE) – Univ. of Lausanne,
Lausanne, Switzerland. CR also at: Centre de Recherches sur les Ecosystèmes d’Altitude CREA, Observatoire du Mont-Blanc, Chamonix, France.
DP also at: Fauna and Flora International, London, UK. WH also at: Konrad Lorenz Inst. for Evolution and Cognition Research, Klosterneuburg,
Austria. AG also at: Inst. of Earth Surface Dynamics (IDYST), Univ. of Lausanne, Lausanne, Switzerland. – R. Engler, Vital-IT group, SIB Swiss
Inst. of Bioinformatics, Génopode, Lausanne, Switzerland. RE and NS also at: SIB Swiss Inst. of Bioinformatics, Génopode, Lausanne, Switzerland.
– J. Pottier, INRA, UR874-Grassland Ecosystem Research Unit, Clermont-Ferrand, France. – L. Pellissier, Landscape Ecology, Inst. of Terrestrial
Ecosystems, Zürich, Switzerland. LP, FTB and CR also at: Swiss Federal Research Inst. WSL, Birmensdorf, Switzerland.

The aim of the ecospat package is to make available novel tools and methods to support spatial analyses and modeling
of species niches and distributions in a coherent workflow. The package is written in the R language (R Development Core
Team) and contains several features, unique in their implementation, that are complementary to other existing R packages.
Pre-modeling analyses include species niche quantifications and comparisons between distinct ranges or time periods,
measures of phylogenetic diversity, and other data exploration functionalities (e.g. extrapolation detection, ExDet). Core
modeling brings together the new approach of ensemble of small models (ESM) and various implementations of the
spatially-explicit modeling of species assemblages (SESAM) framework. Post-modeling analyses include evaluation of
species predictions based on presence-only data (Boyce index) and of community predictions, phylogenetic diversity and
environmentally-constrained species co-occurrences analyses. The ecospat package also provides some functions to
supplement the ‘biomod2’ package (e.g. data preparation, permutation tests and cross-validation of model predictive
power). With this novel package, we intend to stimulate the use of comprehensive approaches in spatial modelling of
species and community distributions.

Background by the increasing availability of numerical geodata and the

tremendous progress in computer and statistical sciences
Spatial ecology lies at the intersection of geography and (Dawson et al. 2013). As a result, many new methods have
ecology and examines the role of space in ecological pro- been developed (Skidmore et al. 2011) and a still increas-
cesses from local to global scales. One of its overarching ing number of tools are becoming available to fit SDMs.
aims is to investigate what factors shape the distribution of Examples are ‘BioMapper’ (Hirzel et al. 2001), ‘Maxent’
species (Franklin 2009, Peterson et al. 2011) and, depart- (Phillips et al. 2006), ‘biomod2’ (Thuiller et al. 2016),
ing from these, to reconstruct communities and ecosystems ‘dismo’ (Hijmans et al. 2016), and ‘sdm’ (Naimi and Araújo
(Levin 1992, Dungan et al. 2002, Guisan and Rahbek 2016). SDMs are increasingly used to address fundamental
2011, Hastings et al. 2011). It is now recognized that four and applied questions in a wide range of domains, including
main factors shape the distribution of species in time and community modeling (D’Amen et al. 2015c), community
space: environmental suitability, colonisation history and phylogenetics (Ndiribe et al. 2014), co-occurrence analyses
ability, biotic interactions, and natural and anthropogenic (D’Amen et al. 2015a), biodiversity patterns (Mateo et al.
disturbances (Pulliam 2000, Lortie et al. 2004, Soberon 2016), anticipating and monitoring biological invasions
2007). In recent years, the use of niche-based, correlative (Guisan et al. 2014), rare species management (Breiner et al.
species distribution models (hereafter SDMs, but other 2015), or conservation phylogenetics (Pio et al. 2011), and
denominations are also used, such as ecological niche mod- the ecospat package offers functions to support analyses
els ENMs, habitat suitability models HSMs, or climatic in these different fields.
envelope models CEM; see Appendix S1 in Guisan et al. Depending on where they take place in the modeling
2013) has risen spectacularly (Guisan et al. 2013), fostered process, methods may be tentatively classified as pre-, core or

Early View (EV): 1-EV

post- modeling (or ‘processing’ in Naimi and Araújo 2016) In addition, two unique core modeling analyses in
analyses. Pre-modeling (or pre-processing) analyses are usu- ecospat are functions to run the very recently proposed
ally based on the observational data themselves, and include ensemble of small models (ESM) approach (Breiner et al.
niche quantification (Fitzpatrick et al. 2007, Medley 2010), 2015) on the one hand, and functions to implement the
co-occurrence (Gotelli 2000), or phylogenetic diversity SESAM framework for spatial predictions of species assem-
(Schweiger et al. 2008, Morlon et al. 2011) analyses. Core- blages (Guisan and Rahbek 2011, D’Amen et al. 2015a, b)
modeling analyses are typically based on the quantification of on the other hand. These functionalities represent new devel-
species–environment relationship, using different techniques opments currently only available in ecospat. The package
to calibrate models and use them to make spatial predictions also uniquely contains functions to calculate the Boyce index
(Guisan and Zimmermann 2000) both at the individual (continuous or based on predefined bins Hirzel et al. 2006)
species and community levels. Finally, post-modeling (or to evaluate species distribution predictions when only inde-
post-processing) analyses make use of model predictions, pendent presence data are available (Petitpierre et al. 2012).
for instance to incorporate a posteriori other filters of spe- A major post-modeling feature is the set of improved func-
cies distributions (Guisan and Rahbek 2011), such as biotic tions for environmentally-constrained co-occurrences analyses
interactions (Wisz et al. 2013, D’Amen et al. 2015b), or (based on Peres-Neto et al. 2001), where the environmen-
dispersal (Dullinger et al. 2012, Engler et al. 2012). Several tal constraint is implemented by using SDM predictions.
software and R packages exist to run core modeling analyses Another post-modeling functionality, nicely complementing
for individual species (see above), but fewer tools are avail- the SESAM framework (e.g. to set phylogenetic constraints on
able to support community-level modeling and the pre- and assemblages), is the function to predict phylogenetic diversity
post-modeling analyses presented here (but see e.g. ‘dismo’ in geographic space from the stacked predictions of species
Hijmans et al. 2016). distribution models (Pio et al. 2011, 2014), which can also be
used for simple pre-modeling analyses based on raw observa-
tions (Ndiribe et al. 2013, 2014, Pellissier et al. 2013a, b).
Aim of ecospat

The aim of the ecospat package is to provide methods Technical features of the package
and utilities for spatial and/or temporal predictions of spe-
cies distributions, communities’ properties and related anal- The ecospat package is written for use in the R language (R
yses (e.g. niche quantification, co-occurrence analyses). It Development Core Team). The latest ver. 2.1.1 was recently
contains a set of miscellaneous functions and methods to fit released in CRAN (The Comprehensive R Archive Network)
and evaluate SDMs and community models in a coherent can be installed and used under R ver. 3.3 or higher (for more
analytical framework. The package includes novel functions details see: < https://cran.r-project.org/package = ecospat > or
(see below) and improvements of existing ones to supple- use the link from < www.unil.ch/ecospat/home/menuguid/
ment those already available in other packages. The inter- ecospat-resources/tools.html >). The ecospat package
est of the package is that it possesses functions not only to depends on a series of other R packages: ‘ade4’ (Dray and
run core modeling analyses at the individual species level, Dufour 2007), ‘sp’ (Qiao et al. 2016), ‘raster’ (Hijmans
but also at the community level, and it additionally provides 2015), ‘biomod2’ (Thuiller et al. 2009), ‘gbm’ (Zhang 2013),
functions using raw observations and predictions to run pre- ‘MigClim’ (Engler et al. 2012) and ‘ape’ (Paradis et al. 2004).
and post-modeling (or processing) analyses respectively. For details on the package dependencies, see Supplementary
material Appendix 1, Fig. A1 and Table A1.
Most functions in the ecospat package can be used
to perform stand-alone analyses (e.g. niche or co-occurrence
Novelty analyses), some functions are specifically designed to be used
The ecospat package includes many novel features. in supplement to other species distribution modeling pack-
A major pre-modeling feature is the set of functions for ages (e.g. ‘biomod2’, ‘dismo’, ‘migclim’; see technical features
niche quantification and tests of observed overlap in con- above), but not to replace any functions in these packages. To
trast to null distributions that were used in Broennimann avoid potential conflicts with other R functions, all the func-
et al. (2012), Petitpierre et al. (2012) and many other stud- tions in the package start with the prefix ‘ecospat.’ (see
ies to quantify climatic niche shifts between the native and Table 1 for the main functions and their descriptions). The
invaded ranges of alien species (Guisan et al. 2014), to anal- package ecospat is available from < http://cran.r-project.
yse niche overlap between different species (Table 1) and to org/web/packages/ecospat/ >, together with all example data
understand temporal dynamics in species niches and distri- sets referred to in the text.
butions (Nogués-Bravo 2009, Maguire et al. 2015). We are It is being constantly updated and we encourage inter-
not aware of any other package (see below, section on niche ested users to look for the latest package version in CRAN or
analyses) – in or outside R – allowing to run the same niche GitHub (< https://github.com/cran/ecospat >) and to report
analyses, including the possibility to measure niche unfilling, any suggestion for further improvement.
stability and expansion (see the COUE framework in Guisan
et al. (2014), a shortening based on centroid shift, overlap, Overview of functions and data sets
unfilling, and expansion of the niche), to define different
environmental spaces, to correct for occurrence densities and The main functions available in the ecospat package
environmental availability, and to run various tests. support either pre, core- or post-modeling analyses, and

2-EV
 Table 1. Main functions, their description and published references of the ecospat package.

Tool Sub-tool Function Description References where functions has been used
Pre-modelling Spatial ecospat.mantel. Investigates spatial autocorrelation of environmental
autocorrelation correlogram covariables within a set of occurrences as a function of
distance.
Variable selection ecospat.npred Calculates the maximum number of predictors to include Petitpierre et al. (2012)
in the model with a desired correlation between
predictors.
Extrapolation ecospat.exdet Assess climate analogy between a projection extent (proj)
detection and a reference extent (cal, used in general as the
background to calibrate SDMs).
Multivariate ecospat.mess Calculates the MESS (i.e. extrapolation) as in Maxent Broennimann et al. (2014a), Mateo et al. (2015)
environmental (Elith et al. 2010).
similarity surfaces
ecospat.plot.mess Plots the MESS extrapolation index onto the geographical
space.
Phylogenetic diversity ecospat.calculate. Calculates phylogenetic diversity measures (Schweiger Pio et al. (2011, 2014), Ndiribe et al. (2013, 2014)
measures pd et al. 2008).
Biotic interactions ecospat. Computes an index of co-occurrences ranging from 0 Pellissier et al. (2010)
co-occurrence (never co-occurring) to 1 (always co-occuring).
Niche quantification ecospat.grid.clim. Using the scores of an ordination, create a grid z of R ⫻ R Boulangeat et al. (2012a, b), Estrada-Pena et al. (2012, 2013),
dyn pixels (or a vector of R pixels when using scores of Monahan and Tingley (2012), Pellissier et al. (2012), Petitpierre
dimension 1 or SDM predictions) with occurrence et al. (2012), Poudel et al. (2012), Wielstra et al. (2012),
densities. Only scores of one, or two dimensions can Ahmadzadeh et al. (2013a, b), Alexander (2013), Di Febbraro
be used. et al. (2013), Fourcade et al. (2013, 2014), Iwamura et al.
(2013), Koch et al. (2013), Maiorano et al. (2013), Martinez-
Cabrera and Peres-Neto (2013), Strubbe et al. (2013a),
Theodoridis et al. (2013), Werner et al. (2013), Blonder et al.
(2014), Broennimann et al. (2014b), DeChaine et al. (2014),
Early and Sax (2014), Glennon et al. (2014), Goncalves et al.
(2014), Grossenbacher et al. (2014), Grundy et al. (2014), Le
Roux et al. (2014), Li et al. (2014), Montecino et al. (2014),
Pena-Gomez et al. (2014), Russo et al. (2014), Saupe et al.
(2014), Silva et al. (2014), Strubbe and Matthysen (2014),
Tingley et al. (2014), Vences et al. (2014), Stiels et al. (2015)
ecospat.niche. Runs a niche equivalency test (Warren et al. 2008) based
equivalency.test on two species occurrence density grids.
ecospat.niche. Runs a niche similarity test (Warren et al. 2008) based on
similarity.test two species occurrence density grids.
ecospat.plot.niche Plots a niche z created by ecospat.grid.clim.dyn().
ecospat.plot.niche. Plots niche categories and species density.
dyn
ecospat.plot. Plots the contribution of the initial variables to the
contrib analysis (i.e. correlation circle). Typically these are the
eigen vectors and eigen values in ordinations.
ecospat.niche. Calculates the overlap metrics D and I based on two
overlap species occurrence density grids z1 and z2 created by
ecospat.grid.clim.dyn().

(Continued)

3-EV
 Table 1. (Continued)

4-EV
Tool Sub-tool Function Description References where functions has been used
ecospat.plot. Plots a histogram of observed and randomly simulated
overlap.test overlaps, with p-values of equivalency or similarity
tests.
ecospat.niche.dyn. Calculates niche expansion, stability and unfilling.
index
Core niche Model evaluation ecospat.cv.glm K-fold and leave-one-out cross validation for GLM. Randin et al. (2006), Pearman et al. (2008), Engler et al. (2009),
modelling Patsiou et al. (2014)
ecospat.permut.glm A permutation function to get p-values on GLM
coefficients and deviance.
ecospat.cv.gbm K-fold and leave-one-out cross validation for GBM.
ecospat.cv.me K-fold and leave-one-out cross validation for Maxent.
ecospat.cv.rf K-fold and leave-one-out cross validation for
randomForest.
ecospat.boyce Calculates the Boyce index (Hirzel et al. 2006). Petitpierre et al. (2012), Falcucci et al. (2013), Strubbe et al.
(2013b), Suárez-Seoane et al. (2014), Yannic et al. (2014),
Maiorano et al. (2015), Mateo et al. (2015), Thalmann et al.
(2015)

ecospat. Calculates several indices of accuracy of community Dubuis et al. (2013), Pottier et al. (2013), D’Amen et al. (2015a)
CommunityEval predictions.
Spatial predictions ecospat.ESM. Calibrates simple bivariate models as in Lomba et al. Lomba et al. (2010), Breiner et al. (2015)
and projections Modeling (2010) and Breiner et al. (2015).
ecospat.ESM. Evaluates and averages simple bivariate models to ESMs.
EnsembleModeling
ecospat.ESM. Projects simple bivariate models into new space or time.
Projection
ecospat.ESM. Projects calibrated ESMs into new space or time.
EnsembleProjection
ecospat.SESAM.prr Implement the SESAM framework to predict community D’Amen et al. (2015a, b)
composition using a ‘probability ranking’ rule.
ecospat.migclim Enables the implementation of species-specific dispersal Bateman et al. (2013), Normand et al. (2013), Gimona et al.
constraints into projections of species distribution (2015), Lurgi et al. (2015)
models under environmental change and/or landscape
fragmentation scenarios, by calling MigClim functions
(Engler et al. 2012).
Post modelling Variance partition ecospat.varpart Performs variance partitioning for binomial GLM based Randin et al. (2009), Pellissier et al. (2010)
on the deviance of two groups or predicting variables.
Spatial predictions of ecospat.cons_Cscore Tests for non-random patterns of species co-occurrence D’Amen et al. (2015a, b)
species and calculates the C-score index for the whole
assemblages community and for each species pair.
are described in Table 1 and summarized in Supplementary within a given background area. It is important to notice
material Appendix 1, Fig. A1, in relation to modeling func- here that all functions in Broennimann et al. (2007, 2012,
tionalities. The ecospat package contains both novel 2014a), and Petitpierre et al. (2012) are uniquely found in
features (see related section above) and more common func- ecospat (see also Guisan et al. 2014). These functions
tions already used or known elsewhere. In the latter cases, were designed to investigate changes in the niche of invasive
however, most functionalities were not yet made available in species, but they can also be used to compare niches between
a wide-ranging R package, or were improved or implemented sister species (Broennimann et al. 2014b).
in ecospat in a different way than in other packages (e.g. This approach consists of the following steps. First, an
spatial autocorrelation measure, multivariate environmental ordination technique – usually a principal component
similarity surface (MESS), standard species co-occurrence analysis (PCA) – is used to transform n correlated variables
analysis). It also includes basic but useful functionalities for into two uncorrelated linear combinations (principal com-
e.g. preparing the data or projecting the models. ponents) of the original variables. The PCA is calibrated
The package provides four test data sets to run examples: using environmental values from all the pixels of both the
ecospat.testData, ecospat.testTree, native and the invaded study areas. The axes of the PCA
ecospat.testNiche.inv and ecospat.test- thus maximize the ecological variance present in the study
Niche.nat. The first two data sets allow illustrating the areas. Then the PCA scores of the two species distributions,
use of all community modeling functions, whereas the third for which the niches must be compared, are projected
and fourth ones allow illustrating niche quantification, here onto a grid of cells bounded by the minimum and maxi-
in the case of a hypothetical invasive species. The first data mum PCA scores in the study areas. A smoothed density
set, ecospat.testData, contains vegetation plots data of occurrences for each species in each cell of the grid is
with presence records of 50 vascular plant species (angio- then estimated using a kernel density function (in a similar
sperms), a set of associated environmental variables (topo- way as in the package ‘hypervolume’). The global overlap
climatic) and SDM predictions for some species in an area between the niches can be calculated using metrics such as
of ∼ 700 km2 in the western Swiss Alps (< http://rechalpvd. Schoener’s D or Hellinger’s I (see Broennimann et al. 2012
unil.ch >). The second data set, ecospat.testTree, for details of the procedure and metrics). Additionally,
is a phylogenetic tree of class ‘phylo’ that contains data for in the case of invasive species, the niche overlap can be
the same 50 plant species. The third and fourth data sets, disentangled into three categories: unfilling, stability and
ecospat.testNiche.inv and ecospat.test- expansion. This decomposition provides more information
Niche.nat, contain data on geographical coordinates of about the drivers of niche dynamic between native and
vegetation plots, a set of associated environmental variables invaded ranges (Petitpierre et al. 2012, Guisan et al. 2014),
(topo-climatic), occurrence sites for a hypothetical spe- or outside biological invasions, about how two sister spe-
cies and the prediction of its distribution in the native and cies have evolved different niches. Note that one can apply
invaded range. all these analyses on each variable individually (see example
In the next sections, we describe the key and most unique 1), a very useful step to investigate which predictors are
classes of methods implemented in the ecospat package. more conserved (Liu et al. 2016), and much more infor-
For the list of main functions, see Table 1 and the ecospat mative than boxplots comparisons (e.g. Fig. 1 in Mandle
reference manual in CRAN. et al. 2010).

Niche quantification, comparison and tests Ensemble of small models (ESM)

Niche quantification, comparison and tests are key com- Another key feature of the package is the implementation of
ponents of ecospat. Methods for quantifying niches the ESM strategy to model species distributions when a lim-
typically rely on either ordination techniques or on niche- ited numbers of presences or occurrences is available (e.g.
based species distribution models (see Fig. 3 in Guisan et al. for rare species, cryptic species), typically with less than 30
2014). Compared to existing packages such as ‘NicheA’ presence records (Wisz et al. 2008). The ESM approach was
(Qiao et al. 2016), ‘spa’ (Zhang 2013) and ‘SDMTools’ shown to reduce overfitting in the models and accordingly
(VanDerWal et al. 2014) that calculate niche overlap, the improve model accuracy in such cases (Lomba et al. 2010,
methods proposed here were initially the most compre- Breiner et al. 2015). It is based on fitting numerous small
hensive ones (Broennimann et al. 2012). Niche functions models, using two (or possibly three or more) variables at
in the ecospat package provide tools to quantify and each time and finally performs an ensemble of the predic-
compare species niches with an ordination approach. The tions weighted by each submodel performance. It avoids
‘hypervolume’ R package (Blonder et al. 2014) allows some overparameterization of the models, following the sample
niche quantification but without tests. Because niche func- size rule-of thumb of 10:1 subjects to predictors in multiple
tions in ecospat allow for direct comparisons of spe- regression (Harrell 2001). For more details on the method,
cies–environment relationships in environmental space, see Lomba et al. (2010) and Breiner et al. (2015).
no assumptions are required for the model-based approach The ESM strategy can now be run and customized with
(Broennimann et al. 2012). Likewise, the direct comparison four simple functions implemented in the ecospat pack-
can be made for a single variable (see example 1) or multiple age. None of the existing packages for modeling species dis-
variables. Most importantly, the package allows performing tributions (e.g. ‘biomod2’ (Thuiller 2003), ‘sdm’ (Naimi and
niche equivalency and niche similarity tests (Warren et al. Araújo 2016), ‘dismo’ (Hijmans et al. 2016)) currently con-
2008), which assess the statistical significance of a measured tains the ESM approach. The ESM functions in ecospat
niche difference against null model niches taken randomly depend on ‘biomod2’ and allow the selection of the eleven

5-EV
Figure 1. A schematic representation of the workflow used in Example 1. The grey boxes indicate the names of the functions. (A) Plot of
the occupancy of the environment by the species in the native and invaded range. Niche overlap values for native and invaded range, in
terms of Schoener’s D, niche similarity p-values and equivalency p-values via randomization. Niche equivalency test and falling outside the
95% confidence interval of the niche similarity test. (B) Model calibration and projection with Ensemble of Small Models (ESM). (C) Plot
of the predicted to expected ratio along the suitability classes.

different modeling techniques implemented there (and oth- Evaluating model predictions with presence-only data
ers if later added). Various evaluation indices to quantify Additionally, the ecospat package provides the currently
model performance are available: AUC, TSS, Kappa, Somers most exhaustive R implementation of the Boyce index, a pres-
D and Boyce. ESMs are therefore widely applicable and one ence-only and threshold-independent evaluator for SDMs’
or a set of preferred modeling techniques can be selected. predictions (Hirzel et al. 2006). The Boyce index measures
The possibility to use multicore computation enables run- the trend in a plot showing the proportion of presences
ning ESMs in ‘parallel’ to speed up computation time. The across classes of model predictions (similar to the POC-plot
strategy proposed with the ESM functions is as follows: 1) of Phillips and Elith 2010) and allows a conceptually cor-
calibration and evaluation of a series of small models with rect way to evaluate predictions from presence-only models
different techniques (i.e. GLM, Maxent, etc.); 2) averaging like Maxent (see Hirzel et al. 2006 for Boyce, and see Tsoar
them into a single ESM per technique, and finally; 3) pro- et al. 2007, Maher et al. 2014 for presence-only models).
ducing a single ensemble prediction through averaging all Therefore, it provides a more adequate evaluation in cases
single-technique ESMs (i.e. running a double ensembling) where absences do not carry reliable information, typically
(Breiner et al. 2015). Step 3 can be omitted if one does not in cases of non-equilibrium situations, such as with ongoing
want to run a double ensembling, because single-technique invasions by exotic species (Petitpierre et al. 2012), in case
ESMs were shown to perform similarly to ESMs based on a of insufficient detectability of a species (Wintle et al. 2012),
double ensembling (i.e. fitting ESMs with different model- or in the case where the model uses background data and
ling techniques and ensembling them; Breiner et al. 2015). not true absences (e.g. Maxent). The Boyce index measures
The final ESM can then be projected in a different time how much model predictions differ from a random distribu-
period or geographic space. tion of the observed presences across the prediction gradients

6-EV
(Boyce et al. 2002). It is the quantitative equivalent of the Environmentally-constrained species co-occurrence
graphical POC plot (Phillips and Elith 2010) and often sug- analyses
gested to be so far the most appropriate metric in the case of In the context of predicting the spatial distribution of species
presence-only models (Hirzel et al. 2006), although further assemblages, a key challenge is to identify significant species
research should be conducted to assess other possible options interactions (Wisz et al. 2013). In this regard, a novel post-
or possible refinements to evaluate presence-only predictions. modeling methodology of the ecospat package allows
The Boyce index is analogous to a Spearman correlation and running environmentally-constrained pairwise species co-
varies between –1 and ⫹1. Positive values indicate a model occurrence analyses based on a null model approach (as in
in which predictions are consistent with the distribution of Peres-Neto et al. 2001). This functionality uses the func-
presences in the evaluation data set, values close to zero mean tion ecospat.cons_Cscore() and allows to identify
that the model is not different from a random model, and significant patterns of co-occurrence in a presence/absence
negative values indicate counter predictions, i.e. predicting species matrix, classifying each individual pair of species, as
areas where presences are more frequent as being highly suit- random, aggregated, or segregated. The strength of the asso-
able for the species (Hirzel et al. 2006). Originally based on ciation is calculated by the C-score index (Stone and Roberts
predefined (and often arbitrary) classifications of continuous 1990), as Cij ⫽ (Ri – D)(Rj – D), where Cij is the C score for
suitability predictions, the Boyce index can also be threshold species pair i and j, Ri is the row total (the number of species
independent if based on a moving window, also known as occurrences) for species i, Rj is the row total for species j, and
continuous Boyce index (Hirzel et al. 2006). The implemen- D is the number of shared sites in which both species are
tation in ecospat proposes the standardized continuous present (Gotelli 2000). C-scores for individual species pairs
index as default but also allows the user to parameterize pre- and for the whole community are calculated and compared
defined bins. to the statistical expectation for a set of environmentally-
constrained null communities. These are generated follow-
Spatial predictions of communities ing the algorithm Ct-RA1 proposed in Peres-Neto et al.
SDMs have been increasingly used over the last years to (2001), where the species presences are reassigned to sites
predict different community properties, such as species rich- in the matrix according to the relative probability values,
ness, taxonomic composition, functional and phylogenetic maintaining fixed species frequencies (see Peres-Neto et al.
patterns, through stacking of individual species predictions 2001 for further explanations). The probability values in the
(S-SDMs, Dubuis et al. 2011, Calabrese et al. 2014, Cord sites can be obtained by fitting species distribution models
et al. 2014, D’Amen et al. 2015c). A novel approach has (SDMs) for each species. This function can also be run as
been proposed to improve S-SDMs predictions by coupling a pre-modeling analysis if the environmental constraints
them with richness predictions (Guisan and Rahbek 2011, applied to the null models are not based on the SDM pre-
D’Amen et al. 2015a, b). This constraint on S-SDM predic- dictions. This approach is expected to maximize the chance
tions was proposed as one step in the larger ‘spatially-explicit of distinguishing between mutually exclusive processes that
species assemblage modeling’ (SESAM) framework. The lat- may shape species distributions and community assembly, as
ter seeks to reconstruct species assemblages by implementing environment is factored out as a possible explanation for the
the series of successive dispersal, habitat and biotic filters, patterns encountered (Peres-Neto et al. 2001).
and – if needed – the previously described macroecological
constraints, as distinct modeling steps (Guisan and Rahbek Evaluating community predictions
2011). SESAM, as implemented in ecospat, had been When modeling species assemblages (Guisan and Rahbek
tested to date by applying three successive steps (D’Amen 2011), it is known that the accumulation of small errors
et al. 2015a, b): 1) calculating the habitat suitabilities of in individual models can produce larger errors in the
individual species presence for each site – these can typically final predictions (Dubuis et al. 2011, Pottier et al. 2013).
be obtained by fitting an SDM for each species (i.e. S-SDM); Therefore, a robust evaluation of community predictions is
this step represents the application of an environmental filter an important step in spatial ecology analyses (Pineda and
to the community assembly; 2) computing richness predic- Lobo 2009, Pottier et al. 2013). The function ecospat.
tions for each site – the richness prediction can be derived in CommunityEval() allows a list of evaluation metrics
different ways, for instance by summing probabilities from (such as sensitivity, specificity, kappa, TSS, etc.) to be cal-
individual species prediction for each site, or by fitting direct culated for each site at the community level, permitting
macroecological richness models (MEM); this step repre- for a better understanding of the prediction accuracy and
sents the application of a macroecological constraint to the differences, for example in diverse habitat types at various
average number of species that can coexist in the consid- elevations (Mateo et al. 2012, Pottier et al. 2013). The func-
ered unit; 3) applying a biotic rule to decide which species, tion can be applied to the outcomes of any community-level
potentially present in the site, should be retained in the final model, e.g. S-SDM or SESAM predictions.
assemblage prediction to match the predicted richness value.
This last step is implemented in ecospat by the ‘probabil- Phylogenetic diversity measures
ity ranking’ rule (PRR), which consists of ranking the spe- The calculation of phylogenetic diversity measures on
cies in decreasing order of predicted probability of presence observed (pre-modeling) or predicted (post-modeling)
(obtained from the SDMs) and selecting species from the coarse species assemblages (as in Pio et al. 2011, 2014) or
most probable species down the list until the sum of selected local communities (Ndiribe et al. 2013, 2014), provides a
species reaches the expected richness value obtained in step 2 comprehensive list of indices and optimized calculations
(D’Amen et al. 2015a, b). for large data sets (which are then important for producing

7-EV
maps). The function ecospat.calculate.pd() uses next illustrate the implementation of the SESAM framework
a phylogenetic tree, a presence or absence (binary) matrix to predict the spatial distribution of communities, applying a
for each species in each location or grid cell (direct obser- biotic rule to estimate the taxonomic composition of species
vations, S-SDM binary predictions or SESAM predictions), in each site. We then calculate several evaluation indices of
and different phylogenetic diversity measures (Schweiger community-level models, in order to evaluate the commu-
et al. 2008). For each grid cell in the study area, the function nity predictions. Finally, we analyze the spatial patterns of
calculates each of the phylogenetic diversity measures listed geographic overlap in the distributions of the species in the
in Schweiger et al. (2008) based on either tree topology, study area to illustrate the application of the environmental-
branch length or minimum-spanning tree measures. Three ly-constrained null model. Complete scripts to run the two
measures are based on topology (node information only): W examples and to illustrate the use of the ecospat package
(standardized taxic weights), Q (basic taxic weights), and P are provided in the Supplementary material Appendix 1.
(normalized measure of Q). Five measures are based on pair-
wise distances: J (intensive quadratic entropy), F (extensive
quadratic entropy), AvTD (average taxonomic distinctness), Example 1: niche quantification and modeling of an
TTD (total taxonomic distinctness), and Dd (pure diver- invasive species
sity). Finally, three so-called minimum-spanning measures
are based on both branch length and node information: The first example illustrates the use of the functions to inves-
PDroot (clade when root ⫽ FALSE; phylogenetic diversity tigate whether the climatic niche of a one species may differ
with basal branches), PDnode (clade when root ⫽ TRUE; between different geographic areas, and the use of the ESM
phylogenetic diversity), and AvPD (species; average phyloge- functions to model the species’ distribution in the two ranges
netic diversity). These measures can be classified on the basis when few occurrences are available. This niche comparison
of whether they sum (Q, P, W, PDnode, PDroot, F, TTD, framework is typically used to quantify and describe environ-
Dd) or average (AvTD, J, AvPD) the evolutionary history mental niche shifts between the native and exotic ranges of
of all species present in an area. Specific examples using this invasive species or between sister species. The ESM functions
function are found in (Pio et al. 2011, 2014). Different R can be used here to model distributions when the species is not
packages allow similar analyses, but ecospat has a few nov- (yet) very frequent, as typically at the start of an invasion in the
elties and differences here. For example, the ‘picante’ R pack- exotic range. An overview of how the tools used in example 1
age does calculate the same phylogenetic diversity measures can be applied to different case studies is given in Table 1 (see
based on topology and minimum spanning measures, but it also Guisan et al. (2014) and Breiner et al. (2015)).
does not refer to them exactly with the same names (Kembel In this example we first compare the niche of a species
et al. 2010). The package ‘Phylomeasures’ (Tsirogiannis and native in North America and introduced in Australia, then
Sandel 2015) also provides functions to calculate the same calculate a SDM using ESM, and finally evaluate the SDM’s
measures. Unlike ‘picante’ and ‘Phylomeasures’, however, ability to predict the distribution with the Boyce index.
the function ecospat.calculate.pd() includes the The script in Supplementary material Appendix 1 gener-
additional option to calculate five measures based on pair- ates occurrence density, niche overlap, niche equivalency
wise distances (Schweiger et al. 2008). and similarity tests, and uses a brand new tool allowing
the calculation and the visualization of the niche dynamic
as objects in the R environment. The occurrence density
Developed examples with ecospat generated by the function ecospat.grid.clim.
dyn() measures the frequency of occurrences of the spe-
To illustrate the use of the ecospat package, we created cies for each combination of conditions (i.e. each grid
two examples of analytical workflows that can be easily cell) of the environmental space using a kernel smoother.
reproduced by users, involving the main set of functions in The function ecospat.niche.overlap() uses
ecospat. In the first example, we illustrate how to com- the differences in occurrence densities between the two
pare the climatic niche of a hypothetical invasive species in species to measure either Schöner’s D or Hellinger’s I
its native and invaded ranges. We then predict its non-native metrics, both of which range from 0 (no overlap) and 1
distribution by applying an ensemble of small models (ESM) (complete overlap). The functions ecospat.niche.
to avoid model overfitting because only few occurrences exist equivalency.test() and ecospat.niche.
for the invasive species in its invaded range. We finally show similarity.test() perform tests of niche equiva-
how to calculate the Boyce index to evaluate the non-native lency and similarity as described in Warren et al. (2008)
prediction to assess model performance in cases like this but applied in environmental space, where the overlap is
where only presence data are available and absences may not better assessed than in geographical space because it better
be reliable (i.e. here the species is still expanding). takes into account climate availability and analogy between
In the second example we perform a spatial analysis of ranges (Broennimann et al. 2012, Guisan et al. 2014; but
a plant community, examining its phylogenetic and co-oc- note that the last version of the ENMtool by Warren et al.
currence patterns and producing spatial predictions of the (2008) implemented the ecospat approach to test niche
assemblages using novel functions in the ecospat package. differences in environmental space). The niche equivalency
We first compute different phylogenetic diversity measures, test assesses, through random permutations of occurrences
which can then be calculated on the predicted species assem- between ranges, whether the two niches are equivalent. The
blage as well, e.g. to estimate the difference in phylogenetic niche similarity test assesses, through random shifts of the
diversity under current and future climatic scenarios. We niches within available conditions in the study area, whether

8-EV
the species’ niches are more or less similar than expected by Example 2: species assemblage’s structure and
chance. Here, we want to test for niche conservatism, so we spatial predictions
use the alternative ⫽ ‘greater’, i.e. the niche overlap is more
similar than random expectations. Note that in the case of The tools used in example 2 can be applied to provide an
invasive species, we fixed the native niche as a reference and evaluation of the probability of occurrence of different spe-
shifted only the invasive niche (rand.type ⫽ 2). In cases cies within a community, and to model the spatial distribu-
where there are no assumption about a reference niche (e.g. tion of richness and composition in a community, obtaining
sister species), both niches can be simultaneously shifted predictions of species assemblages more realistic than those
(rand.type ⫽ 1). For details on niche equivalency and obtained by stacking individual species predictions. For an
similarity tests, see Warren et al. (2008) and Broennimann overview of different questions that can be addressed using
et al. (2012). Figure 1A provides an example scheme of the the tools of this example, refer to Table 1, Pio et al. (2011,
output from the analysis of a species in the invaded and 2014) and D’Amen et al. (2015a, b). Along this second
native range (complete script available in Supplementary example, we explore different community properties of
material Appendix 1). We see that the niche overlap D is plant assemblages in the Swiss Alps. In particular, we analyze
0.22, and the tests further indicate that the niches of the phylogenetic diversity, we predict plant species richness and
example species in the native and invaded ranges are not composition and we use null model analysis to describe the
more equivalent neither not more similar than expected by spatial pattern of pairwise species co-occurrence. First, we
chance with p-value of 1 and 0.079 respectively. Therefore, illustrate the use of the function ecospat.calculate.
we can conclude that there is no significant climatic niche pd() by using presence–absence data with a phylogenetic
conservatism between native and invaded ranges. The niche tree in order to obtain different phylogenetic diversity mea-
dynamics analysis shows that this niche difference is due sures, but the same function can also be applied to a pre-
to the species’ ability to expand into novel climates in the dicted species assemblage. In the code in the Supplementary
invaded range (niche expansion ⫽ 15%) and to the fact that material Appendix 1, we show how to load the phylogenetic
the species has not (yet) colonized all the climate conditions tree and the data for the example. The resulting object is a list
of the native niche (niche unfilling ⫽ 28%). A visual inspec- of phylogenetic diversity values for each of the grid cells in
tion shows also that an important part of the cold climate the presence/absence matrix, which can be plotted to obtain
conditions available in the native range are no longer avail- a graphical representation of the correlation of phylogenetic
able in the invaded range. diversity with species richness (Fig. 2A). Second, we predict
We next built an ESM (ensemble of small models), by community composition with the function ecospat.
fitting with the ecospat.ESM.Modeling() func- SESAM.prr(). This function implements the final step
tion all possible models that contain only two predictors at of the SESAM assemblage, i.e. the ‘probability ranking rule’
a time out of all those available in the initial set of eight (PRR; D’Amen et al. 2015a). The data required for the func-
variables (from ecospat.testNiche.inv), result- tion are two data frames of: 1) continuous probabilities from
ing in 28 bivariate predictor combinations. Then all bivari- SDMs for all species (in columns) in the considered sites
ate models were combined to an ensemble model (Lomba (in rows), and 2) richness value for each site (first column).
et al. 2010, Breiner et al. 2015) with the ecospat.ESM. In cases where the SDMs were fitted with presence/absence
EnsembleModeling() function. The function first data, thus yielding true probabilities, the richness predic-
evaluates each bivariate model by a cross-validated (repeated tion can then be obtained by summing all raw probabilities
split-sample) evaluation index (in our case AUC), and builds (across species) at each site. Alternatively, a similar predic-
then an average of the 28 models weighted by a cross-vali- tion can be obtained by a separate species richness model.
dation score (in our case Somers’ D, which is AUC * 2 – 1). The function returns a data frame of the same structure as
Models with a Somers’ D higher than 0 (i.e. AUC ⬎ 0.5) the input one, containing the community prediction by the
were selected to build the final ensemble. For more details SESAM framework, i.e. binary predictions for the species in
refer to Breiner et al. (2015). For details see Fig. 1B. each site (D’Amen et al. 2015a). For details see example in
Finally, we evaluate the SDM’s ability to predict the non- Fig. 2B. Third, we calculate community-predictions accu-
native distribution with the Boyce index, a presence only racy with different indices with the function ecospat.
evaluator, with the function ecospat.boyce().This CommunityEval() that requires the observed presence–
function returns the Boyce index value: Spearman.cor absence of the species (eval) and their predictions (pred). The
(in our example 0.889), which is the Spearman correlation returned object is a list of evaluation metrics calculated for
of the PE plot, a graphical plot of the predicted to expected each site (see Fig. 2C for details).
ratio within a moving window along the suitability gradi- Finally, we perform a co-occurrence analysis and test for
ent (Fig. 1C). If the rank of the predicted to expected ratio non-random patterns of species co-occurrences applying an
would be completely ordered along the habitat suitability environmentally-constrained null model (Peres-Neto et al.
gradient, then the Boyce index would be 1. With a score of 2001), with the function ecospat.cons_Cscore().
0.889, we can assume that the SDM is useful to indepen- The format required for the input data is a matrix of plots
dently predict the species invasion in Australia (Fig. 1B). (rows) ⫻ species (columns). Input matrices should have col-
Additionally, such PE plots can be used to reclassify maps umn names (species names) and row names (sampling plots).
of continuous suitability values into discrete categories like The function returns the C-score index for the observed
potential presences or absences, which is particularly useful community (ObsCscoreTot), the mean of C-score for
for conservation managers (Hirzel et al. 2006) or to recali- the simulated communities (SimCscoreTot), the p
brate presence-only SDM’s continuous suitability (Phillips values (PVal.less and PVal.greater) to evaluate
and Elith 2010). the significance of the difference between the former two
9-EV
Figure 2. A schematic representation of the workflow used in Example 2. The grey boxes indicate the names of the functions. (A) Correla-
tion of phylogenetic diversity (pd) with species richness (index). (B) Community prediction by the SESAM framework (i.e. binary
predictions for all species for each site). (C) List of evaluation metrics calculated for each site. (D) Frequency histogram for simulated
C-scores in the null community and frequency histogram for standardized effect size in the observed community: number of species pairs
forming perfect checkerboard distributions.

indices, and finally returns the standardized effect size (SES) Citation of the ecospat package
for the whole community (SES.Tot). A SES greater than 2
or smaller than –2 is statistically significant with a tail prob- The latest available version of the ecospat package
ability of less than 0.05 (Gotelli and McCabe 2002). If a in CRAN is 2.1.1 and the license of the package is GPL
community is structured by competition, we would expect (General Public License is a free, copyleft license for software
the C-score to be large relative to a randomly assembled and other kinds of works). It can be installed from CRAN
community (positive SES). Also, to test for patterns in each using the usual code line as follows:
individual pair of species in the matrix, a table is saved in
the specified working directory. This table contains the ⬎ install.packages(“ecospat”)
observed and expected C-score value and the standardized Scientists using ecospat functions in a published paper
effect size for species pairs and it includes only those pairs should cite this article or the ecospat package directly.
with p values lower than 0.05. A further correction for false Citation information can be obtained by typing:
discovery rate may be needed on the basis of the dataset size
(Gotelli and Ulrich 2010). This option is planned but not ⬎ citation(“ecospat”)
yet implemented in the package. Furthermore, the function on the R command prompt.
automatically plots two frequency histograms: the first one Some of the methods implemented in ecospat have
represents the C-score indices from the simulated communi- been published (Table 1) but their scripts were only available
ties, where the observed C-score value is indicated by the in the related papers and so far not in open access. Users
black solid line (Fig. 2D), and the second histogram rep- should refer to the respective publications for details of the
resents the distribution of the SES values calculated for all methods and interpretation. Users of the package who apply
the unique combination of species pairs in the considered these methods should cite as much as possible the original
community (Fig. 2D). In this case, the observed C-score is sources as appropriate, along with ecospat.
significantly lower than expected by chance (in our example To cite ecospat or acknowledge its use, cite this
the p value is 1 e-4), meaning that the community is domi- Software note as follows, substituting the version of the
nated by positive interactions (aggregated pattern; as can be application that you used for ‘version 0’:
expected in alpine plant communities under stressful condi- Di Cola, V., Broennimann, O., Petitpierre, B., Breiner, F. T.,
tions at high elevations; Callaway et al. 2002). The complete D’Amen, M., Randin, C., Engler, R., Pottier, J., Pio, D.,
script to run the example is available in the Supplementary Dubuis, A., Pellissier, L., Mateo, R. G., Hordijk, W., Salamin,
material Appendix 1. N. and Guisan, A. 2016. ecospat: an R package to support

10-EV
 spatial analyses and modeling of species niches and distribu- models: a case study for Mexican pines. – J. Biogeogr. 41:
tions. – Ecography 39: 000–000 (ver. 0). 736–748.
D’Amen, M. et al. 2015a. Using species richness and functional
traits predictions to constrain assemblage predictions from
Acknowledgements – We would like to thank the reviewers, Nick stacked species distribution models. – J. Biogeogr. 42:
Golding, and Matheus Lima-Ribeiro for their insightful comments 1255–1266.
on the paper, as these comments led us to an improvement of the D’Amen, M. et al. 2015b. Predicting richness and composition
work. Functions included in this package were developed in several in mountain insect communities at high resolution: a new
projects mainly funded by the Swiss National Science Foundation test of the SESAM framework. – Global Ecol. Biogeogr. 24:
(SNF; grants nr 31003A-125145, 31003A-152866, 1443–1453.
CR23I2_162754), the SNF National Center of Competence in D’Amen, M. et al. 2015c. Spatial predictions at the community
Research (NCCR) ‘Plant Survival’, and the European Commission’s level: from current approaches to future frameworks. – Biol.
sixth framework project ‘EcoChange’. MD was supported by the Rev. doi: 10.1111/brv.12222
FP7-PEOPLE-2012-IEF Marie Curie Action (project nr 327987). Dawson, M. N. et al. 2013. An horizon scan of biogeography.
RGM was supported by the FP7-PEOPLE-2013-IEF Marie Curie – Front. Biogeogr. 5: fb_18854.
Action (project nr 622620). FB was supported by the Federal Office DeChaine, E. G. et al. 2014. Integrating environmental, molecular,
for the Environment (FTB) (‘Modeling and Evaluation of Species and morphological data to unravel an ice-age radiation of
Distributions’, project number 8T20/10.0042.PJ). arctic-alpine Campanula in western North America. – Ecol.
Evol. 4: 3940–3959.
Di Febbraro, M. et al. 2013. The use of climatic niches in screening
procedures for introduced species to evaluate risk of spread: a
References case with the american eastern grey squirrel. – PLoS One 8:
e66559.
Ahmadzadeh, F. et al. 2013a. Rapid lizard radiation lacking niche Dray, S. and Dufour, A. B. 2007. The ade4 package: implementing
conservatism: ecological diversification within a complex the duality diagram for ecologists. – J. Stat. Softw. 22: 1–20.
landscape. – J. Biogeogr. 40: 1807–1818. Dubuis, A. et al. 2011. Predicting spatial patterns of plant
Ahmadzadeh, F. et al. 2013b. Cryptic speciation patterns in Iranian species richness: a comparison of direct macroecological and
rock lizards uncovered by integrative taxonomy. – PLoS One species stacking modelling approaches. – Divers. Distrib. 17:
8: e80563. 1122–1131.
Alexander, J. M. 2013. Evolution under changing climates: climatic Dubuis, A. et al. 2013. Improving the prediction of plant species
niche stasis despite rapid evolution in a non-native plant. distribution and community composition by adding edaphic
– Proc. R. Soc. B 280: 20131446. to topo-climatic variables. – J. Veg. Sci. 24: 593–606.
Bateman, B. L. et al. 2013. Appropriateness of full‐, partial‐and Dullinger, S. et al. 2012. Extinction debt of high-mountain plants
no‐dispersal scenarios in climate change impact modelling. under twenty-first-century climate change. – Nat. Clim.
– Divers. Distrib. 19: 1224–1234. Change 2: 619–622.
Blonder, B. et al. 2014. The n‐dimensional hypervolume. – Global Dungan, J. L. et al. 2002. A balanced view of scale in spatial
Ecol. Biogeogr. 23: 595–609. statistical analysis. – Ecography 25: 626–640.
Boulangeat, I. et al. 2012a. Accounting for dispersal and biotic Early, R. and Sax, D. F. 2014. Climatic niche shifts between species’
interactions to disentangle the drivers of species distributions native and naturalized ranges raise concern for ecological
and their abundances. – Ecol. Lett. 15: 584–593. forecasts during invasions and climate change. – Global Ecol.
Boulangeat, I. et al. 2012b. Improving plant functional groups for Biogeogr. 23: 1356–1365.
dynamic models of biodiversity: at the crossroads between Elith, J. et al. 2010. The art of modelling range-shifting species.
functional and community ecology. – Global Change Biol. 18: – Methods Ecol. Evol. 1: 330–342.
3464–3475. Engler, R. et al. 2009. Predicting future distributions of mountain
Boyce, M. S. et al. 2002. Evaluating resource selection functions. plants under climate change: does dispersal capacity matter?
– Ecol. Model. 157: 281–300. – Ecography 32: 34–45.
Breiner, F. T. et al. 2015. Overcoming limitations of modelling rare Engler, R. et al. 2012. The MIGCLIM R package – seamless
species by using ensembles of small models. – Methods Ecol. integration of dispersal constraints into projections of species
Evol. 6: 1210–1218. distribution models. – Ecography 35: 872–878.
Broennimann, O. et al. 2007. Evidence of climatic niche shift dur- Estrada-Pena, A. et al. 2012. Occurrence patterns of Afrotropical
ing biological invasion. – Ecol. Lett. 10: 701–709. ticks (Acari: Ixodidae) in the climate space are not correlated
Broennimann, O. et al. 2012. Measuring ecological niche overlap with their taxonomic relationships. – PLoS One 7: e36976.
from occurrence and spatial environmental data. – Global Ecol. Estrada-Pena, A. et al. 2013. Association of environmental traits
Biogeogr. 21: 481–497. with the geographic ranges of ticks (Acari: Ixodidae) of medical
Broennimann, O. et al. 2014a. Contrasting spatio-temporal and veterinary importance in the western Palearctic. A digital
climatic niche dynamics during the eastern and western inva- data set. – Exp. Appl. Acarol. 59: 351–366.
sions of spotted knapweed in North America. – J. Biogeogr. Falcucci, A. et al. 2013. Modeling the potential distribution for a
41: 1126–1136. range-expanding species: wolf recolonization of the Alpine
Broennimann, O. et al. 2014b. Influence of climate on the presence range. – Biol. Conserv. 158: 63–72.
of colour polymorphism in two montane reptile species. – Biol. Fitzpatrick, M. C. et al. 2007. The biogeography of prediction
Lett. 10: 20140638. error: why does the introduced range of the fire ant over-predict
Calabrese, J. M. et al. 2014. Stacking species distribution models its native range? – Global Ecol. Biogeogr. 16: 24–33.
and adjusting bias by linking them to macroecological models. Fourcade, Y. et al. 2013. Confronting expert-based and modelled
– Global Ecol. Biogeogr. 23: 99–112. distributions for species with uncertain conservation status: a
Callaway, R. M. et al. 2002. Positive interactions among alpine case study from the corncrake (Crex crex). – Biol. Conserv. 167:
plants increase with stress. – Nature 417: 844–848. 161–171.
Cord, A. F. et al. 2014. Remote sensing data can improve Fourcade, Y. et al. 2014. Mapping species distributions with
predictions of species richness by stacked species distribution MAXENT using a geographically biased sample of presence

11-EV
 data: a performance assessment of methods for correcting sam- Levin, S. A. 1992. The problem of pattern and scale in
pling bias. – PLoS One 9: e97122. ecology: the Robert H. MacArthur Award Lecture. – Ecology
Franklin, J. 2009. Mapping species distribution: spatial inference 73: 1943–1967.
and prediction. – Cambridge Univ. Press. Li, Y. et al. 2014. Residence time, expansion toward the equator
Gimona, A. et al. 2015. Habitat networks and food security: pro- in the invaded range and native range size matter to climatic
moting species range shift under climate change depends on niche shifts in non-native species. – Global Ecol. Biogeogr. 23:
life history and the dynamics of land use choices. – Landscape 1094–1104.
Ecol. 30: 771–789. Liu, X. et al. 2016. Realized climatic niches are conserved along
Glennon, K. L. et al. 2014. Evidence for shared broad-scale cli- maximum temperatures among herpetofaunal invaders. – J.
matic niches of diploid and polyploid plants. – Ecol. Lett. 17: Biogeogr. doi: 10.1111/jbi.12808
574–582. Lomba, A. et al. 2010. Overcoming the rare species modelling
Goncalves, E. et al. 2014. Global invasion of Lantana camara: has paradox: a novel hierarchical framework applied to an Iberian
the climatic niche been conserved across continents? – PLoS endemic plant. – Biol. Conserv. 143: 2647–2657.
One 9: e111468. Lortie, C. J. et al. 2004. Rethinking plant community theory.
Gotelli, N. J. 2000. Null model analysis of species co-occurrence – Oikos 107: 433–438.
patterns. – Ecology 81: 2606–2621. Lurgi, M. et al. 2015. Modelling range dynamics under global
Gotelli, N. J. and McCabe, D. J. 2002. Species co-occurrence: a change: which framework and why? – Methods Ecol. Evol. 6:
meta-analysis of J. M. Diamond’s assembly rules model. – Ecol- 247–256.
ogy 83: 2091–2096. Maguire, K. C. et al. 2015. Modeling species and community
Gotelli, N. J. and Ulrich, W. 2010. The empirical Bayes approach responses to past, present, and future episodes of climatic
as a tool to identify non-random species associations. – Oeco- and ecological change. – Annu. Rev. Ecol. Evol. Syst. 46:
logia 162: 463–477. 343–368.
Grossenbacher, D. L. et al. 2014. Niche and range size patterns Maher, S. P. et al. 2014. Pattern‐recognition ecological niche
suggest that speciation begins in small, ecologically diverged models fit to presence‐only and presence–absence data.
populations in North American monkeyflowers (Mimulus – Methods Ecol. Evol. 5: 761–770.
spp.). – Evolution 68: 1270–1280. Maiorano, L. et al. 2013. Building the niche through time: using
Grundy, J. P. B. et al. 2014. Testing multiple pathways for impacts 13,000 years of data to predict the effects of climate change on
of the non-native black-headed weaver Ploceus melanocephalus three tree species in Europe. – Global Ecol. Biogeogr. 22:
on native birds in Iberia in the early phase of invasion. – Ibis 302–317.
156: 355–365. Maiorano, L. et al. 2015. Modeling the distribution of Apennine
Guisan, A. and Zimmermann, N. E. 2000. Predictive habitat dis- brown bears during hyperphagia to reduce the impact of wild
tribution models in ecology. – Ecol. Model. 135: 147–186. boar hunting. – Eur. J. Wildl. Res. 61: 241–253.
Guisan, A. and Rahbek, C. 2011. SESAM – a new framework Mandle, L. et al. 2010. Conclusions about niche expansion in
integrating macroecological and species distribution models for introduced Impatiens walleriana populations depend on
predicting spatio-temporal patterns of species assemblages. – J. method of analysis. – PLoS One 5: e15297.
Biogeogr. 38: 1433–1444. Martinez-Cabrera, H. I. and Peres-Neto, P. R. 2013. Shifts in cli-
Guisan, A. et al. 2013. Predicting species distributions for conser- mate foster exceptional opportunities for species radiation: the
vation decisions. – Ecol. Lett. 16: 1424–1435. case of South African Geraniums. – PLoS One 8: e83087.
Guisan, A. et al. 2014. Unifying niche shift studies: insights from Mateo, R. G. et al. 2012. Do stacked species distribution models
biological invasions. – Trends Ecol. Evol. 29: 260–269. reflect altitudinal diversity patterns? – PLoS One 7: e32586.
Harrell, F. E. 2001. Regression modeling strategies: with applica- Mateo, R. G. et al. 2015. What is the potential of spread in invasive
tions to linear models, logistic regression, and survival analysis. bryophytes? – Ecography 38: 480–487.
– Springer. Mateo, R. G. et al. 2016. The mossy north: an inverse latitudinal
Hastings, A. et al. 2011. Spatial ecology across scales. – Biol. Lett. diversity gradient in European bryophytes. – Sci. Rep. 6:
7: 163–165. 25546.
Hijmans, R. J. 2015. raster: geographic data analysis and modeling. Medley, K. A. 2010. Niche shifts during the global invasion of the
– R package ver. 2.5-2, < http://CRAN.R-project.org/ Asian tiger mosquito, Aedes albopictus Skuse (Culicidae),
package = raster >. revealed by reciprocal distribution models. – Global Ecol. Bio-
Hijmans, R. J. et al. 2016. dismo: species distribution modeling. geogr. 19: 122–133.
– R package ver. 1.0-15, < http://CRAN.R-project.org/ Monahan, W. B. and Tingley, M. W. 2012. Niche tracking and
package = dismo >. rapid establishment of distributional equilibrium in the house
Hirzel, A. H. et al. 2001. Assessing habitat-suitability models with sparrow show potential responsiveness of species to climate
a virtual species. – Ecol. Model. 145: 111–121. change. – PLoS One 7: e42097.
Hirzel, A. H. et al. 2006. Evaluating the ability of habitat suitabil- Montecino, V. et al. 2014. Niche dynamics and potential geographic
ity models to predict species presences. – Ecol. Model. 199: distribution of Didymosphenia geminata (Lyngbye) M. Schmidt,
142–152. an invasive freshwater diatom in southern Chile. – Aquat.
Iwamura, T. et al. 2013. How robust are global conservation Invasions 9: 507–519.
priorities to climate change? – Global Environ. Change 23: Morlon, H. et al. 2011. Spatial patterns of phylogenetic diversity.
1277–1284. – Ecol. Lett. 14: 141–149.
Kembel, S. W. et al. 2010. Picante: R tools for integrating phylog- Naimi, B. and Araújo, M. B. 2016. sdm: a reproducible and
enies and ecology. – Bioinformatics 26: 1463–1464. extensible R platform for species distribution modelling.
Koch, C. et al. 2013. Two new endemic species of Ameiva – Ecography 39: 368–375.
(Squamata: Teiidae) from the dry forest of northwestern Peru Ndiribe, C. et al. 2013. Phylogenetic plant community
and additional information on Ameiva concolor Ruthven, 1924. structure along elevation is lineage specific. – Ecol. Evol. 3:
– Zootaxa 3745: 263–295. 4925–4939.
Le Roux, J. J. et al. 2014. Relatedness defies biogeography: the tale Ndiribe, C. et al. 2014. Plant functional and phylogenetic turnover
of two island endemics (Acacia heterophylla and A. koa). – New correlate with climate and land use in the western Swiss Alps.
Phytol. 204: 230–242. – J. Plant. Ecol. 7: 439–450.

12-EV
Nogués-Bravo, D. 2009. Predicting the past distribution of species’ Saupe, E. E. et al. 2014. Macroevolutionary consequences of
climatic niches. – Global Ecol. Biogeogr. 18: 521–531. profound climate change on niche evolution in marine molluscs
Normand, S. et al. 2013. A greener Greenland? Climatic potential over the past three million years. – Proc. R. Soc. B 281:
and long-term constraints on future expansions of trees and 20141995.
shrubs. – Phil. Trans. R. Soc. B 368: 20120479. Schweiger, O. et al. 2008. A comparative test of phylogenetic
Paradis, E. et al. 2004. APE: analyses of phylogenetics and evolu- diversity indices. – Oecologia 157: 485–495.
tion in R language. – Bioinformatics 20: 289–290. Silva, D. P. et al. 2014. Using ecological niche models and niche
Patsiou, T. S. et al. 2014. Topo‐climatic microrefugia explain the analyses to understand speciation patterns: the case of sister
persistence of a rare endemic plant in the Alps during the last Neotropical orchid bees. – PLoS One 9: e113246.
21 millennia. – Global Change Biol. 20: 2286–2300. Skidmore, A. K. et al. 2011. Geospatial tools address emerging
Pearman, P. B. et al. 2008. Niche dynamics in space and time. issues in spatial ecology: a review and commentary on the
– Trends Ecol. Evol. 23: 149–158. special issue. – Int. J. Geogr. Inform. Sci. 25: 337–365.
Pellissier, L. et al. 2010. Species distribution models reveal apparent Soberon, J. 2007. Grinnellian and Eltonian niches and geographic
competitive and facilitative effects of a dominant species on the distributions of species. – Ecol. Lett. 10: 1115–1123.
distribution of tundra plants. – Ecography 33: 1004–1014. Stiels, D. et al. 2015. Niche shift in four non-native estrildid
Pellissier, L. et al. 2012. High host-plant nitrogen content: a finches and implications for species distribution models. – Ibis
prerequisite for the evolution of ant–caterpillar mutualism? – J. 157: 75–90.
Evol. Biol. 25: 1658–1666. Stone, L. and Roberts, A. 1990. The checkerboard score and species
Pellissier, L. et al. 2013a. Phylogenetic alpha and beta diversities of distributions. – Oecologia 85: 74–79.
butterfly communities correlate with climate in the western Strubbe, D. and Matthysen, E. 2014. Patterns of niche conserva-
Swiss Alps. – Ecography 36: 541–550. tism among non-native birds in Europe are dependent on
Pellissier, L. et al. 2013b. Turnover of plant lineages shapes herbivore introduction history and selection of variables. – Biol. Inva-
phylogenetic beta diversity along ecological gradients. – Ecol. sions 16: 759–764.
Lett. 16: 600–608. Strubbe, D. et al. 2013a. Niche conservatism in non-native birds
Pena-Gomez, F. T. et al. 2014. Climatic niche conservatism and in Europe: niche unfilling rather than niche expansion.
biogeographical non-equilibrium in Eschscholzia californica – Global Ecol. Biogeogr. 22: 962–970.
(Papaveraceae), an invasive plant in the Chilean Mediterranean Strubbe, D. et al. 2013b. Niche conservatism in non‐native birds
region. – PLoS One 9: e105025. in Europe: niche unfilling rather than niche expansion.
Peres-Neto, P. R. et al. 2001. Environmentally constrained null – Global Ecol. Biogeogr. 22: 962–970.
models: site suitability as occupancy criterion. – Oikos 93: Suárez‐Seoane, S. et al. 2014. Scaling of species distribution models
110–120. across spatial resolutions and extents along a biogeographic
Peterson, A. T. et al. 2011. Ecological niches and geographic gradient. The case of the Iberian mole Talpa occidentalis.
distributions. – Princeton Univ. Press. – Ecography 37: 279–292.
Petitpierre, B. et al. 2012. Climatic niche shifts are rare among Thalmann, D. J. K. et al. 2015. Areas of high conservation value
terrestrial plant invaders. – Science 335: 1344–1348. in Georgia: present and future threats by invasive alien plants.
Phillips, S. J. and Elith, J. 2010. POC plots: calibrating species – Biol. Invasions 17: 1041–1054.
distribution models with presence-only data. – Ecology 91: Theodoridis, S. et al. 2013. Divergent and narrower climatic niches
2476–2484. characterize polyploid species of European primroses in Primula
Phillips, S. J. et al. 2006. Maximum entropy modeling of species sect. Aleuritia. – J. Biogeogr. 40: 1278–1289.
geographic distributions. – Ecol. Model. 190: 231–259. Thuiller, W. 2003. BIOMOD – optimizing predictions of species
Pineda, E. and Lobo, J. M. 2009. Assessing the accuracy of species distributions and projecting potential future shifts under global
distribution models to predict amphibian species richness change. – Global Change Biol. 9: 1353–1362.
patterns. – J. Anim. Ecol. 78: 182–190. Thuiller, W. et al. 2009. BIOMOD – a platform for ensemble
Pio, D. V. et al. 2011. Spatial predictions of phylogenetic forecasting of species distributions. – Ecography 32:
diversity in conservation decision making. – Conserv. Biol. 25: 369–373.
1229–1239. Thuiller, W. et al. 2016. biomod2: ensemble platform for species
Pio, D. V. et al. 2014. Climate change effects on animal and plant distribution modeling. – R package ver. 3.3-7, < http://CRAN.
phylogenetic diversity in southern Africa. – Global Change Rproject.org/package = biomod2 >.
Biol. 20: 1538–1549. Tingley, R. et al. 2014. Realized niche shift during a global
Pottier, J. et al. 2013. The accuracy of plant assemblage prediction biological invasion. – Proc. Natl Acad. Sci. USA 111:
from species distribution models varies along environmental 10233–10238.
gradients. – Global Ecol. Biogeogr. 22: 52–63. Tsirogiannis, C. and Sandel, B. 2015. PhyloMeasures: fast and
Poudel, R. C. et al. 2012. Using morphological, molecular and exact algorithms for computing phylogenetic biodiversity
climatic data to delimitate yews along the Hindu Kush-Himalaya measures. – R package ver. 1.1, < http://CRAN.R-project.org/
and adjacent regions. – PLoS One 7: e46873. package = PhyloMeasures >.
Pulliam, H. R. 2000. On the relationship between niche and Tsoar, A. et al. 2007. A comparative evaluation of presence-only
distribution. – Ecol. Lett. 3: 349–361. methods for modelling species distribution. – Divers. Distrib.
Qiao, H. et al. 2016. NicheA: creating virtual species and ecological 13: 397–405.
niches in multivariate environmental scenarios. – Ecography VanDerWal, J. et al. 2014. SDMTools: species distribution model-
39: 805–813. ling tools: tools for processing data associated with species
Randin, C. F. et al. 2006. Are niche-based species distribution distribution modelling exercises. – R package ver. 1.1-20,
models transferable in space? – J. Biogeogr. 33: 1689–1703. < http://CRAN.R-project.org/package = SDMTools >.
Randin, C. F. et al. 2009. Land use improves spatial predictions of Vences, M. et al. 2014. New insights on phylogeography and
plant abundances but not occurrences in a mountain landscape. distribution of painted frogs (Discoglossus) in northern Africa
– J. Veg. Sci. 20: 996–1008. and the Iberian Peninsula. – Amphibia-Reptilia 35: 305–320.
Russo, D. et al. 2014. What story does geographic separation of Warren, D. L. et al. 2008. Environmental niche equivalency versus
insular bats tell? A case study on Sardinian Rhinolophids. conservatism: quantitative approaches to niche evolution.
– PLoS One 9: e110894. – Evolution 62: 2868–2883.

13-EV
 Werner, P. et al. 2013. The role of climate for the range limits of Wisz, M. S. et al. 2008. Effects of sample size on the performance
parapatric European land salamanders. – Ecography 36: of species distribution models. – Divers. Distrib. 14: 763–773.
1127–1137. Wisz, M. S. et al. 2013. The role of biotic interactions in shaping
Wielstra, B. et al. 2012. Corresponding mitochondrial DNA and distributions and realised assemblages of species: implications
niche divergence for crested newt candidate species. – PLoS for species distribution modelling. – Biol. Rev. 88: 15–30.
One 7: e46671. Yannic, G. et al. 2014. Genetic diversity in caribou linked to past
Wintle, B. A. et al. 2012. Designing occupancy surveys and inter- and future climate change. – Nat. Clim. Change 4: 132–137.
preting non-detection when observations are imperfect. Zhang, J. 2013. spaa: species association analysis. – R package ver.
– Divers. Distrib. 18: 417–424. 0.2.1, < http://CRAN.R-project.org/package = spaa >.

Supplementary material (Appendix ECOG-02671 at < www.

ecography.org/appendix/ecog-02671 >). Appendix 1.

14-EV

View publication stats