Balancing structural complexity with ecological insight in Spatio‐temporal species distribution models

The potential for statistical complexity in species distribution models (SDMs) has greatly increased with advances in computational power. Structurally complex models provide the flexibility to analyse intricate ecological systems and realistically messy data, but can be difficult to interpret, reducing their practical impact. Founding model complexity in ecological theory can improve insight gained from SDMs. Here, we evaluate a marked point process approach, which uses multiple Gaussian random fields to represent population dynamics of the Eurasian crane Grus grus in a spatio‐temporal species distribution model. We discuss the role of model components and their impacts on predictions, in comparison with a simpler binomial presence/absence approach. Inference is carried out using Integrated Nested Laplace Approximation (INLA) with inlabru, an accessible and computationally efficient approach for Bayesian hierarchical modelling, which is not yet widely used in SDMs. Using the marked point process approach, crane distribution was predicted to be dependent on the density of suitable habitat patches, as well as close to observations of the existing population. This demonstrates the advantage of complex model components in accounting for spatio‐temporal population dynamics (such as habitat preferences and dispersal limitations) that are not explained by environmental variables. However, including an AR1 temporal correlation structure in the models resulted in unrealistic predictions of species distribution; highlighting the need for careful consideration when determining the level of model complexity. Increasing model complexity, with careful evaluation of the effects of additional model components, can provide a more realistic representation of a system, which is of particular importance for a practical and impact‐focused discipline such as ecology (though these methods extend to applications for a wide range of systems). Founding complexity in contextual theory is not only fundamental to maintaining model interpretability but can be a useful approach to improving insight gained from model outputs.


