Distribution modeling is a process which aims to understand and visualize, in a spatial context, the past, present or future distribution of a study object by relating it to predictor (explanatory) variables (Guisan and Zimmermann, 2000). The predictors can be stand-alone variables, interactions between variables, or even abstract complex-gradients (Simensen et al., 2020). This places distribution modeling among gradient analysis techniques (Halvorsen, 2012) as a correlative procedure for modeling the realized distribution of the study object (Wiens et al., 2009). Distribution modeling methodology has proliferated strongly over the last two decades (Franklin, 2010; Yates et al., 2018), and is now used to predict the distribution of single species, species assemblages, species richness, communities, and habitats (e.g., Drew et al., 2011; Moriondo et al., 2013; Horvath et al., 2019).

The diversity of research disciplines in ecology, questions asked, data collected and methods used for distribution modeling, has led to the development of several parallel distribution modeling practices (Araújo et al., 2019). These practices compete with and/or encompass a variety of terms, such as ecological niche modeling, habitat suitability modeling, species distribution modeling and occupancy modeling (see e.g., Gogol-Prokurat, 2011; Sales et al., 2021, and references therein). The practices differ in many respects, but they build on a common theoretical basis that focus on correlative predictions of distributions in general (e.g., Guisan and Zimmermann, 2000; Austin, 2007; Halvorsen, 2012). Although there are many ways of presenting distribution modeling (see Brown and Yoder, 2015; Mazzoni et al., 2015; Rodríguez-Rey et al., 2019; Zurell et al., 2020), for the purpose of this synthesis we have generalized the process into a five-step conceptual framework illustrated in Figure 1: (step 1) ecological understanding, assumptions and problem formulation; (step 2) data collection and preparation; (step 3) modeling method: choice, tuning (determining options and settings) and parameterization; (step 4) model evaluation; and (step 5) implementation and use, e.g., presentation of modeling results in the form of maps of an object’s predicted distribution, in a given area (space) and/or at a given time.

A series of important decisions are made throughout the distribution modeling process. The reliability, or uncertainty, of the output will be determined by the errors, biases and inaccuracies that have accumulated throughout the process, from problem formulation to model evaluation and use (Lindner et al., 2014). Identification of potential uncertainties at each step (Figure 1) challenges the modeler’s knowledge of the study object, the dataset and of modeling methodology in general. While the modeler is responsible for providing end users with information about the reliability of the model, the responsibility for appropriate use of the modeling outcome (which typically is a map) lies with the end user. Although the attention given to reliability in distribution modeling is increasing (e.g., Beale and Lennon, 2012; Mouquet et al., 2015; Yates et al., 2018), most studies addressing this have focused on specific parts of the distribution modeling process, such as the implication of missing environmental variables or low spatial resolution (Suárez-Seoane et al., 2014), lack of map layers that represent biotic interactions (Svenning et al., 2014; Kass et al., 2021), quality of input data (Aubry et al., 2017), the relative performance of statistical methods (Elith et al., 2006; Lobo et al., 2008), sampling bias (Varela et al., 2014; Støa et al., 2018), transferability (Yates et al., 2018) or model evaluation (Halvorsen et al., 2016; Fernandes et al., 2018). The many sources of uncertainties in distribution models strongly call for a systematic approach to avoid, reduce, or at least understand, the major sources of uncertainty at each step in the modeling process. In order to improve the quality of distribution models and facilitate model interpretation, such a systematic approach should inspect the effects of uncertainties accumulated throughout the process. The purpose of this synthesis is therefore to identify the main sources and main implications of accumulated errors, biases, and inaccuracies for the reliability of distribution models, to describe how uncertainties may be reduced, and to come up with recommendations, in the form of step-by-step guidelines, for dealing with these uncertainties. We accomplish this aim by presenting and discussing each of the five steps in the distribution modeling process (Figure 1). Our ambition is that this synthesis will be useful for both modelers and end users of distribution modeling results, such as planning authorities, nature managers and industry stakeholders. The focus of this synthesis is limited to correlative distribution models and will neither address dynamic (e.g., Horvath et al., 2021) nor mechanistic (e.g., Dormann et al., 2012), i.e., process-based, models.

