The effects of climate change on Australia’s only endemic Pokémon: Measuring bias in species distribution models

Species distribution models (SDMs) are frequently used to predict the effects of climate change on species of conservation concern. Biases inherent in the process of constructing SDMs and transferring them to new climate scenarios may result in undesirable conservation outcomes. We explore these issues and demonstrate new methods to estimate biases induced by the design of SDM studies. We present these methods in the context of estimating the effects of climate change on Australia's only endemic Pokémon. Using a citizen science dataset, we build species distribution models for Garura kangaskhani to predict the effects of climate change on the suitability of habitat for the species. We demonstrate a novel Monte Carlo procedure for estimating the biases implicit in a given study design, and compare the results seen for Pokémon to those seen from our Monte Carlo tests as well as previous studies in the same region using both simulated and real data. Our models suggest that climate change will impact the suitability of habitat for G. kangaskhani, which may compound the effects of threats such as habitat loss and their use in blood sport. However, we also find that using SDMs to estimate the effects of climate change can be accompanied by biases so strong that the data themselves have minimal impact on modelling outcomes. We show that the direction and magnitude of bias in estimates of climate change impacts are affected by every aspect of the modelling process, and suggest that bias estimates should be included in future studies of this type. Given the widespread use of SDMs, systemic biases could have substantial financial and opportunity costs. By demonstrating these biases and presenting a novel statistical tool to estimate them, we hope to provide a more secure future for G. kangaskhani and the rest of the world's biodiversity.


| INTRODUC TI ON
Anthropogenic climate change has been demonstrated to negatively impact the distribution of suitable habitat for many species, leading to local extirpations, range shifts and extinctions (Bellard et al., 2012;Chen et al., 2011;Parmesan, 2006). Managers attempting to mitigate the effects of climate change often have to make decisions that require some information about the tolerance of species for given combinations of temperature, precipitation or other environmental variables in order to decide whether a given area is currently suitable, or will be so given expected environmental change (Hutchinson, 1978;Soberon & Nakamura, 2009). However, experimental approaches to estimating this 'environmental niche' are often not feasible for reasons of cost or due to other practical concerns. As such, a number of correlative approaches have been developed that attempt to derive estimates of species' tolerances by correlating species occurrence data with the geographic distribution of environments. These species distribution models (SDMs, also known as environmental niche models or ENMs) are often one of the key lines of evidence used to make decisions to mitigate the effects of climate change. Unfortunately there are many methodological and conceptual issues involved in the construction and application of SDMs that remain poorly understood.
One major concern in SDM studies is that the methodological choices made in the SDM process may bias the predictions that the resulting models make (Barbet-Massin et al., 2010;Beaumont et al., 2016;Datta et al., 2020). Using SDMs requires decisions to be made regarding the choice of algorithm, predictor variables, emissions scenario, climate model and study area among other parameters. In turn, these choices may affect the distribution of environments available for modelling as well as the distribution of environments that models will be projected onto, and as such each of these aspects of an SDM study's experimental design has the potential to introduce bias. Clearly, more approaches directed towards quantifying and mitigating biases in SDM studies are needed.
In order to better understand these issues, we use SDMs to examine the effects of climate change on Australia's only endemic Pokémon: kangaskhan (Garura kangaskhani, Figure 1). To estimate the extent to which our results might be driven by methodological biases, we compare our predictions of changing habitat suitability for G. kangaskhani to similar models built from virtual species (Warren et al., 2020), and to a recent study of 220 species of Australian mammals (Beaumont et al., 2016).
Garura kangaskhani is the only Pokémon endemic to the geographic region used for these two studies, and is therefore the Pokémon most comparable to those studies when examining the effects of SDM study design. We also perform a post hoc analysis on models from Warren et al. (2020) to illuminate which aspects of SDM study design are responsible for inducing biases when predicting the effects of climate change. Finally, we develop a novel Monte Carlo method that can be used to estimate the bias in any given SDM study design, and compare our predictions for G. kangaskhani to those expected by bias alone.

| Study system
We provide a detailed discussion of the current (woefully inadequate) state of knowledge about G. kangaskhani in supplement S1, but for the sake of brevity we restrict ourselves here to those aspects of its biology most relevant to conservation. Garura kangaskhani is an omnivorous (Bulbapedia, 2020a), diurnal (PokéBase, 2020) vertebrate endemic to Australia. It is primarily associated with grasslands and savanna habitats, although frequent sightings by field researchers in highly disturbed metropolitan areas indicate that they are also fairly resilient to anthropogenic disturbance; one study of 'Normal' type Pokémon (which includes G. kangaskhani) found that parking lots and universities were the second and third most common sources of sightings, with sightings more frequent during partly cloudy weather

