Can habitat suitability estimated from MaxEnt predict colonizations and extinctions?

MaxEnt has been widely used to model species’ geographic distributions as functions of environmental variables and to predict changes in distributions in response to environmental change. Here, we test the predictive ability of MaxEnt models through time by modelling colonizations and extinctions.


| INTRODUC TI ON
Species distribution models (SDMs), fundamentally, are correlative models that relate species' presence at particular geographic locations to environmental characteristics of those locations (Warren, 2012). Most commonly, those characteristics are variables that plausibly define aspects of a species' Hutchinsonian niche Warren, 2012) such as temperature and precipitation.
It is frequently proposed that species' ranges shift in response to climate change in order to keep the species within their climatic niches Menéndez et al., 2006;Parmesan & Yohe, 2003). SDMs have been used to predict those range shifts (Bellard et al., 2012;Gallon et al., 2014). But to infer that an SDM reflects a species' niche (albeit imperfectly), one must assume that species' environmental niches are constant through time and space (i.e. the assumption of stationarity) (Pearman et al., 2008;Wiens et al., 2010). One must assume further that ranges are near equilibrium with their driving environmental variables (i.e. suitable habitats are mainly occupied, while unsuitable habitats are unoccupied).
The assumptions of climate-based SDMs might not be met under climate change (Yackulic et al., 2015). Environmental conditions within the geographic range currently occupied by a species only approximately reflect the species' physiological tolerances (Araújo et al., 2013;Boucher-Lalonde et al., 2014a;Hargreaves et al., 2014; but see Lee-Yaw et al., 2016). Geographic ranges potentially can be affected by other factors such as biotic interactions (Davis et al., 1998;Svenning et al., 2014), dispersal (Marsico & Hellmann, 2009), macroevolutionary history (Roy et al., 2009), stochastic factors (Rosindell et al., 2011), etc. Insofar as non-climatic factors influence species' ranges, SDMs based on realized distributions of species in climate space will be poor predictors of species' responses to climate change (Anderson, 2013). In principle, laboratory-measured physiological tolerance should be better predictors of responses to climate change (Phillips, 2008); however, this is not necessarily true (e.g. reptiles in Figure 5 in Araújo et al., 2013). For want of a better alternative, environmental niche models have often been used to predict species' responses to climate change (Matthews et al., 2011;Reese & Skagen, 2017) up to several decades in the future (Thomas et al., 2004;Thuiller, 2004).
In this study, we ask how well the SDM MaxEnt (Elith et al., 2011;Phillips et al., 2006) predicts changes in species' ranges as climatic variables change. With > 6,000 published applications (Phillips et al., 2017), MaxEnt is among the most used (and misused : Yackulic et al., 2013) SDM techniques, principally because it can be used with either presence-only data or presence-absence data. MaxEnt has been applied not only to potential effects of climate change on species' ranges (Yates et al., 2010;Brotons 2014) but also the potential ranges of introduced species (Wang et al., 2007). Although many authors prefer SDMs based on presence-absence data, Rapacciuolo et al. (2012) found that, over a wide range of species, MaxEnt performed as well, on average, as presence/absence methods. See Appendix S1 for further advantages of MaxEnt.
To predict changes in species' ranges when climate changes, an SDM must capture a causal link between ranges and climate. Yet, SDMs may appear to perform well, even without that link. Spatially autocorrelated variables such as environmental variables and species ranges tend to be correlated, even when no causal link exists (e.g. in simulated data : Chapman, 2010;Currie et al., 2020;Fourcade et al., 2018). The performance of SDMs is typically assessed by cross-validating through space: often, data are randomly partitioned into a calibration subset and a subset against which the predictions of the calibrated model are tested. Spatial autocorrelation unavoidably inflates measures of model performance in these spatial crossvalidations. Using test data from independent geographic regions provides stronger tests (Wenger & Olden, 2012;Radosavljevic & Anderson, 2014; but see Currie et al., 2020). Better still, to test whether a correlation is causal, one can use a natural experiment: observe whether the dependent variable changes as predicted when the independent variable changes through time. If SDMs capture causal relationships, then the SDMs should accurately predict temporal changes in species' ranges when hypothesized driving variables change (e.g. Kharouba et al., 2009). It should be possible, first, to relate geographic variation of a species' presence to variation in climatic variables using an SDM and then to predict where colonizations and extinctions will occur when climate changes. Sofaer et al. (2018) recently asked whether species distribution models accurately predict which species will gain or lose range area under climate change. They examined the ranges of 190 songbirds species within the conterminous United States between 1977 and 2014, using the North American Breeding Bird Survey (BBS) (Pardieck et al., 2017). They found that SDMs calibrated on data from 1977-1979 showed a "striking inability" to predict change in range size in 2012-2014. Yet, poor predictions may have resulted from two methodological issues. First, Sofaer et al. used only the portions of species ranges that fell within the conterminous United States. This potentially excluded parts of some species' environmental niches. Second, the study included small-ranged species whose ranges are unlikely to be constrained by climatic variables (Proosdij et al., 2016;Stockwell & Peterson, 2002).
The purpose of the present study is to test whether changes in climatic suitability, as estimated from an SDM, accurately predict local colonizations and extinctions. To do this, we examined the colonization and extinction dynamics from 1979 to 2009 of 21 passerine bird species whose ranges were reasonably large and were entirely sampled in the BBS. We tested the predictions that MaxEntestimated changes in habitat suitability resulting from changes in environmental variables are 1) positively related to changes in colonization probability, and 2) negatively related to changes in extinction probability over a relatively short time scale (one to four years), and over 30 years. Then, we tested whether the ability of MaxEnt models to predict temporal variation in occupancy relates to how well they describe spatial variation in occupancy. Finally, we compared the predictive ability of MaxEnt to that of a simple spatial model based on occupancy of neighbouring sites by conspecifics. We ask how much temporal variability in species occupancy is related to environmental variables, beyond what can be explained by the spatial structure of species' distributions.