A major challenge for distribution modeling has been, and to some extent still is, the lack of a unified theoretical basis and consensus on best-practice standards (but see Halvorsen, 2012; Araújo et al., 2019; Zurell et al., 2020). Incomplete knowledge of the natural variation of the study object and its environment, and lack of understanding of the ecological processes responsible for the study object’s relationship with its surroundings, may result in suboptimal decisions at some point in the distribution modeling process. Species’ distributions are controlled by four main structuring processes: (1) the inherent biological adaptations of the species, which set limits for their ecophysiological tolerance; (2) biotic interactions with other organisms, such as competition, amensalism, mutualism, and predation (e.g., Svenning et al., 2014); (3) demographic processes such as dispersal, establishment, survival and space limitation (e.g., van Groenendael et al., 2000); and (4) feedback mechanisms with the abiotic surroundings (e.g., Bonan, 2008) or engineering of the ecosystem by the species themselves (Linder et al., 2012). Lack of knowledge about causality and the object’s abiotic or biotic constraints may lead to misleading or even meaningless models.

A fundamental assumption of distribution models is the existence of an equilibrium between the study object and the surrounding environment, i.e., that the object occurs in all or a predictable fraction of suitable sites and is absent from all unsuitable ones (Guisan and Zimmermann, 2000). This is, however, not necessarily true, even for species typical of relatively stable habitats. Species may, for instance, maintain sink populations in some areas for some time or moderate the ecosystem to their specific needs (Linder et al., 2012). Moreover, environmental change causes lack of equilibrium in species abundances and time-lagged responses, in the forms of extinction debts or immigration credits (Jackson and Sax, 2010). The time-lags involved in community dynamics may vary from species to species and from area to area, depending on properties such as generation time, the mode and capacity of dispersal and the rates of environmental change (Svenning and Skov, 2005). The more extensive the environmental change, the larger extinction debt and immigration credit can be expected. This is particularly relevant for distribution models developed for evaluation of range shifts, such as altitudinal or latitudinal responses to climate change (e.g., tree lines, Bryn et al., 2013). Invasive species that have not yet reached all suitable habitats also violate the assumption of equilibrium with the environment (Barbet-Massin et al., 2018; Dimson et al., 2019).

Distribution modeling methods can be applied for a variety of reasons. The most basic purpose of distribution modeling is to produce a model that fits a study object’s distribution within an area at a specific time-point (or time interval) in the best possible way (Franklin, 2010). Distribution modeling methodology can also be used to transfer model predictions to another spatial and/or temporal setting than the one represented by the available data (Yates et al., 2018). Examples are efforts to predict effects of future climate change, for cost-efficient search for rare species, to identify areas suitable for habitat restoration, or in search for optimal locations for an activity, such as aquaculture.

Most purposes of distribution modeling can be placed along one single gradient of distribution modeling purposes. In accordance with Halvorsen (2012), we use the terms spatial prediction modeling (SPM) and ecological response modeling (ERM) for the end-points of this gradient. In SPM, the fit between model predictions and the true distribution of the study object in a specific geographical space is optimized, while in ERM the object’s response to major environmental variables, i.e., a relationship in a general ecological space is modeled. While SPM represents a purely descriptive purpose, ERM comes with an intention of understanding the study object’s distribution. The extremes along the gradient of distribution modeling purposes call for different sampling designs and methods for data collection (step 2), as well as different choices at the subsequent steps in the modeling process (steps 3–5). Because lack of a clearly formulated modeling purpose at step 1 may severely impair the reliability and applicability of distribution modeling end products, the specific problem formulation and the position of the modeling purpose along the SPM-ERM gradient could guide the choices made during all steps in the modeling process.

The lack of a unified theoretical basis and consensus on best-practice standards poses both advantages and challenges. Advantages have been creative proliferation, leading to rapid mobilization of data, development of new methods and operationalization of these methods in user-friendly software. Challenges are a large diversity of distribution modeling concepts built on different assumptions, with different demands on data quality and methods. The reliability of the resulting plethora of methods and models obviously varies, but lack of a unified theoretical basis makes them hard to compare. A statement of model assumptions and a clear formulation of the problem are necessary for transparency in the communication of modeling results.

