Species distributions models may predict accurately future distributions but poorly how distributions change : A critical perspective on model validation

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Species distributions models may predict accurately future distributions but poorly how distributions change : A critical perspective on model validation © 2023 The Authors. Diversity and Distributions published by John Wiley & Sons Ltd. Published version Piirainen, Sirke; Lehikoinen, Aleksi; Husby, Magne; Kålås, John Atle; Lindström, Åke; Ovaskainen, Otso Piirainen, S., Lehikoinen, A., Husby, M., Kålås, J. A., Lindström, Å., & Ovaskainen, O. (2023). Species distributions models may predict accurately future distributions but poorly how distributions change : A critical perspective on model validation. Diversity and Distributions, 29(5), 654-665. https://doi.org/10.1111/ddi.13687 2023 Diversity and Distributions. 2023;00:1–12. | 1wileyonlinelibrary.com/journal/ddi Received: 4 May 2022 | Revised: 16 December 2022 | Accepted: 8 February 2023 DOI: 10.1111/ddi.13687 RESEARCH ARTICLE Species distributions models may predict accurately future distributions but poorly how distributions change: A critical perspective on model validation Sirke Piirainen1,2 | Aleksi Lehikoinen1 | Magne Husby3,4 | John Atle Kålås5 | Åke Lindström6 | Otso Ovaskainen7,8,9 1Zoological Unit, Finnish Museum of Natural History, University of Helsinki, Helsinki, Finland 2Arctic Centre, University of Lapland, Rovaniemi, Finland 3Section of Science, Nord University, Levanger, Norway 4BirdLife Norway, Trondheim, Norway 5Norwegian Institute for Nature Research, Trondheim, Norway 6Department of Biology, Lund University, Lund, Sweden 7Department of Biological and Environmental Science, University of Jyväskylä, Jyväskylä, Finland 8Organismal and Evolutionary Biology Research Programme, University of Helsinki, Helsinki, Finland 9Department of Biology, Centre for Biodiversity Dynamics, Norwegian University of Science and Technology, Trondheim, Norway This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2023 The Authors. Diversity and Distributions published by John Wiley & Sons Ltd. Correspondence Sirke Piirainen, Zoological Unit, Finnish Museum of Natural History, University of Helsinki, P. O. Box 17, Helsinki 00014, Finland. Email: [email protected] Funding information Academy of Finland, Grant/Award Number: 275606, 309581 and 326338; BiodivScen ERANetCOFUND programme; H2020 European Research Council, Grant/Award Number: grant agreement No 856506; ERCsynergy project LIF; Jane ja Aatos Erkon Säätiö; Koneen Säätiö, Grant/Award Number: 201903886; National Science Foundation, Grant/Award Number: CLO, ICER1927646; Norges Forskningsråd, Grant/Award Number: Centres of Excellence Funding Scheme (223257) and NINA: 295767; Suomen Kulttuurirahasto; Svenska Forskningsrådet Formas, Grant/ Award Number: 201802441 Editor: Boris Leroy Abstract Aim: Species distribution models (SDMs) are widely used to make predictions on how species distributions may change as a response to climatic change. To assess the reliability of those predictions, they need to be critically validated with respect to what they are used for. While ecologists are typically interested in how and where distributions will change, we argue that SDMs have seldom been evaluated in terms of their capacity to predict such change. Instead, typical retrospective validation methods estimate model's ability to predict to only one static time in future. Here, we apply two validation methods, one that predicts and evaluates a static pattern, while the other measures change and compare their estimates of predictive performance. Location: Fennoscandia. Methods: We applied a joint SDM to model the distributions of 120 bird species in four model validation settings. We trained models with a dataset from 1975 to 1999 and predicted species' future occurrence and abundance in two ways: for one static time period (2013– 2016, ‘static validation’) and for a change between two time periods (difference between 1996– 1999 and 2013– 2016, ‘change validation’). We then measured predictive performance using correlation between predicted and observed values. We also related predictive performance to species traits. 14724642, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ddi.13687 by University Of Jyväskylä Library, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2 | PIIRAINEN et al. 1 | INTRODUCTION Understanding how species will respond to climate change is one of the current key challenges in ecology and nature conservation (Malhi et al., 2020; Sofaer et al., 2018). To gain such understanding, various types of species distribution models (hereafter SDMs) are widely used to make predictions on where the species occurs now, and projections on how species distributions are likely to change as a response to predicted climatic change (Huntley et al., 2007; Leathwick et al., 1996; Massimino et al., 2017). Thus, SDMs are important tools for ecologists and can have large implications for management and conservation (Kujala et al., 2013). SDMs can for example help to identify conservation needs, to define alternative conservation actions and to evaluate the effects of such actions (Guisan et al., 2013). The utility of predictions clearly depends on how much they can be trusted, and to estimate reliability the predictions need to be critically validated with respect to what they are used for. In this paper, we seek to simulate a situation where an SDM is intended to be used in projecting community changes to a future time period approximately 20 years from the current moment. As Vaughan and Ormerod (2005) claim, the ideal way to measure a model's generalizability is to test the model by incorporating situations that resemble its future applications as closely as possible. We take this idea as read that if a model is intended to predict, for example 20 years ahead in future, then its validation should reflect this pattern by using data from the past to make predictions into present (a ‘future’ time period in relation to model fitting data), which lies approximately a similar length of time ahead. This means ignoring information from the latest part of the data for model fitting, but simulates the reallife data gap that an ecologist making projections into a more distant future is always faced with. Here, we argue that while ecologists are typically interested in how species distributions will change in a changing climate, SDMs have seldom been critically evaluated in terms of their capacity to predict such change. This has led to a situation where our ability to make reliable range change forecasts may often be overestimated. Even though the benefits of using temporally independent data for testing SDM's predictive performance have been documented already decades ago (Araújo et al., 2005), still a typical SDM used for forecasting is not validated by means of such data (UribeRivera et al., 2022). Typically, in those cases when temporally independent data exist and there has been an attempt to validate the model by means of retrospective forecasting, this has been done by splitting the dataset into two temporal blocks: the earliest records of the dataset form a training data block and the latest a test data block. There might or might not be a temporal gap between the training and testing datasets. Depending on the length of the gap, the test data can be seen as temporally independent data with respect to the training data. The model is then calibrated with training data and predictions are made to one static point in time (either to data from 1 year only or to a mean of multiple years). Lastly, these predictions are compared with observations (‘static validation method’ in Figure 1), and the predictive performance of the model is considered high if the predictions match with the observations. Several authors have used temporally independent static validation, for example on butterflies (Eskildsen et al., 2013), plants (Dobrowski et al., 2011) and birds (Araújo et al., 2005; Regos et al., 2019). However, we argue that if the aim was to evaluate how well a model can predict changes in the patterns of occupancy or abundance of the species, then using a static validation method that only considers the spatial distribution of occupancy or abundance (predicts and validates a single static point in time) is not adequate. While the static method is able to evaluate how well the model predicts a future distribution on a general level, it does not evaluate whether or not the model is able to predict if and how the species distribution changes. This is because the species ranges may remain mostly stable and thus there can be a high correlation between locations where the species has occurred in the past and where it will occur in the near future (Rapacciuolo et al., 2014). Results: Even though static validation method evaluated predictive performance as good, change method indicated very poor performance. Predictive performance was not strongly related to any trait. Main Conclusions: Static validation method might overestimate predictive performance by not revealing the model's inability to predict change events. If species' distributions remain mostly stable, then even an unfit model can predict the near future well due to temporal autocorrelation. We urge caution when working with forecasts of changes in spatial patterns of species occupancy or abundance, even for SDMs that are based on time series datasets unless they are critically validated for forecasting such change. KEYWORDS birds, climate change, Fennoscandia, forecasting, land use, model validation, prediction, species distribution modelling, species traits, temporal transferability 14724642, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ddi.13687 by University Of Jyväskylä Library, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 3 PIIRAINEN et al. So even if the model fails to predict change events, which usually happen especially at the range margins, it may still predict a major part of the unchanged distribution correctly and thus achieve an overoptimistic evaluation of predictive performance when only a static validation method is applied (Rapacciuolo et al., 2012; Sofaer et al., 2018). As a solution on how to effectively test the performance of SDMs intended for predicting biodiversity changes in future, we propose that the model validation method should be explicitly based on measuring and validating a change over a given time period. By ‘change’, we mean the difference in a population measure between two points in time, that is increase, decrease or stability in occurrence probability or abundance. Ideally, the time period over which a change is calculated should be independent in relation to the model training data (‘change validation method’ in Figure 1). We base our change validation method on the work of Rapacciuolo et al. (2014) where the authors developed a temporal validation plot to visualize and measure the agreement between predicted and observed changes and applied it to two bird species from the Great Britain and a virtual species. The temporal validation plots solved the issue of generally static distributions by focussing the assessment of model performance only on locations where range change events (a grid cell either gaining or losing occupancy) over time were either observed or predicted. By focussing on change events only, the authors were able to reveal aspects of the relationship between species' range change and model variables that might not have been identified through rangewide measures (including also grid cells with stable presences or absences). While the need to validate the predicted change rather than validating the predicted future may sound selfevident, to our knowledge our work is one of only a few papers to apply such an approach in the context of climate change predictions (Johnston et al., 2013; Rapacciuolo et al., 2014). A greater proportion of studies have investigated models' ability to predict to new geographical areas by applying spatial validation methods (Bahn & McGill, 2013; Charney et al., 2021; Journé et al., 2020; Rousseau & Betts, 2022). While there is a certain analogy between spatial and temporal block crossvalidation through the spacefortime substitution in distribution modelling (Blois et al., 2013), we, however, think that spatial validation methods are not best suited for testing a model, which is not intended to predict to new areas, but to a new time, and therefore, results obtained from spatial validation studies do not directly apply FIGURE 1 Schematic picture of the two validation methods compared in the study. The evaluated model uses data from 1975 to 1999 for model training (indicated by the grey bar). Static validation method uses the model to predict a static pattern, species occurrence or abundance in only one, independent time period t2 (mean of 2013– 2016, indicated by a circle P2). Predictive performance is calculated as the correlation between the prediction and the corresponding observed value (diamond O2) over study locations. Change method predicts both time periods t1 (1996– 1999, circle P1) and t2 (circle P2) and then counts the predicted change (ΔP), that is difference in species occurrence or abundance between these two separate time periods, and compares it to the corresponding observed change (ΔO). Δ Δ Δ 14724642, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ddi.13687 by University Of Jyväskylä Library, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4 | PIIRAINEN et al. in the context of temporal validation and making predictions to future time. In multispecies studies, there is usually a high level of variability in predictive performance among the species (Araújo et al., 2005; Dobrowski et al., 2011; Sofaer et al., 2018; Venne & Currie, 2021), meaning that some species are easier to model than others. One reason for this variation can be ecological differences between species (Tessarolo et al., 2021). For instance, species traits such as range size, migratory behaviour, rarity and body size have been linked with model performance in birds (McPherson & Jetz, 2007). Studying the connection between species' traits and predictive performance can give us an understanding on how those traits affect modelling and how to interpret results. In this study, we compared our ‘change validation method’ with the more routinely used ‘static validation method’ using largescale monitoring data on 120 bird species from Finland, Sweden and Norway. We studied whether predictive performance differs between these two validation methods and whether the variation in predictive performance is linked with certain species' traits. Concerning the influence of species traits, we expected that common, largesized, nontropical migrant species which prefer naturally patchy habitats show higher predictive performance than rare, smallsized, longdistance migratory species, which prefer broad and structurally more varied habitats such as forest, or are more influenced by factors outside their breeding range (Johnston et al., 2014; Laaksonen & Lehikoinen, 2013; McPherson & Jetz, 2007). 2 | METHODS 2.1 | Bird data We gathered data from national bird monitoring surveys performed in North Europe: Finland, Sweden and Norway (Lehikoinen et al., 2014; Lindström et al., 2015) covering circa 1900 km latitudinal gradient. In total the data consist of 2591 bird survey sites (Norway 486, Sweden 1197 and Finland 908, Figure 2). Data come from line transects (Finland and Sweden) and/or point count routes (Norway and Sweden), which cover the countries systematically at approximately 20– 40 km distance of each other, as well as from more randomly located surveys. Survey sites have been visited once per year during the start of the breeding season, but not necessarily each year. All the data have been gathered by volunteers during 1975– 2016. The earliest observations are from 1975, 1978 and 1991 from Sweden, Finland and Norway, respectively. In each survey, abundances of observed birds were recorded either as pairs (Finland and Norway) or as individuals (Sweden). We approximated that these measures correspond because in the Finnish and Norwegian survey FIGURE 2 Bird survey sites in Norway, Sweden and Finland. (a) Survey sites which were visited between 1975 and 1999 and thus contribute to model training data. (b) Survey sites which were visited at least twice in each of the time periods t1 (1996– 1999) and t2 (2013– 2016) (which contribute to change validation method for occurrence) are indicated with filled triangles (n = 264). Survey sites which were visited at least twice in time period t2 (which contribute to static validation method for occurrence) are indicated with both filled triangles and open circles (n = 1214). (a) (b) 14724642, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ddi.13687 by University Of Jyväskylä Library, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 5 PIIRAINEN et al. methods, observations of solitary males (the by far most common records) and females are in any case, interpreted as pairs. In total, 264 bird species were recorded. In Finland, waterbirds were included in the surveys from 2006 onwards. For each survey, a method (point or line) and effort (number of points belonging to a point count survey or length of line in metres for a line transect survey) of survey was recorded, as well as the starting point or centroid coordinates of the survey site. Traits of each bird species (prevalence, body mass, migratory behaviour and habitat preference) are listed in Table S1.2: Appendix S1. Source data for species traits and more precise hypotheses concerning species traits are presented in Table S1.1: Appendix S1 2.2 | Environmental variables To describe the habitat where birds were observed, we drew a 300mwide buffer around each line transect and group of points in point count route and determined the land cover inside the buffer from the 100m resolution Corine land cover data (© European Union, Copernicus Land Monitoring Service 2022, European Environment Agency, EEA). In our study area, the Corine land cover survey had been conducted in the years 2000, 2006 and 2012. We extrapolated information from this data to cover the years 1975– 2003, 2004– 2009 and 2010– 2016, respectively. We acknowledge that in the earliest period there is a maximum of 25years of mismatch between land cover and bird data. As the long time series from the 1970s provides a strong calibration dataset for the future predictions, and as the land cover variables have changed less than the climatic variables that are of our primary interest, we accepted this mismatch. We classified land cover types into six classes and calculated their proportions (value range 0– 1) on the buffer. The classes were urban and agricultural areas (Urb), broadleaved and mixed forests (Br), coniferous forests (Co), mountainous areas and shrublands (Op), marine water and beaches (Ma) and wetland, peatland and inland water areas (We). For climatic variables, we chose only mean daily temperature measurements, since previous studies have shown that temperature has much larger impact on North European bird populations compared with, for example rainfall (Tayleur et al., 2015). For each survey site, we determined mean temperatures of the current year spring (April and May), previous winter (December– February) and previous summer (June and July of the preceding year) from EOBS weather data (version 15, 0.25degree resolution) by the European Climate Assessment & Dataset project, ECA&D (Cornes et al., 2018). We selected these weather parameters because they reflect the conditions which the birds and their breeding environment have experienced by the time of the observation, after the past summer and winter and approximately at the start of the breeding season (PearceHiggins et al., 2015). Earlier studies have shown that climate variables, particularly temperature, tend to be more important in explaining abundances of bird species in Northern Europe compared with land use variables (Howard et al., 2015), which supports our goal to study especially the impacts of changing climate. 2.3 | Species distribution models In the analysis, we used the Hierarchical Modelling of Species Communities approach (HMSC; Ovaskainen et al., 2017), which belongs to the class of hierarchical Bayesian joint species distribution models (JSDMs; Warton et al., 2015). We considered each visit to each survey site (line transect or point count survey) as one sampling unit. We included in the HMSC analyses data on species abundances (matrix Y of sampling units times species) and environmental covariates (matrix X of sampling units times the included covariates). As the response variable, that is components of the matrix Y, we considered the abundance of species, that is the vector of species counts for each sampling unit. To account for the zeroinflated nature of the data, we applied a hurdle approach, in which we first modelled species presence/absence data with a probit regression and then modelled abundance conditional on presence with a lognormal regression. We included all species that had the prevalence (fraction of surveys in which the species was observed) of at least 0.05, which limited the data to 120 species. As fixed effects, as components of the matrix X, we included the mean (1) winter and (2) spring temperatures of the present year, (3) the mean summer temperature of the previous year, (4) an indicator variable describing the survey method and (5) survey effort (logtransformed), and (6– 11) proportions of habitat classes within the 300m buffer. We fitted the HMSC model with the Rpackage Hmsc (Tikhonov et al., 2020) assuming the default prior distributions (Ovaskainen & Abrego, 2020, Chapter 8). We sampled the posterior distribution with four Markov Chain Monte Carlo (MCMC) chains, each of which was run for 1,500,000 iterations, of which the first 500,000 were removed as burnin. The chains were thinned by 1000 to yield 1000 posterior samples per chain and so 4000 posterior samples in total. We examined MCMC convergence by examining the potential scale reduction factors (Gelman & Rubin, 1992) of the model parameters that measure species responses to the covariates included in the X matrix. Since differences between modelling approaches are typically less impactful than other factors (especially modelling strategy, included covariates, and the study system) when evaluating temporal predictive performance (Rapacciuolo et al., 2012), we use a single model, instead of using consensus modelling approaches which put more emphasis on precision rather than accuracy (Dobrowski et al., 2011). We selected to use specifically HMSC because it performed generally the best in terms of predictive power in a recent comparison among singlespecies and JSDMs (Norberg et al., 2019). One of the benefits of this multispecies approach is that when the model is estimating the environmental responses of each species, through its hierarchical structure it can borrow information across all species (Ovaskainen & Abrego, 2020). This is likely to increase 14724642, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ddi.13687 by University Of Jyväskylä Library, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 6 | PIIRAINEN et al. predictive performance, especially for rare species for which the data alone is not sufficiently informative for accurate parameterization. Accordingly, Norberg et al. (2019) showed that HMSC models that included such feature had better predictive performance compared with corresponding singlespecies models. We report our data and modelling process according to the ODMAP protocol (Zurell et al., 2020) in Appendix S2. We emphasize that the details regarding model algorithms were secondary in this study as the focus was on the comparison of validation methods. 2.4 | Model validation methods To test the ability of the model to make predictions for the future, we first parameterized the model using data until the year 1999 only (3486 visits to 901 unique routes). We then used the model to make predictions yP for each survey site (see below for how the predictions were generated for the two components of the hurdle model), both for the period 1996– 1999 (denoted by the vector yP 1999 ) and for the period 2013– 2016 (denoted by yP 2016 ). Then, we computed analogous variables from the raw data to record the observed yO (denoted by vectors yO 1999 and yO 2016 ). For the static validation method, we measured the model's ability to predict the observed distribution for each focal species as the Pearson correlation between yP 2016 and yO 2016 over the survey sites. For the change validation method, we computed the predicted change (denoted by the vector cP ) for each survey site as the difference between the predictions for the two time periods, that is cP =y P 2016 −y P 1999 . To record the observed change (denoted by the vector cO ), we computed the difference between the true observations, that is cO =y O 2016 −y O 1999 . We then computed the model's predictive performance for each species as the correlation between cO and cP . To evaluate the presence– absence part of the hurdle model, we used as the measure y the probability of occurrence, that is proportion of years in which the survey site was occupied. For evaluating the abundance part of the hurdle model, we used the mean of logtransformed counts over the years. Predicted mean abundances were counted only for years corresponding to species' presence. We chose correlation as our measure of predictive performance because we were interested in identifying locations or environmental conditions which will influence the species especially favourably or unfavourably, that is on measures of discrimination. Evaluating predictive performance in terms of accuracy, calibration and precision would clearly be relevant as well (Norberg et al., 2019), but for the above reason we considered discrimination to be of the highest relevance and focussed on that. For occurrence data, we chose only survey sites that had been visited at least twice, either during only the latter time period (static method) or during both time periods (change method). This resulted in 1214 and 264 sites, respectively. For abundance, the data were further restricted to only those survey sites where the species was observed at least once. This resulted on average 422 sites per species (range 67– 1182) for the latter time period (static method) and 83 sites (range 3– 262) for both time periods (change method). Despite the smaller sample size and geographical coverage of validation data for the change method (Figure 2), the validation data still cover both the northern range edge and the core breeding areas for over 90% of study species (Keller et al., 2020). Species with observations from <10 sites were excluded from the validation data. Thus, in the abundance data for the change method, hazel grouse (Bonasa bonasia) and rustic bunting (Emberiza rustica) were excluded. To quantify the overall amount of change per species, we calculated speciesspecific mean observed and predicted changes from across all the study sites, which contributed to the change validation datasets (n = 10– 264 depending on species and type of data). To ask whether the model was able to predict which species increased and which decreased, we counted a correlation coefficient between the mean predicted and observed changes over the species. 2.5 | Effect of validation method To investigate the differences in predictive performance between the two validation methods, we fitted a linear mixedeffects model, separately for occurrence and abundance datasets. For both occurrence and abundance data, we used the speciesspecific Pearson's correlation coefficient, that is predictive performance, as the response variable. As a fixed effect, we entered the validation method (categorical variable: static or change). Species was included as a random effect and species' phylogeny was modelled in the variance structure to account for if closely related species showed similar predictive performance. Phylogeny tree was acquired from birdt ree.org (Jetz et al., 2012). The modelling was conducted using the Rpackage and function MCMCglmm (Hadfield & Nakagawa, 2010) with Gaussian error distribution using 1,003,000 iterations, where first 3000 were used for ‘burning in’ and thinning interval was 1000. 2.6 | Effect of species traits To investigate how the predictive performance for each validation method was related to species prevalence (continuous variable), migration strategy (long or short migration or resident), body mass (continuous) and habitat preference (forests, wetlands, cultural environments, mountains and mires), we fitted a set of 17 competing models. Prevalence and body mass variables were logtransformed prior to analysis to reduce nonnormality. Our full model included each species trait and its twoway interaction with validation method as fixed effects. The remaining 16 models included 14 models with reduced combinations of the twoway interactions between method and traits, one model with only the main effects of all variables and one model with only the main effect of validation method (all models are listed in Table S1.6: Appendix S1). Species was included as a random effect and species' phylogeny was modelled in the variance structure. 14724642, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ddi.13687 by University Of Jyväskylä Library, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 7 PIIRAINEN et al. The models were fitted using the function MCMCglmm (Hadfield & Nakagawa, 2010) with Gaussian error distribution (1,003,000 iterations, 3000 warmup iterations, thinning interval of 1000). For model selection, we used the multimodel inference approach (Burnham & Anderson, 2002) by calculating the corrected deviance information criterion (DIC) for each model and selecting the one with the lowest DIC (Spiegelhalter et al., 2002). 2.7 | Robustness test We assessed the sensitivity of our results to the effect of using a single species instead of JSDM. To do so, we considered otherwise identical HMSC models as in the case of joint species distribution modelling but fitted the models separately for each of the 120 study species. We then calculated the same predictive performance metrics with static and change methods for each species. To investigate the effect of modelling approach (single vs. joint SDM), we fitted a linear mixedeffects model where both validation method and modelling approach were used as fixed effects to explain predictive performance. Species was included as a random effect and species' phylogeny was modelled in the variance structure. The model was fitted with the MCMCglmm - function (Hadfield & Nakagawa, 2010) with Gaussian error distribution (1,003,000 iterations, 3000 warmup iterations, thinning interval of 1000). 3 | RESULTS Both HMSC models, for presence/absence data and for abundance data, reached good MCMC convergence and thus their posterior distributions were adequately sampled (Table S1.3: Appendix S1). Species occurrences and abundances were more often statistically connected with climatic variables rather than land use variables (Figures S1.1 and S1.2: Appendix S1). The overall mean speciesspecific observed and predicted changes across all study sites varied almost symmetrically around zero, thus some species declining and some increasing in their overall distribution and abundance. The observed mean changes varied between −0.11 and 0.15 in the occurrence probability and between −0.74 and 0.95 in the logtransformed abundance. The predicted mean changes varied much less, between −0.05 and 0.10 in the occurrence probability, and between −0.30 and 0.25 in the logtransformed abundance (Figures S1.3 and S1.5: Appendix S1). For most species, observed changes were more negative (less increase or greater decrease) than what was predicted (Figures S1.4 and S1.6: Appendix S1). The correlation coefficient between mean observed and mean predicted changes was 0.21 (p = .023) for occurrence data and 0.25 (p = .0056) for abundance data (Figures S1.3 and S1.5: Appendix S1), meaning that the model predicted better than by random which species declined and which increased in general. 3.1 | Predictive performances For occurrence data, the mean correlation coefficient for the static validation method was high, 0.54 (SD 0.18), whereas for the change method it was very low, 0.024 (SD 0.13; Figure 3, speciesspecific correlations in Table S1.2). For abundance data, the corresponding mean correlations were 0.31 (SD 0.19) for the static method and 0.003 (SD 0.18) for the change method (Figure 3, Table S1.2: Appendix S1). Thus, for both occurrence and abundance data, the static validation method resulted in much higher predictive performance compared with change method. When focussing on the static validation method, the mean predictive performances were higher for occurrence compared with abundance data, but for the change method no such difference was detected (Figure 3). FIGURE 3 Predictive performance, measured as the correlation between observed and predicted values, for each species (n = 118– 120). The two panels correspond to occurrence and abundance analyses, and for each panel, results using static and change validation methods are shown. The dot and the error bar show the overall mean for all species and its 95% confidence interval. Occurrence Abundance Static Change Static Change −0.5 0.0 0.5 Validation method Predictive performance (Pearson's r) 14724642, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ddi.13687 by University Of Jyväskylä Library, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 8 | PIIRAINEN et al. The posterior mean estimates for the effect of validation method were −0.51, 95% CI [−0.55, −0.48] for occurrence data and −0.31, 95% CI [−0.35, −0.27] for abundance data (Tables S1.4 and S1.5: Appendix S1). 3.2 | Species traits Habitat preference and prevalence were the most important traits affecting predictive performance in occurrence data (Table S1.6: Appendix S1). However, these effects were strongly moderated by the validation method used, as evidenced by the significant interactions between the terms (Table S1.7: Appendix S1). Species which prefer cultural (urban and agricultural areas) and mountain and mire habitats were associated with higher predictive performance but only with the static method, while no such connection was detected with the change method (Figure S1.7: Appendix S1). Higher prevalence was also associated with higher predictive performance with the static method, whereas with the change method the positive association was very weak (Figure S1.8: Appendix S1). In abundance data, habitat preference, prevalence and migration behaviour were the most important traits affecting predictive performance (Table S1.8: Appendix S1). More common, shortmigrating species preferring cultural habitats were associated with higher performances when conducting validation with the static method (Table S1.9, Figures S1.9– S1.11: Appendix S1). However, again these connections were much weaker or nonexistent with the change method. In summary, the effects of species traits were rather weak and less visible for change method, which suffers from low predictive performance in general. 3.3 | Single versus joint species distribution modelling When comparing the results from single SDMs to joint SDMs, for static method the predictive performances were much lower (but still higher than change method), whereas for the change method the performances remained near zero (Table S1.10, Figure S1.12: Appendix S1). The model results showed that the singlespecies modelling approach had a negative effect on predictive performance. The posterior mean estimate for the effect of modelling approach was −0.13, 95% CI [−0.15, −0.10] (Table S1.11: Appendix S1). In summary, our robustness test indicated that the use of joint instead of single SDM does not influence the relative order of estimates of model performances between static and change validation methods. 4 | DISCUSSION Our study aimed at developing the concept of SDM validation by comparing two different validation methods for assessing temporal predictive performance, one that measures and predicts a static distribution in a population, that is only 1 year or mean of several years, and one that measures a change, that is difference in distribution between two separate time periods. 4.1 | Static validation ignores change events Our results showed a major plummet in the estimate of predictive performance when change validation was applied, compared with the static method. This means that the interpretation of model performance changes drastically depending on which type of validation is applied, static or change. This highlights that, even if validation is done as recommended, with temporally independent data (Harris et al., 2018), the estimates of model's predictive performance might still be overly optimistic if evaluating future static distribution instead of change in distribution. The main reason for this difference lies behind the fact that species ranges are often rather static, expanding or retracting mainly at their edges. This stability increases the possibility that, even as time passes, a species is still present, and in somewhat same quantities, in most parts of its previous range. Spurious species– environment associations during model calibration may not be detected by a static performance metric if the temporal validation compares only mostly unchanged areas. Our model, on average, predicted only very moderate changes across most of the species' range. It therefore failed especially at those locations where more pronounced changes had happened, but at the same time succeeded in correctly predicting little change in most locations where distributions remained rather stable. Several studies have found similar results where static temporal validation has shown good estimates but closer inspection of occurrences at sites where changes have happened, usually at the species' range margins, has revealed that predictions of change events have failed (Eskildsen et al., 2013; Rapacciuolo et al., 2012; Sofaer et al., 2018). 4.2 | Reasons for failing to predict change events Why we failed to predict change could be that there were only minor changes happening in either bird distributions or in land cover or climate variables in the study area within the forecasted time period (approximately 16 years). The smaller the variation in the calibration data, the harder it will be for the model to estimate the underlying species– environment associations and thus to predict the spatial pattern of change events correctly. However, our bird data and climate records from Fennoscandia show that poleward shifts in both species' abundance (1.5 km/year) and climate (temperature, 7.4 km/ year) have occurred in our study area (Lehikoinen & Virkkala, 2016). In addition, logging volumes have for instance increased c. 20% in the 21st century within the study area (Virkkala et al., 2020). Therefore, we cannot conclude that our SDMs would have failed to predict changes simply because nothing had changed. 14724642, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ddi.13687 by University Of Jyväskylä Library, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License