| Threats
The vulnerability of animal and human populations is often a complex function of multiple interacting stressors (O'Brien et al., 2004).
Garura kangaskhani were previously hunted to the brink of extinction (Bulbapedia, 2020b), and one of the primary threats to the persistence of G. kangaskhani populations continues to be the poaching of adults and eggs. While kangaskhan may sometimes be eaten (Gilbert, 2020), poaching is primarily motivated by the acquisition of specimens intended for blood sport; captured animals are used in contests in which the maternal defence of the young is exploited to motivate adult G. kangaskhani into fighting both conspecifics and heterospecifics (Molloy, 2013). Individuals deemed to be unfit for combat are frequently ground into 'candy' which is fed to combatants (D'Anastasio, 2016). This destructive practice is not limited to G. kangaskhani, nor is it geographically restricted to Australia. Current estimates indicate that tens of millions of 'trainers' and 'breeders' are participating in these matches worldwide (Wikipedia, 2020), with hundreds of millions of dollars being spent on supplies alone (Phillips, 2018). With the recent discovery of a rare 'shiny' morphotype of G. kangaskhani that can add a visual enhancement to a trainer's 'battle league party', we anticipate that collection pressure on this species will only increase in the future. We see little room for optimism in this scenario given that the culture around Pokémon harvesting so explicitly favours overexploitation, as seen in the common refrain of 'gotta catch 'em all' (Prof. Y. Ohkido, pers. comm.).
Compounding the threats already posed by exploitation, anthropogenic climate change may constitute an additional challenge for G.

kangaskhani.
Little is known about the species' climate niche beyond a marked preference for partly cloudy weather (PokemonGoHub, 2020;Pokémon GO Wiki Community, n.d.). It is therefore difficult to estimate a priori the response of G. kangaskhani to environmental change, or how this additional stressor might combine with the effects of poaching to impact the species' long-term survival. Moreover, their apparently broad habitat tolerances and the fact that most weather can be considered at least 'partly' cloudy, suggest that G. kangaskhani has few strong climate associations within Australia. This lack of strong associations may challenge our ability to make precise forecasts for this species, but may serve to make this species more useful for understanding the effects of bias; if a method shows a tendency to make strong predictions when there is no actual relationship between environmental predictors and the occurrence of the species, it suggests that those predictions may be driven by methodological biases.
In this study, we use a citizen science dataset of occurrences to estimate the effects of climate change on populations of G. kangaskhani. These data were initially recorded by hobbyists and professional trainers seeking out G. kangaskhani for exploitation, and as such the use of these data for conservation modelling does pose some ethical issues. However, given the uncertain conservation status of these noble creatures and the paucity of georeferenced localities of pokémon in natural history collection databases, few other options exist. We also acknowledge that the small sample size used here (37 points) may be cause for some concern, but note that this is not an unusually small sample size for SDM studies (Galante et al., 2018;Loe et al., 2012). Further, we compare our results to other SDM studies with N = 100 (Warren et al., 2020) and N ranging from 10 to 7,137 (Beaumont et al., 2016), and find broadly concordant patterns (see Section 4). We also note that seminal contributions have been made by investigators using data sources of similar quality (Hurlbert, 1990;Lozier et al., 2009;Michelot et al., 2016), and as such we feel that the use of these data are justified.

| MATERIAL S AND ME THODS
We obtained occurrence data for G. kangaskhani from online databases maintained by individuals seeking to collect them for use in blood sport or for personal gains that exploit the lack of international trade regulations for this species (PKMNGoTrading, 2016;KOPlayer, 2018;The Road and Community, 2016). Given how little is known about G. kangaskhani, many decisions in the modelling process were necessarily made based on common practice in the literature rather than a priori knowledge of the underlying biology (Araújo et al., 2019). In order to reduce the impact of spatial sampling bias we removed duplicate occurrences, as is common in SDM studies (Peterson et al., 2011). We also removed any occurrences falling in regions with no environmental data, resulting in 37 unique data points. Climate data were obtained from CliMond (Kriticos et al., 2012) at 10-km resolution. We used the 19 Bioclim layers for both present and future projections. Using the ENMTools r package (Warren, Matzke, et al., 2021), we constructed SDMs using six different algorithms; generalized additive models (GAM), generalized linear models (GLM), maximum entropy (Phillips et al., 2006), random forests, Bioclim (Nix & Busby, 1986) and Domain (Carpenter et al., 1993). As this study is based on a citizen science dataset and limited research funds are available for systematic sampling of Pokemon distributions, true absence data are not available for G. kangaskhani. We therefore adopt a presence/background modelling approach that is common in SDM studies, in which environmental conditions at occurrence points are contrasted with the distribution of available environments within the species range (Peterson et al., 2011). This approach is common in SDM studies without true absence data (Araújo et al., 2019).
Background points were extracted from a 250-km circular buffer around each occurrence. Models were built using 10,000 background points. For each model, 30% of the presence data were randomly withheld for model validation. Code for model construction is given in Appendix S3 .
In order to summarize the effects of climate change on the distribution of suitable habitat, we projected models onto each future climate scenario and then calculated the difference between the current and future suitability at each grid cell in the study area. We then average those projected changes to get the overall mean change in habitat suitability projected by a given combination of distribution model, emissions scenario and climate model. Change in habitat suitability was evaluated at two spatial scales; within the species' current range, and at a continental scale including all of Australia and Tasmania.
To evaluate the extent to which predictions concerning the future of G. kangaskhani may have been driven by methodological bias, we compare our empirical predictions to predictions derived from two other sources. First, we examine the predicted effects of climate change on a set of models for 220 simulated species in Australia from Warren et al. (2020). Second, we compare our observations for G. kangaskhani to a distribution of predictions from a novel Monte Carlo approach that randomly chooses occurrence points from the study area.

| Post hoc analysis of simulations
In the simulations from Warren et al. (2020), virtual species' ranges were derived from the availability of suitable habitat based on an underlying simulated niche, with a variety of niche breadths and range sizes. Within these ranges, species occurrences were sampled in proportion to the product of the suitability of habitat and F I G U R E 2 Graphical sketch of Monte Carlo method for estimating bias in model construction and transfer. For simplicity we demonstrate these effects using a simple bioclim model, but a similar argument holds for any algorithm. In panel (a), we choose our study area (dashed lines) based on the distribution of our occurrence points (Gangura kangaskhani silhouettes). This affects the distribution of environments available during model construction and evaluation (dashed area, panel b). Using the occurrence points, we build a model that parameterizes the distribution of the species' occurrences in environment space (dark rectangle, panels b and c). Given an emissions scenario and climate model, the distribution of available environments may differ in a future time period (dashed area, panel c). Hence, almost any bioclim-style model built using points from the available environment in panel (b) will estimate as suitable some set of environmental conditions that are no longer present in panel (c).
Similarly, models will lack information to determine the suitability of the environmental conditions available in panel (c) that were not present in panel (b). Climate envelope methods such as bioclim do not tend to extrapolate much beyond the conditions they were trained in, so the change in available environments creates a bias towards predictions of declining habitat suitability and/or range contractions. Other modelling algorithms may have similar built-in biases, but the direction and magnitude may differ with their individual tendencies to extrapolate. It is worth noting that the bias inherent in a given study design will be affected by the change in available environments (which is affected by choice of study area, emissions scenario and climate model) as well as the shape and extent of the region of environment space predicted to be suitable (which are affected by study area, sample size and modelling algorithm). As such, all of these choices may affect our baseline expectation for what models will predict in a given study design independent of the actual locations of the species occurrences. To measure this bias, we repeatedly permute the locations of our data points within the study area and construct models using these random data (panel d), with every other aspect of the modelling process identical to the study design used for our empirical data. This eliminates any biological connection between the occurrence of our species and the environmental gradients, allowing us to estimate a distribution (panel e) of expected changes under the hypothesis that our occurrence data contain no information, and our predictions are consequently driven entirely by the bias implicit in the other choices made in our study design (sample size, study area, emissions scenario, climate model, etc.) a model of spatial sampling bias of varying strength (see Warren et al., 2020 for details). Models from these datasets were constructed using methods similar to those used to analyse G. kangaskhani, and were projected onto the same set of future scenarios. In order to understand the sources of bias and uncertainty in these predictions we also conducted a set of post hoc GLM analyses in which the average suitability change and associated errors from these models were regressed against year, with interaction terms for SDM algorithm, emissions scenario and climate model. Significant interactions from these regressions indicate that a given component of the study design modifies the direction or magnitude of the bias in predictions of future habitat suitability.
See Appendix S2  for details and R code for these analyses.

| Monte Carlo estimates of methodological bias
To estimate biases in our predictions for kangaskhan, we developed a novel Monte Carlo approach that is derived from a test developed by Raes and ter Steege (2007). As originally designed, this approach is used to generate null distributions of predictive performance for species distribution models. The null distribution is constructed by selecting localities at random from within the study area until a dataset of the same sample size as the empirical dataset is obtained.
These random localities are then used to generate SDMs using the same methods as those used for the empirical model, and the predic- where N is the number of occurrence points used for modelling in the empirical dataset. Using these random points as 'occurrences' we build SDMs with all other modelling choices being identical to those for the empirical dataset. We then project our replicate models onto the different combinations of climate model and emissions scenario, and calculate the predicted changes in habitat suitability for those models in the same way as with the empirical data. This allows us to estimate a distribution of expected changes in habitat suitability given the same study area, SDM algorithm, emissions scenario and climate model as our empirical data, under the hypothesis that the actual occurrence of the species within the study area is unrelated to the environmental variables used for model construction.

| Interpretation
Our post hoc analyses of simulations from Warren et al. (2020) is intended to quantify the expected behaviour of our modelling approaches given a dataset in which there is an underlying (simulated) biological mechanism relating its probability of occurrence to a set of environmental predictors, while the novel Monte Carlo approach is intended to estimate the behaviour of a model built under identical conditions to our empirical model (including study region), but for which there is no biological mechanism underlying the distribution of occurrence points within that study area. The distributions of predictions produced by these two approaches are not necessarily intended to represent a null hypothesis that must be rejected in order to trust model results, but rather as estimates of the behaviour of our modelling approaches independent of our data. By comparing our empirical results to these distributions, we can determine to what extent our predictions are driven by our data, over and above the biases induced by our analytical framework.

| RE SULTS
We find that the predicted change in suitability of habitat for G. kangaskhani is relatively consistent across spatial scales (continental vs. within species' current range), but depends strongly on the method of analysis. Models built using Bioclim, Domain or GAM algorithms predict that habitat suitability will decline on average, but disagree substantially on the magnitude of that decline (  Year Predicted change within current range

| D ISCUSS I ON
In this study, we used a citizen science dataset in conjunction with climate data to estimate the effects of climate change on Australia's only endemic Pokémon, and find that our conclusions about changes in habitat suitability are profoundly affected by methodological choices. By examining multiple simulated data sources, we highlight the important possibility that the results seen here might not be indicative of any direct effect of climate change on the suitability of habitat for G. kangaskhani, but instead represent the effects of bias that is imposed by the methods used to construct these models. The predicted effects of climate change on the suitability of habitat for G. kangaskhani generally reflect predictions we derive from models for simulated species as well as predictions from Monte Carlo tests based on randomly chosen data points within the range of G.

kangaskhani.
It is interesting also to compare our results with Beaumont et al. (2016), a study that used models built with real data for 220 Australian mammal species to examine which modelling approaches were more prone to making predictions of range contractions or expansions ( Figure 6). Those models were based on real data for species for which the climate niche is not known, and as a result it was not possible from that study to say which of those models better estimates the true impacts of climate change. However, we note that many of the model behaviours seen in Beaumont et al. (2016) (Figures 4-6). The similarity between the results for real species, simulated species and Monte Carlo replicates strongly suggests that the effects of the biases inherent in the design of SDM analyses are a ubiquitous feature of studies of this type.
In this study, our Monte Carlo tests focus on how a given set of modelling conditions may bias our estimates of the effects of climate change on average habitat suitability. However, we note that this framework can easily be extended. There are many different metrics that may be used to quantify the effects of climate change on species (e.g. change in habitable area based on thresholded suitability scores, geographic shifts in ranges or centroids, changes in climate niche breadth, etc.). Applying this approach to those metrics is simply a matter of applying the chosen metric to predictions from each Monte Carlo replicate to obtain the distribution of that metric given randomly sampled occurrence data. Similarly, this approach could easily be extended to obtain distributions for models being trans- Year Predicted change within current range design more than by the data used in modelling (Figure 7). Finally, we note that in the current study occurrence points for Monte Carlo replicates were sampled with uniform probability across the study with no spatial autocorrelation or sampling bias. While this approach has been widely used for constructing null distributions of model performance (Bohl et al., 2019;Raes & ter Steege, 2007;Warren, Dornburg, et al., 2021), other methods that attempt to maintain similar spatial autocorrelation to empirical data when sampling may yield different null distributions (Beale et al., 2008;Nunes & Pearson, 2017;Roxburgh & Chesson, 1998). We consider this a very promising area for future study, as it may help to better understand which aspects of data collection and curation are most instrumental in reducing bias in predictions.
Species distribution models are widely used with the goal of guiding management and policy. However, the fact that we produce consistent predictions from meaningless data indicates that there are some combinations of emissions scenario, climate model, study area and species distribution modelling method that are strongly biased towards a given prediction. Crucially, in some modelling conditions we find that all Monte Carlo replicates agreed on the direction of change in habitat suitability in the future, despite the fact that the points for these replicates were chosen at random from within Kangaskahan's range and therefore represented no underlying relationship between the environment and the probability of occurrence. At the very least this suggests that those study designs should be chosen only when there is a compelling reason to believe that they accurately represent the underlying biology. In the case of G. kangaskhani we have insufficient information to make such a decision; in addition to our scant knowledge of the species' ecology we anticipate that there are likely to be serious issues with data quality due to spatial sampling bias and small sample size due to an unfortunate lack of scientific attention to wild Pokémon populations. In particular, we would like to point out that there may be some room for optimism given our level of ignorance of the distribution of G. kangaskhani; we do not currently know how far into the interior of Australia these animals might be F I G U R E 6 Projected proportional range change for 220 Australian mammal species under two different climate change scenarios, modified from Beaumont et al. (2016) Figure 3. In that study models were built using fourteen different algorithms, but here the figure is modified to include only the algorithms used for Gangura kangaskhani to facilitate comparison with simulation and Monte Carlo results. Proportional range change is calculated as the future area of suitable habitat divided by the present area of suitable habitat based on binary predictions made from thresholded models. As such it is not directly comparable to the change in average suitability presented here for the other analyses, but the two are expected to be highly correlated (i.e. declining suitability scores will typically lead to a predicted reduction in the area of suitable habitat) Proportional range change F I G U R E 7 Using the 100 Monte Carlo replicates for the Maxent models for Gangura kangaskhani, we recorded the mean change in suitability (panel a) and the standard deviations of those predictions (panel b) when predicted to all four combinations of emissions scenario and climate model for the year 2100. As these models represent the behaviour of the study design given uninformative occurrence data, the desired behaviour is for the mean predicted change (panel a) to be near zero (indicating no directional bias), and the standard deviation around that mean (panel b) to be high (indicating that model predictions are sensitive to changes in the data). Although there are some regions that approximately fit this set of desired outcomes, we note that there are areas along the northern and southern coasts where predicted changes are largely insensitive to the data (SD is low, panel b), and these areas include predictions of both increased (northern coast and Tasmania) and decreased (southwestern coastal areas) habitat suitability from the present day (panel a). These indicate regions where Maxent models are biased towards predicting a specific direction of change in habitat suitability, and where the data have little leverage with which to counteract these biases found and so are unlikely to have fully characterized its environmental tolerances. Currently very little conservation funding is allocated for land acquisition and captive breeding programs for Pokémon, but increased investment in this area could be broadly beneficial; charismatic megavertebrates can often be used as flagship species to increase community involvement in conserving whole ecosystems (Cummins, 2020). In order to effectively manage this majestic species, we suggest that future funding efforts be directed towards sending the authors of this study into remote areas of Australia to collect more reliable data, preferably equipped with a large supply of the 'incense' (presumably some sort of pheromone) that Pokémon hunters use to attract animals for capture.
It is reasonable to view with heavy scepticism the results of any study in which the analytical design is only capable of producing one answer regardless of the data. This may not mean, however, that the predictions from such a study are necessarily wrong. To illustrate this point, it is worth considering an extreme case in which the entire study area becomes unsuitable for any organism currently living within it (e.g. in future scenarios all currently occupied habitat is constantly on fire). In this case, any reasonable model built from any set of random occurrence data in that region would show complete extirpation of the organism, which would also coincide with the results seen from any carefully vetted empirical dataset for a real species. The conclusions of those models would be correct, but they would be wholly unaffected by the data; they would be an inevitable result of the modelling conditions.
Any dataset would likely support the same conclusions regardless of whether those data represent any underlying biological phenomenon.
In the real world, however, the effects of climate change are rarely so absolute, and are rarely known in advance. Finding that we are strongly biased towards making specific predictions irrespective of our data is therefore a cause for great concern; it implies that we may have essentially selected our conclusions by selecting our study design. In such a situation, the resulting outputs should at minimum be interpreted in the context of the biological plausibility of those methodological choices, the biases they induce, and their consequences for decisions made based on those models. In particular, if we find ourselves in a situation where different approaches show strong opposing biases, interpretation of those results must be made in the context of how biologically realistic the alternative modelling approaches are. Arguably one should go further than that, however, and simply not trust the outcome of a modelling approach that produces the same result regardless of the data provided to it. In these cases, additional lines of evidence are needed to adequately inform conservation decisions. Understanding when we are in such a situation is the key to making these decisions, and we believe that Monte Carlo tests such as the one presented here may be a valuable tool for measuring the biases inherent in the design of empirical studies.

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