To be able to formulate a modeling purpose precisely, a certain degree of ecological understanding of the relationship between the study object, environmental predictors and other important structuring factors (such as biotic interactions) is needed. Furthermore, it should be noted that the purpose of projecting modeling results to a different area or future climate scenarios, requires models that addresses the general relationship of the study object with the environment, i.e., an ERM model. Furthermore, knowledge of relationships in the new area or assumptions about processes that drive future development is crucial. A complete understanding of causality is, however, difficult, even often impossible, to achieve; distribution modeling basically addresses correlative relationships (see Shipley, 2016). Biotic interactions, time-lags (e.g., dispersal processes) and historical constraints (e.g., isolation due to ice age barriers) are particularly difficult to take into account in distribution modeling, as these mechanisms are difficult to represent by predictor variables. Acknowledging the lack of information on, for example, biotic interactions, as a source of reduced model performance are therefore important (Dormann et al., 2018a). We therefore recommend a conscious selection of a scale (both spatial and temporal) for the study that accord with its purpose. Of relevance in this context is that although biotic interactions are basically neighbor (local) phenomena that affect coexistence of species at finer scales, biotic interactions may occur so regularly or over so large areas that they add up to patterns that are recognizable on regional scales (Giannini et al., 2013). Moreover, when genetically and ecologically different sub-populations are to be modeled, knowledge on the ecological response at sub-population level is needed.

Violation of the assumption of equilibrium between the distribution of the modeled object and its environment reduces the reliability of the model (Chen and Leites, 2020). Violating this assumption will likely lead to under- or overpredictions by the model, depending on whether the object is a “leading edge” or a “rear end” species (e.g., Rumpf et al., 2018) and on the environment’s direction of change. Accordingly, a complex pattern of environmentally conditioned under- and overpredicted sites may result for one and the same study area. In such cases, the distribution model provides a map of the potential distribution of the study object given a specific set of assumptions (see Wiens et al., 2009). Even then, the map may be useful for specific purposes, but it should be specified that the model is predicting suitable habitats rather than a map of the realized distribution (see e.g., Gogol-Prokurat, 2011).

The properties of the modeled object and the size of the study area relative to the extent of the study object’s occurrence influence the expected level of uncertainty in distribution models (Stokland et al., 2011): weak relationships with the environment result in distribution models with lower predictability, as expected for instance for short-lived compared to longer-lived species (Hanspach et al., 2010) and for fast growing and early successional tree species compared to slow-growing species with high competitive abilities (Guisan et al., 2007). Furthermore, dispersal constraints are expected to be rare when the size of the study area is small, or if the study object is widely distributed (Soberón and Peterson, 2005).

In theory, although not yet tested in practice, compensating for violation of the assumption of equilibrium between the object’s distribution and its environment should be possible, given that the following conditions are satisfied: First, the main reason for lack of equilibrium must be recognized, for example that climate warming causes upward treeline shifts. Second, an idea about the response time or distributional time-lag of the study object is needed. For example, the time needed for treeline trees to disperse and grow tall at higher elevations. Third, historic and up-to-date climate data as well as the response of the study object must be available at the resolution relevant for distribution modeling, for example high-resolution data layers on temperature and monitoring data of treeline dynamics. If these assumptions hold true, time-lag corrected predictions can be obtained by parameterizing a model using response data that represent the present distribution and appropriate historical climate data, and project to the present climate. An example of how the response- or predictor data can be used to compensate for non-equilibrium situations is given in Figure 2. For species with rapidly advancing ranges, such as invasive species, the equilibrium assumption may be tested by building retrospective models (Barbet-Massin et al., 2018; Dimson et al., 2019).

According to Osborne and Leitão (2009), sampling bias is the biggest potential source of error in response variable data for a study object. Other errors are related to imprecise or wrong georeferencing and direct errors such as misidentification of species and typing (punching) errors. All of these shortcomings introduce uncertainties into distribution models. While the accumulated impact of direct errors on the modeling result is simply proportional to the number of such errors, georeferencing errors/inaccuracies and sampling bias have more complex effects.