| MATERIAL AND ME THODS
An outline summarizing the following methods can be found at the beginning of Appendix S1.

| Species data
Species range data came from the North American Breeding Bird Survey (Pardieck et al., 2017). The data consist of observations of bird species' presence and abundance on ~ 4,900 routes in the United States and southern Canada during birds' breeding season (May to June). Most routes are sampled repeatedly. Sampling intensity is strongly related to human population density and road density.
We used BBS data from 1979 to 2009 because our environmental data only covered that range.
The BBS contains observations on > 700 species, but we used only 21 species for the following reasons. First, we only included Passerines, the species most reliably detected using the BBS protocol (Cunningham et al., 1999). Second, we excluded non-native species. Third, to observe distribution shifts in all directions, we included only species whose entire breeding range fell within the area sampled by the BBS. Finally, we excluded species that were observed on < 20 routes during any given year because very small ranges are unlikely to be constrained by climate (Proosdij et al., 2016;Stockwell & Peterson, 2002), and because SDMs based on < 20 points have low precision. This left 21 species (Appendix S2). Range shifts have been reported over the study period for those 21 species (Currie & Venne, 2017).
Following Currie and Venne (2017), we first defined each species' potential breeding range as the union of its breeding range map depicted by BirdLife International (http://www.birdl ife.org/dataz one/info/spcdo wnload), and the BBS routes at which the species was observed at least once between 1979 and 2009, plus a 100 km buffer around each of those routes (not all BBS presences fall within the BirdLife ranges). The 100 km buffer is arbitrary, but it was intended to capture areas that would be within dispersal distance of occupied routes. Paradis et al. (1998;their Figure 1) found that 82%-99% of the breeding and natal dispersal distances of bird species in Britain were within 100 km.