| INTRODUC TI ON
The continuing increase and the improvement both of the availability and detail of ecological information, and of computational resources allows realistically complex and flexible statistical models to be fitted to ecological data. However, increasing structural complexity through the inclusion of additional model components without ecological justification will not improve our ability to understand fundamental ecological responses and predictions. Hence increased model complexity is only relevant if the more complex models provide us with additional insight into ecological processes than simpler and often computationally cheaper models (Bolker, 2009).
For instance, more complex models may improve the accuracy of predictions or abundance estimates: they may have less restrictive assumptions, thus avoiding incorrect inference and model interpretation, and they may have the ability to explicitly reflect complex ecological processes.
In the specific context of species distribution models, which seek to provide insight into the distribution of species in geographic and environmental space, complex spatio-temporal models are often utilised. Increased model complexity is supported by the fact that those complex models may reflect the complex ecological processes driving the species distribution, such as demography (Pagel & Schurr, 2012), dispersal (Elith & Leathwick, 2009;Iverson et al., 2004), and physiology (Buckley et al., 2011) and hence provide additional ecological insight. Increasing model complexity through including model structures that reflect spatio-temporal dynamics not only accounts for autocorrelation not explained by covariates, but may also reflect ecological processes such as dispersal limitation or site fidelity. However, caution is needed to avoid such random effects being overfitted and concealing the impact of model covariates (Sørbye et al., 2019).
Model complexity is occasionally influenced by the maximisation of predictive performance, not ecological theory, as has been noted for correlative species distribution modelling (Austin, 2002(Austin, , 2007.
In many studies of species distribution, the research aims include interpretation of effects or mechanisms, or extrapolation to spatial areas or time periods outside of those studied. In these instances, model interpretability and ecological realism should be considered alongside predictive performance during model selection.
In this paper, we highlight that the level of complexity in a spatiotemporal ecological model has to be carefully chosen and tailored to the system under investigation. To this end, we carefully inspect the complexity and limitations of a marked spatio-temporal point process model related to a model that has previously been fitted (Soriano-Redondo et al., 2019) with the aim of understanding and predicting the spatial distribution of a reintroduced species. The modelling in Soriano-Redondo et al. (2019) is facilitated through a Bayesian hierarchical modelling approach where spatio-temporal dynamics in behaviour are represented by a flexible stochastic partial differential equation (SPDE) (Lindgren et al., 2011). This allows the model to be fitted in continuous space (Simpson et al., 2017) and within a realistic timeframe through the computational efficiency of integrated nested laplace approximation (INLA) (Rue et al., 2009).
We revisit this approach here, fitting a similar model using the user-friendly software package inlabru, with the specific aim of reflecting on the role of complex model components and how they represent ecological processes in species distribution models. We show that model complexity has to be carefully aligned with the complexity of the ecological processes that can be observed in the data, to provide adequate model interpretability and hence relevant ecological insight.

| Approaches to species distribution modelling
Species distribution models (SDMs) link information on the presence/absence or abundance of a species to environmental variables to predict where (and how much of) a species is likely to be present in unsampled locations or time periods (Martinez-Minaya et al., 2018).
Existing methods for species distribution modelling include: approaches developed to deal with presence-only datasets (such as maximum entropy algorithm, distance sampling, similarity, and envelope methods such as MAXENT, Gower metric, Mahalanobis distance, and ecological niche factor analysis); machine-learning algorithms that are iterative in ecological systems (such as artificial neural networks and deep learning); and classic regression models (such as generalised linear models (GLM), their semiparametric extension, generalised additive models (GAM), and a related method, multivariate adaptive regression splines (MARS)).
Several of the above methods are based on the assumption that the observations being modelled are independent. This is not always the case because ecological (and particularly species distribution) data usually present residual spatial or spatio-temporal correlation structures (Kneib et al., 2008). These structures can be the effects of processes that are unaccounted for in simple habitat models, such as social interaction with conspecifics, dispersal limitations or interspecific competition (Lichstein et al., 2002;Miller, 2012;Storch et al., 2003). Spatial data are also subject to spatial autocorrelationthe phenomenon that two points are more similar the closer they are in space, (Tobler, 1970).
Failing to account for unexplained spatial or temporal correlation structures in species distribution modelling can lead to spurious significance and an increase in the rate of Type I errors (false positives) in the interpretation of explanatory variables (Dormann, 2007;Dormann et al., 2007;Miller, 2012;Segurado et al., 2006). In some INLA, marked point process, spatio-temporal model, species distribution model cases, not including these effects can lead to inference of an incorrect direction in the relationship between environmental covariates and response (Kühn, 2006). However, the inclusion of model components to account for unexplained spatial or temporal structures can rectify the above issues, as well as allowing insight into otherwise overlooked effects (Lichstein et al., 2002) and improving model fit (Dormann, 2007). Although, careful prior choice is advised to reduce the risk of correlation structures incorrectly lowering the estimated significance of fixed effects (Hodges & Reich, 2010;Sørbye et al., 2019).
GAM and MARS can account for spatio-temporal autocorrelation by taking a semiparametric approach through including smoothing spline terms. A very powerful and flexible alternative is to use a parametric approach in which the spatial and temporal correlation structures are represented by a (Gaussian) random field and approximated by an SPDE (Lindgren et al., 2011;Miller et al., 2020). This allows us to flexibly model species distribution data in continuous space as georeferenced data (for count or presence-absence data) or as point patterns (for presence only or effort-corrected presence data; Yuan et al., 2017). This approach is taken in Soriano-Redondo et al. (2019) as well as in this paper.

| Point process methodology and INLA
As model complexity increases-for example, when spatiotemporal correlation structures are included as complex model components-an advantageous approach can be found in Bayesian hierarchical modelling (Cressie et al., 2009). Hierarchical models allow parameters to vary at more than one level via the introduction of random effects. Constructing such models in a Bayesian framework enables full inference through the quantification of parameters and uncertainty; a feature of great practical benefit in applied conservation contexts (Wade, 2000;Wintle et al., 2003).
Bayesian modelling also enables the inclusion of informative prior information, which can help to avoid overfitting of random effects (Sørbye et al., 2019). Several authors have applied a Bayesian approach to analyse different aspects of species distribution, such as comparing Gaussian processes against other popular approaches to modelling (Golding & Purse, 2016), visualising spatial patterns with Bayesian hierarchical logistic regression , analysing species distribution using presence/absence data , working with occupancy models (MacKenzie et al., 2002), or analysing distance sampling data using a Bayesian GAM (Sigourney et al., 2020).

Hierarchical Bayesian models have traditionally relied on
Markov chain Monte Carlo (MCMC) simulation techniques, which are computationally expensive and technically challenging, consequently limiting their application. However, INLA methodology and its powerful application to modelling complex datasets has opened hierarchical Bayesian modelling to a wider range of applications (Illian et al., 2013). As opposed to MCMC simulations, INLA uses an approximation for inference and hence avoids the intense computational needs, convergence, and mixing problems sometimes encountered by MCMC algorithms. Moreover, INLA can be used in combination with the SPDE approach, which approximates Gaussian random fields with a Matérn covariance structure, allowing for computationally efficient modelling in continuous space (Lindgren et al., 2011;Simpson et al., 2012). Implementation of these methods in R (R Core Team, 2021) is made accessible through the package inlabru (Bachl et al., 2019), which was designed as a user-friendly wrapper and extension to R-INLA (Rue et al., 2009), providing easier syntax for point process model fitting as well as reflecting ecological observation processes. These statistical methods and their packages for implementation constitute an efficient and reliable framework for complex spatio-temporal modelling, with a range of applications in ecology and beyond (Williamson et al., 2021;Yuan et al., 2017).
Here, we focus on the application of spatio-temporal marked point process models to analyse complex ecological systems, made a Gaussian random field. Given this intensity, the point locations are assumed to be independent, that is, to follow a Poisson process. As a result, these are a special case of latent Gaussian models, the class of models that may be fitted with INLA (Cox & Isham, 1980;Illian et al., 2012;Møller et al., 1998).
An extension of the simple point process model is the marked point process model, in which information on the location of objects or events is modelled along with the qualitative or quantitative properties (commonly referred to as 'marks') of these events or objects.
If a dependency between the values of a mark and the point distribution is assumed, these can be jointly modelled with two or more dependent likelihoods. Simultaneously modelling the distribution of points (as a point process) and some data feature collected at the point locations greatly increases the model complexity and computational cost. The computational efficiency of INLA, however, makes fitting these models feasible (Illian et al., 2012) and inlabru facilitates the fitting of both simple as well as marked point processes, which can be rather cumbersome in plain R-INLA.
When fitting a marked point process model in inlabru, the spatio-temporal structure of the marks (independent of the point distribution) and the spatio-temporal structure of the points can be represented with different Gaussian random fields in a shared representation of continuous space and discrete time. In this way, multiple data features and different sources of spatial clustering can be incorporated into a single model of species distribution.
In this paper, we reflect on a modelling structure first developed by Soriano-Redondo et al. (2019), which follows the population spread of the Eurasian crane Grus grus across England. Here, we construct several models based on this work, to assess the role of complex model components in capturing spatial and spatiotemporal clustering dynamics in the species distribution. We take advantage of the new and advanced analytical approaches developed in inlabru, to explore the spatiotemporal structures and AR1 temporal correlation component in Soriano-Redondo et al. (2019). Here, the point pattern of habitat patches also reflects the observation process, as data-collectors used prior information on species habitat preference to influence the sampling design. Accounting for spatially varying detection probability is a particular strength of inlabru, which was developed specifically for (ecological) datasets with complex observation processes.

| Data
We investigate the spatial distribution of a resident breeding population of Eurasian crane in England following the return of the species to the UK in 1979 (Stanbury, 2011), with the aim of predicting the distribution of the population in future years. Breeding pairs of cranes are only found in wetland habitats, which creates a spatial restriction on the suitable habitat available to the species, although little is known about its dispersal ability. The distribution of wetland locations in England is heterogeneous, with various underlying factors creating a clustered structure of habitat patches. Thus, the distribution of the crane population throughout the UK and hence its future spread is dependent on the availability of suitable wetland habitat. Due to this known species habitat preference, data on crane presence were collected only at wetland locations, and not in interstitial areas. We demonstrate how this type of data can be modelled as a marked point process in inlabru, in order to infer habitat preferences and predict species distribution, while accounting for known species habitat requirements and observation  (Stanbury & Sills, 2012). Wetland patches smaller than 5 ha were also removed, as estimates show that cranes required at least 8 ha to nest (Johnsgard, 1983). The fine grain of this classification allows us to treat the wetland areas as points, since extensive wetland regions are fragmented by other types of habitats, roads or other human infrastructure, and are indeed composed by multiple smaller wetland areas. Each wetland area has associated information: its coordinates, data on the presence or absence of breeding pairs of cranes, wetland extent, perimeter-to-area ratio, the proportion of urbanised areas in a 10 km surrounding terrestrial buffer.
To evaluate the structure of the model in Soriano-Redondo et al. (2019), we present two modelling approaches with the aim of estimating the spatio-temporal distribution of cranes in England.
The first approach makes use of a single Gaussian random field to represent spatio-temporal correlation in the distribution of the observed population of cranes. The second is a marked point process approach, which incorporates a second Gaussian random field, accounting for the distribution of wetlands. Within these modelling frameworks, we explore differences arising from the inclusion of a temporal correlation structure via a first-order autoregressive (AR1) process, where estimates for each year are dependent on the immediately previous year, compared to treating each year of data as independent and identically distributed (IID).

| Single-field models
In order to explore what level of model complexity is needed to answer relevant ecological questions based on the crane data, we start with a relatively simple spatio-temporal model in continuous space.
To improve our understanding of the spatio-temporal distribution of cranes across England, we initially ignore the fact that cranes will only nest in wetland habitats, and we first construct a model containing a single spatio-temporal Gaussian random field. This is essentially a model of geo-referenced (binomial) data (often referred to as a geostatistical model) with a single likelihood, which does not take the observation process into account and assumes that crane presence is equally likely across England. Here, crane presence is modelled as a binomial distribution: The occurrence of cranes O s,t at each location s, in each year t can be considered an individual Bernoulli trial with probability P(s, t) .
Here, the location s is a vector representing coordinates s 1 and s 2 in two-dimensional space. Note that for larger-scale models, location could also be represented on the sphere (Simpson et al., 2016).
The probability of crane presence P in a location s at time t is dependent on the observed distribution of the existing population in space and time, as well as the effect of environmental covariates: (1) O s,t ∼ Bernoulli(P(s, t)).
Here, 0 represents an intercept term and x i (s, t) is the value of each

| Marked point process models
The binomial presence/absence model in Equations (1) and (2) captures the spatial correlation structure in the occurrence of breeding pairs of cranes, but ignores the distribution of suitable habitat patches (wetlands). However, sampling only took place in locations where cranes would be expected to nest, this observation process can also be interpreted as preferential sampling (Diggle et al., 2010;Pennino et al., 2019). In this example, crane detectability is assumed to be known, as due to the large size of the animals, they are easy to detect. However, with rare or small organisms that are more difficult to detect, it may be of interest to use models to predict and estimate abundance across space. In

| RE SULTS
In both the binomial presence/absence model (Equation 2) and marked point process model (Equation 4), which incorporated an AR1 temporal correlation structure, the correlation parameter for the AR1 process, t , (which is bounded between −1 and 1) was estimated to be extremely high ( Table 1). This indicates a very strong temporal correlation across years. Due to this strong temporal correlation, the model estimated the same spatial structure in the mark random field M(s, t) and predictions of crane probability of presence across all 5 years considered in the study (see Appendix A, Figure A2). These estimates constituted a spatial structure reflecting an average of all years of data, as opposed to accurately representing changes in spatial structure each year. Such averaged estimates may be inaccurate when interpreting within-year spatial effects, and so graphical representations of results are provided here for the models which treat each year of data as IID only.
In the marked point process and binomial presence/absence models, the mark random field M(s, t) represents the spatial distribution of the crane population each year (Figure 1, centre plot). This constrains predictions to near where nesting pairs of cranes have been observed, and without it, the species could be predicted to occur in areas of suitable environmental conditions, but that are an unrealistic distance from the established population.
The marked point process models also contain a point random field, G(s), which models the spatial distribution of wetlands ( Figure 1, left plot; Figure A4). The strength and direction of the interaction between the wetland density and probability of crane presence is estimated through the scaling parameter . In the models containing the point random field, is estimated as a significant positive effect (Table 1). This means that probability of crane presence is estimated to be higher in areas of high wetland density, as opposed to low density areas. Without the inclusion of the point random field (as in the binomial presence/absence model), predictions of species distribution could be correlated across wide expanses of interstitial habitat and are not assumed to depend on the density of suitable habitat patches.
Estimated regression coefficients were similar across all models for both wetland perimeter-to-area ratio and wetland extent. For all models, wetland perimeter-to-area ratio had a significant negative effect on crane presence, and wetland extent had a significant positive effect (Table 1). However, there was a slight difference in estimated effect of density of surrounding urbanised areas between the IID and AR1 models. Although density of surrounding urbanised areas was found to have a negative effect in all models, this was estimated as not being significant in the binomial presence/absence and marked point process models which incorporated an AR1 temporal correlation structure, but was estimated as being significant in those models with an IID temporal structure ( Table 1).
Probability of presence of breeding pairs of cranes across space in 2015 was predicted using the binomial presence/absence and marked point process models with IID temporal structure.
Predictions from the binomial presence/absence model (Figure 2, left plot) show a wider spatial range in areas of high probability of presence, compared to predictions from the marked point process model (Figure 2, right plot).

TA B L E 1
Posterior mean and 95% credible intervals for: Regression coefficients of environmental covariates; scaling parameter ( ) representing the interaction between G(s) and the probability of crane presence; temporal correlation parameter from the AR1 process; parameters of the spatial, G(s) , and spatiotemporal, M(s, t), fields. The covariate coefficients and scaling parameter values with credible intervals, which do not cross zero are in bold. Model run times and Watanabe-Akaike information criterion (WAIC) scores are also given. The covariates included in the model were: Density of surrounding urbanised areas (within a 10 km terrestrial buffer), wetland perimeter-to-area ratio and wetland extent. The values of these covariates were standardised using the Z-score formula prior to inclusion in the model. All values given are rounded to 2 decimal places

| DISCUSS ION
In this paper, we have fitted four different models of varying complexity. The simplest model is a spatio-temporal model with a single likelihood with an IID assumption between years. We will now com- F I G U R E 1 Estimated point random field (G(s), left) and mark random field (M(s, t), centre) for 2015 from the marked point process model with IID temporal structure. Colour scale is given in low-high intensity as we are interested in relative differences across space and not absolute values. Distribution of wetland locations and observed crane presences in England in 2015 (right). Purple background colour is used for visual clarity and does not represent a value on the colour scale.

F I G U R E 2
Mean predicted probability of the presence of a breeding pair of cranes at each wetland location in 2015. Predictions were made using the binomial presence/absence model (left) and the marked point process model (right) with IID temporal structure. Colour scale is given in low-high probability as we are interested in relative differences across space and not absolute values. 'Low' represents a value of mean predicted probability of crane presence less than or equal to the 99th quantile of mean predicted probabilities across all four models analysed. 'High' represents a value of mean predicted probability of crane presence greater than this cutoff point. Purple background colour is used for visual clarity, and does not represent a value on the colour scale.
The scaling parameter represents the strength and direction of the interaction between the point field G(s) and the probability of crane presence. In the marked point process model with an AR1 temporal structure, the significant positive effect of wetland density on the probability of crane presence is of greater magnitude than that of the model with IID temporal structure ( Table 1). This is due to the fact that the AR1 model borrows strength across years, meaning that the repeated presence of cranes in high-wetland-density areas over time indicates a stronger effect than that observed in the IID model, which considers each year of data independently.
The spatiotemporal field M(s, t) in the AR1 models has a smaller posterior spatial range and larger standard deviation than that of the IID models ( Table 1). This is due to the fact that the AR1 models borrow strength across years, meaning that the spatial structure is influenced by the repeated occurrence of relatively isolated crane presences over time. Similarly, the standard deviation is high, as this allows sufficient variation in the field for the estimated species intensity to range from low to very high values. On the other hand, the IID models consider each year of data independently, and so are more strongly influenced by the priors, which suggest a larger spatial range and lower standard deviation, with the intensity of the species correlated across larger intervals and lowering smoothly away from areas of higher intensity. While the priors are biologically realistic, they may contrast with the trends observed in the data due to the early stage of reintroduction at which the data was collected.
In the models with AR1 temporal structure, the temporal correlation, t , is estimated to be extremely high (Table 1) and the spatial structure of M(s, t) is predicted to be the same for each year of data (Appendix A, Figure A2). This indicates that the spatiotemporal field in these models captures a spatial structure in the data that is constant over time. This strong trend of a static spatial structure may be due to the early stage of reintroduction when observations were made. At this early stage, the population is small and so may not be subject to strong dispersal pressures such as habitat density dependence, meaning that movement between habitat patches may be less common than it would be in a densely populated area. The population is also only observed in a small proportion of its potential range, meaning that a large amount of the study area is taken up by consistent zero counts. Limited movement, large numbers of consistent absences, and the binary nature of presence/absence data all lead to few observed changes in the dataset between years (Appendix A, Figure A1). Therefore, the static spatial structure picked up by the model may not reflect the behaviour of the study species, for example, strong site fidelity, and instead could be a feature of the stage of the reintroduction process the population was in when the data were collected. A longer period of observation may be required to accurately capture and predict the dispersal dynamics of the species.
Since the AR1 models detect the strong trend of a static spatial structure over time, they estimate no change in the spatial representation of the spatio-temporal field M(s, t) between years. This means that the spatio-temporal field must represent the most likely spatial correlation structure for all years, and so produces an average from the full dataset. This averaging can produce misleading estimates of species distribution, particularly when the species is observed to 'appear' in a new area of the study region part-way through the study, as it can be predicted to have a high probability of presence in this area before the appearance event. For example, breeding pairs of cranes were first observed in Somerset (South West) in 2013 as the result of a successful reintroduction project (Appendix A, Figure A1).
However, the spatio-temporal field from the AR1 models shows a high intensity of cranes in this area in the years prior to this event, when no breeding pairs of cranes had been observed in the area (Appendix A, Figure A2). In this example, the covariates included in the model are also stationary across the interval of observation.
As such, predictions from the models with AR1 temporal structure are very similar across all years in the study period (Appendix A, Figure A5). This not only conceals fine scale interannual differences in spatial distribution but also removes novel insight from future predictions. Such an effect may be avoided by including an extra, purely spatial, field in the mark likelihood and setting priors to manually restrict this to capturing the static spatial effect, separating this strong signal from the changes in spatial structure over time which could then be picked up by the spatiotemporal field. As mentioned above, the effect observed here is a result of the data structure, and so a larger dataset collected over a longer duration may provide sufficient information to avoid these issues.
The issues created by the strong, static spatial trend picked up by the spatiotemporal field in this example demonstrate that inclusion of complexity (here, an effect accounting for temporal correlation) may not necessarily improve ecological insight gained from a model. An even more complex model, containing another purely spatial field and strong priors, may improve output quality when aiming to predict future distribution of the population.
Of the four models considered in this paper, the most complex (the marked point process model with AR1 temporal structure) had the longest running time, and the simplest (the binomial presence/absence model with IID temporal structure) had the shortest (Table 1).
While running time is an important methodological consideration, each of these models is made extremely computationally efficient through the use of the INLA approach, allowing model choice to be more strongly determined by ecological insight gained, as opposed to computational efficiency.
Watanabe-Akaike information criterion (WAIC) scores were computed for each of the four models (Table 1), although it is important to note that the binomial presence/absence and marked point process models cannot be compared using this criterion.
However, within each model type, it is interesting to note that the AR1 models have a lower score than the IID models, indicating a better fit to the data. This is unsurprising when considering that the AR1 models are more strongly influenced by the data, as information is borrowed across years, whereas the IID models are more reliant on the priors. The true meaning of such criteria in the area of point process modelling is not well understood, due to the historically theoretical nature of the topic, and as demonstrated in this manuscript, these criteria do not necessarily provide a definitive answer as to selecting the 'best' model for interpretation and ecological insight.
Species distribution data is often a combination of the true underlying distribution of the species, and the observation process used to collect data, creating a need to disentangle the true drivers of species distribution from sampling effects. In the early stages of a species invasion or reintroduction, there is a time lag, which means that the observed distribution may not match the true potential distribution of the species once the population has been established (De Marco et al., 2008). Making observations at this initial stage could bias inference of the effects of environmental covariates, as the distribution of the species is restricted by dispersal limitations, which remain unaccounted for. Including complex model components to represent spatial correlation structures, as is demonstrated here via the inclusion of G(s) in the marked point process model, can account for the fact that the species distribution is limited by the early stage of its introduction to the environment, giving a more accurate prediction of the range being constrained close to where the species has been observed. Accounting for spatial correlation can also aid in avoiding Type I errors when inferring the significance of covariate effects . However, caution is advised when incorporating additional structural complexity into models of this type of data, due to its zero-inflated and relatively static nature (observed here through the issues associated with inclusion of an AR1 temporal correlation structure). Accounting for the observation process, and spatially varying detection probability is a particular strength of the inlabru package, which can also be used to model subsampled and distance sampling data (Jullum et al., 2020;Yuan et al., 2017).
Complex model components such as spatial and spatio-temporal random fields can be used to represent population dynamics in species distribution models. In this example, density of wetland

ACK N OWLED G EM ENTS
The authors thank Andrew Seaton (University of Glasgow) for assistance in the development of Appendix B on PC priors.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/2041-210X.13957.

DATA AVA I L A B I L I T Y S TAT E M E N T
A transformed version of the data used in this analysis is available on Dryad Digital Repository https://doi.org/10.5061/dryad.2z34t mpps (Laxton et al., 2022b). These data have been randomly transformed to prevent an exact identification of nest locations and avoid potential disturbance to cranes during the breeding period. R code which can be used to fit the models developed in this manuscript to