Georeferencing inaccuracies. We often assume that the study object is recorded without georeferencing errors (Osborne and Leitão, 2009). However, positioning inaccuracies are, to some degree, present in all georeferenced data, e.g., due to imprecision in the equipment used for recording positions. This is the case for museum records (Bloom et al., 2017), but in particular for data from citizen science projects, where considerable variation in the quality of the equipment (GPS) used for georeferencing must be expected. When mapping objects at sea, weather conditions, such as strong winds, waves and currents, will impact the precision of the recordings (Figure 3). Overall, positioning accuracy and reliability of georeferenced data has improved considerably over the last decades, after GPS equipment came into common usage. This will introduce differences in georeferencing accuracy among old and new observations. Also, the effect of positioning inaccuracies may vary in space, e.g., be higher in areas with large variation in environmental conditions over short distances, such as in steep terrain, or within the narrow tidal zone along the coast.

Spatial autocorrelation is the systematic spatial variation in a variable, i.e., the tendency of observations from closely situated sites to be more similar (positive spatial autocorrelation) or less similar (negative spatial autocorrelation) than observations made further apart. Positive spatial autocorrelation is most common in data representing ecological phenomena. Spatial autocorrelation in data used for training as well as evaluation of distribution models poses challenges, and has been shown to impact both coefficients and inference in statistical analyses (e.g., Dormann et al., 2007 and references therein). This is exemplified by Tognelli and Kelt (2004) who found that the relative importance of explanatory variables may depend on the presence or absence of spatial autocorrelation in the data. Halvorsen et al. (2016) argue that challenges with spatial autocorrelation in the context of distribution modeling arise from autocorrelated residuals in the response variable and point to the importance of accounting for the spatial autocorrelation by predictor variables in the model.

Sampling bias is a broad term that encompasses many types of biases. We find it useful to distinguish between spatial, temporal and study-object bias; see discussions in Robertson et al. (2010) and Støa et al. (2018). Spatial bias is very common and occurs when some parts of the geographical and environmental space are sampled more intensively than others or when object detections are imperfect (see e.g., Lahoz-Monfort et al., 2013; Tobler et al., 2019 and references therein). That is, where the sampled frequency distribution deviates from the objects’ true distribution in the geographical and environmental space (Støa et al., 2018). This typically applies to deep oceans, high mountains, extremely strongly wave-exposed sites, war zones or other sites that are difficult to access. Spatial bias in the form of oversampling is likely to occur in particularly easily accessible locations (e.g., Wollan et al., 2008), for example as a result of sampling by citizen science (Hughes et al., 2021). Temporal bias occurs if an object is more intensively sampled during some periods (e.g., summer or daytime) than others or if, e.g., new museum material is digitized and made available before older material. Study-object bias results from specific properties of the study object (e.g., species or habitat). Thus, rare and/or inconspicuous species may be represented by too few observations to give an appropriate picture of their distribution (e.g., Varela et al., 2011) while more charismatic, “iconic” or commercially interesting ones are more prone to oversampling. Study-object bias influences the reliability of model comparison studies.

Contrary to the response data, which are sampled as point data, wall-to-wall data layers of the whole study area is needed for wall-to-wall predictive maps (Franklin, 2010). This is an important constraint which typically precludes many relevant predictors from being available for distribution modeling. Also, a predictor layer may itself be the output of a modeling process, e.g., climate data interpolated from observations at weather stations or composite layers (Simensen et al., 2020; Tessarolo et al., 2021). Depending on their origin, explanatory (predictor) data used for distribution modeling are subject to many sources of inaccuracies and errors, such as imprecise or erroneous variable values and inappropriate resolution (in space or time). Too often, the predictor variables used in distribution models come with unknown degrees of uncertainty. Typical examples are WorldClim data for terrestrial modeling (Fick and Hijmans, 2017), Bio-ORACLE for marine modeling (Assis et al., 2018) and variables derived from digital elevation models (DEM), as well as land-cover and land-use maps. When many such variables are used to build distribution models, uncertainties accumulate to an unknown magnitude.