| Environmental data
Monthly environmental data came from the National Center for Environmental Prediction (Saha et al., 2010) at a spatial resolution of 0.5° (~ 50 km). We calculated the average of each environmental variable during the breeding season (May to June) for each year from 1979 to 2009. We used seven environmental variables in the MaxEnt models: average temperature, average precipitation rate, average transpiration, soil type, vegetation type, average per cent vegetation cover, and elevation (Appendix S3).

| MaxEnt models
MaxEnt estimates the relative suitability of sites for a given species by comparing environmental conditions at known occupied sites to the available environmental conditions in the study region (the "background"). We fitted MaxEnt models (Phillips et al., 2017) using the dismo package (Hijmans et al., 2015), in R version 3.4.0 (R Core Team, 2018) with the default settings (i.e. all feature classes except threshold). MaxEnt output (hereafter, habitat suitability) is proportional to the density of individuals per unit area (raw output), and it is monotonically related to the probably of presence (logistic output) (Elith et al., 2011;Yackulic et al., 2013). Some studies convert suitability into a binary variable (suitable/unsuitable); here, we treated suitability as a continuous variable.
The background should include "the area that has been accessible to the species of interest over relevant time period" (Barve et al., 2011). The area included in the background can strongly influence MaxEnt model fit (Merow et al., 2013;Royle et al., 2012).
As species' potential ranges were intended to define the area that a species could potentially colonize during our study, we sampled that area for background conditions. Given that BBS sample sites are not distributed randomly in space, we sampled the environmental conditions at all BBS routes within the potential breeding range, following Thibaud et al. (2014). This way, presence and background points are subject to the same sampling bias (Elith et al., 2011;Phillips et al., 2009).
For purposes of fitting MaxEnt models, we used species' presences in the BBS data. For purposes of evaluating the models' predictive ability (and for defining colonizations and extinctions), we followed Phillips et al. (2006) and Elith et al. (2011), and we treated the BBS data as presence/absence. Thus, we test a model's ability to distinguish between sites where a species is observed, versus sites where it is not observed, either because it is absent, or too rare to be detected. We assessed MaxEnt model fit using AUC (Area Under the Curve). AUC measures the ability of the model to discriminate between sites where a species is present, versus sites where it is absent. Random association of habitat suitability and species presence would yield AUC = 0.5, while a perfect model would produce AUC = 1. See Appendix S1 for discussion of limitations of AUC.
Following earlier literature, we used cross-validation to calculate model AUC. The ability of SDMs to predict distributions through space is typically tested by (randomly) allocating sampling locations to calibration and testing subsets (Araújo et al., 2005). The hypothesis that species' ranges track climatic conditions through time assumes stationarity and equilibrium (see above); the year when a presence was observed should not be important to the model fit. Therefore, in this study, we test predictive ability through time by randomly allocated sampling years to calibration and testing subsets. When a species was present on a given route during ≥ 1 year between 1979 and 2009, we randomly selected one presence-year, and we recorded the environmental conditions that year. We similarly generated a sample of environmental conditions during absence-years.
To produce the MaxEnt background sample (i.e. training background data), we combined the samples of absences and presences because the background represents the range of conditions available in the potential range. AUC was calculated using all the presences and absences that were not used during model fitting. This produced a testing dataset that was > 2X larger than the calibration dataset for all species. Because the number of presences and absences can affect AUC (Yackulic et al., 2013), we resampled the data to have equal numbers of absences and presences before calculating AUC. We repeated this process 10 times, and site suitability was then estimated as the average of the 10 models' predictions, weighting each model based on its AUC: where w i is the weight of the i th model, and a i is AUC value of that model.