Lack of explanatory data needed to predict an object’s distribution accurately may give rise to negative results (models with no or low predictability). This is a type of “omission error” that is likely to be underreported because of the tendency for a positive publication bias, i.e., that negative results fail to get published (Kicinski, 2013). Factors known to be important for the distribution of a species cannot be included as a predictor in the analysis unless available as wall-to-wall maps.

Lack of ecological relevance of the explanatory variables included in a distribution model is common, attributed to the fact that when variables that adequately represent the factors likely to govern the distribution of the study object are not available, we use proxies or substitutes instead. Examples of predictors that are hard or impossible to obtain wall-to-wall coverage for, are nutrient availability, oxygen level or substrate type in marine systems (Bekkby et al., 2008; Rinde et al., 2014) and water deficiency in terrestrial systems (Slette et al., 2019). Lack of relevant predictors will inevitably reduce the explanatory power of a model. Although use of proxies or substitutes that are correlated with the causal factors may improve the predictive ability of spatial prediction models, they are inappropriate for ERMs unless their relationship to the putatively causal factors is well known. Moreover, inclusion of redundant predictors may lead to errors in parameter estimates (because of collinearity), reduce precision in parameter estimation and lead to overfitting of SPM models (see e.g., Beale and Lennon, 2012; Dormann et al., 2013; Cade, 2015) and make ERMs inapplicable for their purpose.

Imprecise or erroneous predictor variable values may, for instance, be introduced by spatial interpolation of environmental models (Peters et al., 2009) based on a limited number of observations, of which each comes with its own uncertainties. The quality and relevance of modeled predictor variables used in distribution modeling studies are rarely discussed, and deserve more attention.

Inappropriate spatial resolution of important explanatory variables may impact the final model. Contrary to the response data, which are typically registered at a fine resolution in a restricted study area, predictor layers are often generated from DEM, remote sensing products or interpolated from climate data, and therefore usually available as broad-scale, low-resolution maps. Mismatch between the spatial resolutions of predictor and the response variable(s) often results in a modeling of fine-scale distributions by use of coarse-scale explanatory variables.

Inappropriate temporary resolution is exemplified by predictor variable data that do not represent the environmental conditions at the time when the distribution of the modeled object was shaped. Averaged or composite climate variables have low explanatory power when the object’s distribution is restricted by past weather extremes (e.g., Wernberg et al., 2012). The same may be the case when historic distributional data are modeled as a response to predictors that describe present environmental conditions.

Georeferencing inaccuracies should always be reduced to a minimum. This can be accomplished by using high-quality, state-of-the-art GPS equipment. However, an off-the-shelf GPS with an accuracy of 2–10 m will be sufficient for most distribution modeling purposes. If historic data (with low georeferencing quality) or data sampled under suboptimal conditions (e.g., bad weather or turbulent sea) have to be used, reliable estimates of the uncertainty should be taken into account when the distribution map is interpreted. This applies both to the model builder and the model user. The impact of georeferencing uncertainty on distribution modeling results can be estimated by simulating random location errors (see Graham et al., 2008). A general recommendation is to choose a pixel resolution that is large enough to reduce the impact of georeferencing uncertainty while at the same time accord with the modeling purpose (see e.g., Gábor et al., 2020). Finally, standardized protocols for correction of inaccurate geographic coordinates may be used, for example the SAGA protocol proposed by Bloom et al. (2017) for museum data. Such protocols are useful for correction of already sampled data, but a more provident solution is to develop rigorous protocols for sampling of new data so that the problem itself diminishes with time.

Although recent studies have indicated that spatially autocorrelated response data for the study object influences modeling results less than previously assumed (Halvorsen et al., 2016), spatial autocorrelation may still represent a challenge. Several measures can be taken to reduce the effects of this phenomenon. For instance, a number of studies enforce minimum distances between presence observations (e.g., Bryn et al., 2013) or buffering methods to avoid spatial autocorrelation (e.g., Simensen et al., 2020), whereas other studies apply systematically sampled (area-representative) training data (e.g., Horvath et al., 2019). An alternative approach, which is based on the assumption that spatial autocorrelation is less problematic than previously assumed, is “background thickening” (Vollering et al., 2019a). By this method, the density of randomized uninformed background observations to be used in presence-only distribution modeling is increased or decreased locally to match the density of presence points.

In order to reduce sampling bias, knowledge about the study object’s autecology is crucial. A basic requirement for good distribution models is that the sampling design provides accurate information about the frequency of presence (FoP) of the modeled object along all important environmental gradients (Merow et al., 2013; Halvorsen et al., 2015; Støa et al., 2018) and that the sampling matches the detectability of the study object. Irregularities in the sampling, i.e., parts of important environmental gradients that are under- or overrepresented in the sample, may be detected by comparison between observed and theoretical FoP curves (Støa et al., 2018). Identifying areas or contexts where the model is not valid, i.e., parts of gradients not adequately represented by data, will help to identify the problem. When already sampled data, for instance from museum collections or public databases (e.g., GBIF), are used for distribution modeling, all available information about sampling design and other methodological aspects related to the sampling should be used to assess potential sampling bias inherent in the dataset. Ideally, only external data that are described properly should be used. Fortunately, the number of peer reviewed data papers are increasing rapidly, e.g., in GBIF from 1 in 2011 via 13 in 2014 to 27 in 2018 (GBIF, 2019).

Small datasets come with more stochasticity, which normally entails higher difficulties in obtaining significant and interpretable patterns by modeling. Ideally, data should be collected specifically for each modeling study, based on a tailor-made sampling design that matches the specific modeling purpose. The spatial pattern of objects with “clumped distributions” and specific environmental demands can generally be modeled with greater precision than patterns of widely distributed objects (Carrascal et al., 2006; Hernandez et al., 2006; Stokland et al., 2011). In order to model “ubiquitous” objects, the spatial resolution of the predictor variables should be high enough to capture variation along relevant local environmental variables (Halvorsen, 2012). MacKenzie et al. (2002) and Tyre et al. (2003) provide methods for taking imperfect detection of small, rare and inconspicuous organisms into account when building models.

Lack of access to explanatory variables that may directly represent factors known to be important for an object’s distribution can be dealt with by using available proxies. This could be ocean depth as a proxy for solar irradiance at the sea floor, or altitude as a proxy for terrestrial temperature. However, the use of proxies reduces the interpretability of the studied object’s relationship with the environment, as many proxies are strongly correlated with several predictor variables. Global environmental models, such as those provided by WorldClim (Fick and Hijmans, 2017) and Bio-ORACLE (Assis et al., 2018) for terrestrial and marine environments, respectively, have rapidly become very popular, and are used in many bioclimatic studies. Still, they should be used with caution, for several reasons. First, such models are often developed from data with unknown levels of inaccuracy. Second, generic data may not match the appropriate scale for a specific study and hence be more or less irrelevant, or even direct the study toward a suboptimal choice of spatial resolution. Advantages of using easily accessible, modeled predictor data are that they have typically been described, tested and analyzed with respect to uncertainty in many studies, and that much experience has been gained by their use for practical distribution modeling (e.g., Bedia et al., 2013). Furthermore, in general we think that the use of remote sensing based explanatory variables could be exploited more in ecological distribution models (see e.g., Tang et al., 2020).

One way to deal with imprecise predictor variables, or biased response data, is to create a confidence map, which provides the user with information about where a distribution model is judged or estimated to be reliable and where it is less so (e.g., Hemsing and Bryn, 2012). A confidence map may, for instance, contain information about the number of data points that back up each pixel value, areas with poor coverage of satellite data, areas where comparable datasets produced by others deviate, or areas where predictions from models obtained by different methods deviate.