| Colonization and extinction models
After fitting the MaxEnt models, we used logistic regressions to model the probabilities of colonization and of extinction for each species as functions of the habitat suitability estimated by MaxEnt.
Extinction models predict the probability of absence of a species in a particular year, given that it was present before, while colonization models predict the probability of presence of a species, given that it was absent before.
The dependent variable in the colonization models indicates if a BBS route was colonized or not colonized. Colonization occurs when a species was absent (not observed) on a BBS route for k consecutive years and then became present the following year. Non-colonization occurs when a species was absent on a BBS route for k consecutive years and remained absent the following year.
The dependent variable in the extinction models indicates if a species went extinct from a BBS route, or if it persisted. We defined extinction as occurring when a species that was present on a BBS route became absent (not observed) the following year and remained absent for k consecutive years. Persistence occurred when a species that was present on a BBS route is still present in at least one year during the next k consecutive years.
BBS sampling is subject to detection errors: a species is present on a route, but undetected. A false absence is more likely than a false presence (Lobo et al., 2008). Some studies have attempted to . This procedure has several problems, in our view (Appendix S1). Instead of attempting to correct for detectability, we varied the value of k. Using a higher value of k (i.e. a species is absent only if it has not been observed on a given route for k consecutive years) reduces the risk of false absence, but it also reduces sample sizes. When a route was not sampled every year in a given period of k years, occupation status is unknown and the observation was not used. We repeated all analyses with k = 1, 2, 3 and 4 years (Appendix S4 shows the resulting sample sizes). Below, we present the results for k = 2 years unless stated otherwise. Results are similar for k = 1 to k = 3. When k = 4, there are fewer significant relationships between colonization or extinction probabilities and habitat suitability (see the sensitivity analysis in Appendix S5).
The independent variable in colonization and extinction models is habitat suitability for a given species, estimated by MaxEnt models, at each site in the potential range. We separated habitat suitability into two components: 1) previous suitability, that is, suitability in the year(s) before the colonization or extinction event, and 2) ∆suitability, that is, the difference between suitability in the years(s) prior to the colonization or extinction event and the suitability in the year(s) of the event. When k > 1, suitability is averaged over the period of k years. Δsuitability > 0 indicates increased habitat suitability; Δsuitability < 0 indicates decreased suitability.
Most estimated ∆suitability are small, and the error on ∆suitability must be, on average, twice the error on suitability (because ∆suitability combines two measures of suitability). This is a problem since logistic models assume that independent variables are measured without error. To reduce this issue, for each species, we calculated mean Δsuitability for the 10 MaxEnt models produced during the cross-validation procedure. For colonization and extinction analyses, we included a site for a given species only if all 10 MaxEnt models indicated that suitability had changed in the same direction.
This eliminated ≈45% of cases for analyses of colonization versus non-colonization and 50% for extinction versus. persistence.
Non-colonizations are much more frequent than colonizations, and persistences are more frequent than extinctions. Logistic regression models underestimate the probability of rare events (colonizations or extinctions in this case) (Komori et al., 2015;Maalouf & Trafalis, 2011). We therefore resampled the data such that the numbers of colonizations and non-colonizations (or extinctions and persistences) were equal.
We modelled colonization and extinction in relation to previous suitability and ∆suitability. If MaxEnt models capture a causal relationship between species occurrence and environmental variables, then the expected effects of these variables are the same because they measure variation of the same variable-habitat suitability-through space (previous suitability) and through time (∆suitability). This prediction provides a particularly strong test of the causal link.
Next, we tested whether the apparent strength of MaxEnt models (through space) is indicative of their capacity to predict temporal changes in species' distributions. To do this, we determined the correlation between the AUC of MaxEnt models and two measures of the predictive accuracy of the colonization and the extinction w i = a i ∑ 10 1 a i models (AUC and Brier score, estimated with cross-validation) and one measure of goodness of fit (McFadden's R 2 ). As the three measures led to the same qualitative conclusions, we focus on the AUC results below. Bahn and McGill (2007) showed that "many explanatory variables with suitable spatial structure can work well in species distribution models" and that "the predictive power of environmental variables is not necessarily mechanistic" (see also Beale, 2008;Chapman, 2010;Fourcade et al., 2018;Currie et al., 2020). Because habitat suitability predicted by MaxEnt is always spatially autocorrelated, we added a spatial autocovariate (SA) to the colonization and extinction models.
This allows us to ask: How much of the variability in colonization and extinction can be explained by habitat suitability, after accounting for proximity of conspecifics? To answer this question, we partitioned the deviance in the colonization and extinction models related to previous suitability, ∆suitability and SA using standard variance partitioning methodology (Ozdemir, 2015). proportional to the square of the distance from the focal route.