The use of correlative (i.e., not causative) predictors in space or time can be justified if the aim is to obtain the model that best predicts patterns of presence in a certain area in a specific time interval, i.e., SPM. In these cases, the correlative predictors may serve as useful proxies for missing relevant variables. However, if the aim is ERM, in order to understand ecological relationships or to transfer (extrapolate) predictions to another spatial and/or temporal setting than those represented by the data, great caution should be taken when using non-causative variables in the modeling.

To avoid biologically unrealistic or otherwise unreliable modeling results and interpretations, a deep ecological understanding of the study object is needed (Figure 5). Thus, the choice of modeling method and the tuning of different parameters to optimize predictive accuracy should be guided by the modeling purpose and the properties of the study object and the types of response and predictor data used. The type of response data, such as PO, PA, ordered factor, or continuous, call for different choices of method. Methods for PO data call for particular attention to the sampling design for uninformed background or pseudo-absence observation, which may heavily impact modeling results (Stokland et al., 2011; Støa et al., 2018). Several alternative sampling strategies have been proposed, for example “target-group background,” “presence thinning,” and “background thickening” (Phillips and Dudík, 2008; Støa et al., 2018; Vollering et al., 2019a). In general, model performance is more affected by false positives than false negatives (Fernandes et al., 2018).

Different distribution modeling approaches rely on different assumptions and have different advantages and disadvantages (see discussion in Shabani et al., 2016). For instance, some methods lean on a priori assumptions about the residual distribution (e.g., normal, binomial or poisson) or response curve shapes (e.g., linear or non-linear); some methods are suitable for capturing non-additive behaviors and complex interactions; some produce complex and often overfitted models that open for spurious interpretations while others may tend to produce too simple models. Some methods are robust to noise, some require large datasets, whereas others may perform well also with small samples (see Heikkinen et al., 2006). Furthermore, although many modeling methods are flexible with respect to data properties, available data and data formats will reduce the number of appropriate modeling methods. For instance, some methods handle presence-only data well (e.g., MaxEnt), others are developed for analyzing presence-absence data (e.g., GAM), whereas others again handle ordinal response variables (e.g., CLM). Explorative statistical analyses of data, such as calculation of descriptive statistics and analysis of frequency-of-presence (FoP) curves, should be standard elements in all distribution modeling studies, conducted before any definite decisions about modeling setup is made.

Many distribution modeling methods come with considerable flexibility with respect to tuning for complexity. In general, the choice between simpler or more complex models should be guided by modeling purpose; a complex SPM model is not overfit if its predictive ability is better than any simpler model, while ERM models are overfit if they include predictors that are not generally important for the distribution of the study object (Halvorsen, 2012). The tuning options of, for example, the regularization multiplier (Warren and Seifert, 2011) or the significance level of the F-ratio test in MaxEnt (Halvorsen, 2013) or the smoothing function in GAM (Guisan et al., 2002) are examples that control the fit of the model to data and modeling purpose. Empirically based advices on how to make these choices are, among others, given by Merow et al. (2013) and Halvorsen et al. (2015). Building models with an appropriate degree of complexity (fit) is critical for robust inference (Merow et al., 2014). Complexity indicators such as number of parameters and shapes of predicted relative FoP curves (Støa et al., 2018) are important tools for assessment of model reliability.

Model averaging or ensemble modeling, with the purpose of generating a unified (and agreed-on) model are increasingly used in distribution modeling studies (Buisson et al., 2010; Dormann et al., 2018b). Ensemble modeling may reduce the effect of biases, inaccuracies and errors in each component model. However, the merging of results from different approaches does not resolve fundamental challenges such as modeling with biased training data, missing explanatory variables, model misspecification or under- and overfitted models. Merging the outputs from many models may make detection of such pitfalls more difficult, among others by signaling degrees of precision and quality that are not supported by the models. A general advice is therefore to run distribution models with different methods, not with the purpose of obtaining and presenting one unified model, but as a means for comparing results and testing model robustness. For averaging and ensemble modeling, it would be beneficial to establish protocols for how to present the material, methods and outputs. An example of this is the World Climate Research Programme (WCRP), where the Coupled Model Intercomparison Project (CMIP) was established to standardize multi-model output formats of Earth System Models, as well as other goals (Eyring et al., 2016).

Evaluation by use of fully independent data is the only way to completely circumvent pitfalls related to sampling bias in a dataset (Austin, 2007; Veloz, 2009). Most likely, an independently collected evaluation dataset will not include the same errors, biases and inaccuracies as the training dataset, but it may come with its own shortcomings. The ideal dataset for evaluation of distribution models is collected by a sampling strategy that avoids sampling bias of any kind. Many uncertainties of models that are based upon PO data can only be circumvented by access to independent PA data for model calibration (Merow et al., 2013; Halvorsen et al., 2015); i.e., assessment of the numerical accuracy of model predictions. According to Harrell et al. (1996), independent presence-absence data are required for any kind of model calibration, i.e., the conversion of model predictions from a relative (RPPP) to a predicted probability scale (PPP).

Independent P-A evaluation data can be obtained by a wide spectrum of different approaches, e.g.: (1) Stratified random sampling, using modeled RPPP values as a criterion for stratification (Halvorsen et al., 2015). (2) Environmentally stratified sampling, by which all levels along all environmental gradients important for the study object is deliberately captured. (3) Spatially systematic sampling involves the selection of observations according to a pre-designed spatial sampling. This method avoids over- or underrepresentation (spatial clustering or the converse) of data points (Horvath et al., 2019). These three sampling procedures have different advantages and disadvantages. While approach (1) may be optimal for testing one particular model, it is less useful for testing other models. Also, it implies a two-phase sampling procedure. Environmentally stratified sampling (2) may be theoretically optimal because it provides balanced data for calculation of frequency-of-presence (FoP) curves (e.g., Halvorsen et al., 2015), but difficult to carry out in practice when models are complex with many predictor variables. Both of approaches (2) and (3) may fail to provide enough data for valid evaluation of models for study objects with low prevalence due to spatially aggregated and/or restricted occurrences. Furthermore, these approaches will be practically unmanageable for study areas with large extent and/or contains a large fraction of remote or otherwise inaccessible areas.

When the modeling purpose is predictions in time, and evaluation by independent data is not possible, comparing models trained with different data and parallel use of several modeling methods is pivotal (e.g., Elith et al., 2006).

To make a distribution model and model-derived results useful for other purposes and other users than originally intended, these auxiliary users should be provided the basic information needed to make well-informed judgments of the appropriateness of the product for their intended use. Models should therefore be accompanied by a fact sheet that includes understandable information on the original modeling purpose, the data and the methods used, and, if relevant, specific restrictions to the generality of the model. Confidence maps and other documentation or judgments of biases, errors and uncertainties should be standard elements in the fact sheet. The ODMAP protocol suggested by Zurell et al. (2020), or the INSPIRE data infrastructure for mapping and modeling in EU, which also includes standards for metadata¹, may be relevant.

The users of distribution models need to be reminded that a distribution model is not a true map of the object’s distribution or abundance, but rather mapped predictions of the (relative) potential for a given object to be present in each grid cell based on the available explanatory variables and response data. The study object may well be absent from a grid cell with high predicted probability of presence, or present in a cell with low predicted value. This uncertainty calls for special care when models are used for prioritization for conservation purposes (Schuetz et al., 2015). Information about the predictive performance of the model (e.g., Mouton et al., 2010), preferably obtained by evaluation on independent data, is therefore essential.

Maps resulting from distribution modeling can be tailored to different purposes by, for instance, selecting different cut-off values, i.e., threshold values used to transform predictions from a continuous scale to binary predictions (presences and absences, Figure 6; Scherrer et al., 2020). It is often assumed that binary predictions makes prediction maps easier for managers and planners to understand and use than continuous maps. Binary maps based upon a high threshold value may be useful for users who want to identify areas of high probability of finding the study object (Figure 6). However, for users who want to identify areas that maximize the chance of covering the range of a study object, a low cut-off value is more appropriate (see Freeman and Moisen, 2008). It should be noted, however, that conversion from continuous to binary predictions inevitably implies loss of information. Another possibility is to refine or post-process the distributions models afterward (e.g., Kass et al., 2021), so that the models are improved according to the specific needs.