| 30-year change in suitability and occupancy
The analyses above test whether colonization and extinction probabilities, calculated over two to four years, are related to environment suitability or neighbourhood occupancy. Changes in route occupancy may be more strongly related to changes in suitability

| RE SULTS
Extinctions (versus persistence) and colonizations (versus noncolonizations) showed considerable spatial variability across species' ranges (see maps for each species in Appendices S6 and S7). The spatial patterns are qualitatively similar for different values of k (the period used to define extinction or colonization).

| MaxEnt models
By the standard criterion used to judge MaxEnt models (viz., AUC > 0.7 represents acceptable quality: but see Lobo et al., 2008;Raes & Ter Steege, 2007), the models for the 21 species in this study described their geographic distributions well. AUC ranged from 0.72 to 0.90 (median = 0.84; mean = 0.83). It therefore seems reasonable to use these MaxEnt models to calculate habitat suitability for each species, at each site, in each year. Per cent vegetation cover, temperature and elevation were the most important predictors, based on permutation importance (Searcy & Shaffer, 2016).

| Extinction models
Local extinctions occur fairly frequently (median = 408 local extinctions per species over the study period; mean = 470). Persistence is more frequent than extinction, but the ratio of extinctions to persistences varies greatly among species, from 0.027 to 0.472 (median ratio = 0.106, mean ratio = 0.084). For most species, extinctions occurred throughout the species' range, rather than mainly at a range margin (Appendix S6).
The variation in extinction probability among sites is related to habitat suitability. In regressions predicting extinction probability from previous suitability and ∆suitability, the former is always significant, and the latter is significant for most species (15 out of 21). Significant coefficients are always in the predicted direction (Table 1). However, extinctions are also strongly related to the occupancy of neighbouring sites. Adding previous SA (the spatial autocovariate in the year preceding extinction or persistence) to the model reduces the importance of previous suitability, and it causes ∆suitability to become non-significant for most species (Table 1). In models including all three predictors, previous SA is significant for 15 of 21 species (Table 1).
The largest amount of deviance is explained by the collinear variation of previous suitability and previous SA (Figure 1). On average, the amount of deviance uniquely explained by previous suitability is small, and it is similar to that of previous SA (paired t test, two-tailed MaxEnt models make better predictions of extinction/persistence when k (the time window used to define absence) is larger (Appendix S5, Tables E.17 -E.32). This indicates either that longerterm extinctions (a species is absent for ≥ 4 years) are generally better predicted by habitat suitability than short-term extinctions (e.g. a species is absent for 1 year), or that there is significantly less error in the dependent variable.
Variation in colonization probability among sites is related to habitat suitability (Figure 3). In the suitability-only model (i.e. ∆suitability + previous suitability), the effects of ∆suitability and previous suitability on colonization probability are always in the predicted direction when they are significant (Table 2). Previous suitability was significant for all species, whereas ∆suitability had F I G U R E 1 Partitioning of the explained deviance in models of extinction for 21 species of passerine birds. Extinction was modelled as a function of previous habitat suitability (Suit prev ), change in habitat suitability (∆Suit ) and previous spatial autocovariate (SA prev ). For each species, the total explained deviance is divided into the proportion uniquely explained by each predictor and the proportion that is jointly explained by different combinations of predictors. The symbol ∩ indicates the proportion of explained deviance that is shared between predictors. For example, Suit prev ∩ SA prev shows the deviance that is related to the collinear variation of Suit prev and SA prev . The boxplots represent the variation of explained deviance among the 21 species TA B L E 1 Direction and statistical significance of the coefficients of the predictor variables in two different extinction models  (15) Note: The first model includes the effects of ∆suitability and previous suitability. The second model includes the effects of ∆suitability, previous suitability and previous spatial autocovariate (SA).

| 30-year change in suitability and occupancy
Extinction and colonization were, at best, weakly related to longerterm change in environmental suitability. This is, in part, because suitability rarely changed monotonically over 30 years for a species on a given route. Few changes in suitability over 30 years were significantly different from zero ( Figure 5). Changes were equally likely to be positive or negative. Generally speaking, climatic conditions neither systematically improved nor worsened for most of our study species on any given route.
Changes in the logit probability of occupancy over 30 years were significantly related (at α = 0.05) to change in habitat suitability over 30 years for only six of the 21 species ( Figure 6). However, all of the significant changes were in the direction expected if changing suitability induces ranges to shift. Moreover, for 15 of 21 species, the sign of the change in occupancy was in the expected direction of the hypothesized causal relationship (binomial P(≥15 successes in 21 trials)=0.039).

| Does Maxent capture a causal relationship?
Extinction and colonization probabilities both relate to previous suitability in the expected way (i.e. the probability of colonization is higher, and the probability of extinction lower, at sites with higher previous suitability). Extinction and colonization also relate to the presence/absence of conspecifics at neighbouring sites. The direction of causation is unclear: a site may be surrounded by conspecifics because regional habitat is suitable. Or, a poor-quality site surrounded by neighbours may be quickly recolonized and therefore appear to be suitable in an SDM. Because of this collinearity, SDMs provide only a weak test of the hypothesis that species distributions are causally related to "habitat suitability." If the relationship is causal, then changing suitability should lead to changes in occupancy, that is colonization and extinction.
Thus, relationships between ∆suitability and the probabilities of extinction and colonization provide stronger tests of the causal link. Soroye et al. (2020) suggest that, for bumblebees, changes in extinction and colonization are strongly related to interannual fluctuations in temperature and precipitation, which cause sites to flip between being climatically suitable versus. unsuitable. We found that the effect of ∆suitability is very small, often undetectable. In fact, once neighbour effects are accounted for, extinctions and colonization are only rarely detectably related to ∆suitability. can be explained by interannual variation in abiotic conditions. For most of our species, either the causal link is very weak, and/or the effect of ∆suitability over one to four years was too small to be detected.
Integrated over 30 years, changes in occupancy were rarely significantly related to change in suitability, but those changes in occupancy were in the predicted direction for most species. Lack of significance is not due to low statistical power: there was a median of 546 routes per species, sufficient to detect relationships as weak F I G U R E 2 Does habitat suitability, estimated from MaxEnt models, make better predictions of extinctions when the MaxEnt models are better (i.e. have higher AUC)? Here, we show one measure of the strength of the extinction ~ habitat suitability relationship (extinction AUC) as a function of the AUC of the MaxEnt model for 21 species of Passeriformes. The relationships using Brier score or McFadden R 2 instead of AUC were nearly identical to this one as r 2 = 0.007. Instead, despite climate change from 1979 to 2009, net change in habitat suitability for most species was very small. Sitelevel occupancy changed frequently, but < 2% of the variability in occupancy through time was related to change in (mainly) climatedriven suitability.
In this light, are projected changes in species' geographic ranges based on SDM-estimated suitability and projected climate change likely to be accurate? Many studies have argued that species' geographic ranges track changing climate, shifting towards the poles and/or towards higher elevation. Yet, that hypothesis is rarely explicitly tested (Parmesan & Yohe, 2003;Hitch & Leberg, 2007;Virkkala & Lehikoinen, 2014; but see Taheri et al., 2016;Currie & Venne, 2017;Zhu et al., 2012). Our results (and those of Sofaer et al., 2018) suggest that changes in occupancy are largely independent of changes in site suitability. MaxEnt predictions of range shifts due to climate change over a few years or decades (e.g. Matthews et al., 2011;Reese & Skagen, 2017) are unlikely to be accurate, for reasons discussed below.
Might SDM algorithms other than MaxEnt might make better predictions? For some species, perhaps. However, Rapacciuolo et al. (2012) observed that, averaged over many species, the mean differences in prediction accuracy between MaxEnt and several other algorithms, including ensemble models, were small.
Or, might lags between changing environmental suitability and species' responses (Devictor et al., 2008(Devictor et al., , 2012La Sorte & Jetz, 2012) conceal a relationship between the two? Perhaps, but many extinctions and colonizations occurred despite environmental conditions changing very little. Thus, most changes in occupancy F I G U R E 3 Partitioning of the explained deviance in models of colonization for 21 species of Passeriformes. Extinction was modelled as a function of previous habitat suitability (Suit prev ), change in habitat suitability (∆Suit) and previous spatial autocovariate (SA prev ). For each species, the total explained deviance is divided into the proportion uniquely explained by each predictor and the proportion that is jointly explained by different combinations of predictors. The symbol ∩ is not an interaction; rather, it indicates the proportion of explained deviance that is shared between predictors. For example, Suit prev ∩ SA prev shows the deviance that may be explained by the collinear variation of Suit prev and SA prev . The boxplots represent the variation of explained deviance among the 21 species in the observed range-may change through time (function g, above).
In this light, it is not surprising that the common measure of SDM performance, AUC, is related to MaxEnt models' ability to predict extinction ( Figure 2), but it is not related to the models' ability to predict colonization (Figure 4).

| How does Maxent perform compared to a simple spatial model?
Like Eaton et al. (2014), we found that estimates of neighbourhood occupancy (SA) relate positively to colonization probability and negatively to extinction probability. The effect of SA was generally similar to, or greater than, the effect of habitat suitability.
To a degree, the presence of conspecifics at neighbouring sites is probably an indicator of regional habitat suitability. However, the strong partial effect of SA in both colonization and extinction models when controlling for environmental suitability indicates that species distributions exhibit spatial structure beyond that owing to the environmental variables in the MaxEnt model (Phillips et al., 2017).  , 2014). If the effect of SA was due to missing habitat variables, we would expect the effect of SA to be similar in the colonization model and extinction models, or perhaps smaller in the former, because a BBS route that was unoccupied could remain unoccupied if there are no nearby colonists, despite a suitable environment.
In contrast, low suitability should have a more immediate impact on an occupied route. We found the opposite: the effect of SA is usually greater for colonization than extinction (when SA is based on presence). This, and the lack of correlation between the predictive ability of colonization models and MaxEnt models, suggests that colonization is more strongly driven by endogenous processes that are not captured by the SDMs. Clearly, colonization always requires dispersal. Extinction, in contrast, only happens when there is no successful dispersal. We propose that this difference leads to a greater effect of dispersal on colonization than extinction.
Given that birds are among the most mobile organisms, we expect that dispersal limitation will be even more important for other taxa responding to climate change. Thus, when projecting species' responses to climate change, the spatial structure of the species' distributions should be considered. The use of a spatial autocovariate (SA) had been proposed as a way to account for dispersal limitations (Yackulic et al., 2015). Although it is unclear whether SA measures dispersal limitation or some other process, one thing is clear: including the SA generally increases the predictive ability of species distribution models.
In conclusion, we have shown that, MaxEnt models apparently capture a real, but small, effect of environmental conditions on changes in species' distributions. Predictions of extinctions over years to decades are weak, and even weaker for colonizations. Over short time scales, occupancy of neighbouring sites by conspecifics predicts changes in occupancy as well as, or better than, MaxEnt.
Future models of species' responses to climate change should 1) recognize that the processes of colonization and extinction are not equally well predicted by species distribution models and 2) account for the spatial structure of species distribution.

ACK N OWLED G EM ENTS
We thank the Natural Sciences and Engineering Research Council of Canada for funding and the editors and four anonymous reviewers for constructive comments on earlier drafts of this study.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13238

DATA ACCE SS I B I LIT Y
All of the data in this paper were downloaded from publically accessible websites cited in the main text.