Estimating the abundance of a group‐living species using multi‐latent spatial models

Statistical models use observations of animals to make inferences about the abundance and distribution of species. However, the spatial distribution of animals is a complex function of many factors, including landscape and environmental features, and intra‐ and interspecific interactions. Modelling approaches often have to make significant simplifying assumptions about these factors, which can result in poor model performance and inaccurate predictions. Here, we explore the implications of complex spatial structure for modelling the abundance of the Serengeti wildebeest, a gregarious migratory species. The social behaviour of wildebeest leads to a highly aggregated distribution, and we examine the consequences of omitting this spatial complexity when modelling species abundance. To account for this distribution, we introduce a multi‐latent framework that uses two random fields to capture the clustered distribution of wildebeest. Our results show that simplifying assumptions that are often made in spatial models can dramatically impair performance. However, by allowing for mixtures of spatial models accurate predictions can be made. Furthermore, there can be a non‐monotonic relationship between model complexity and model performance; complex, flexible models that rely on unfounded assumptions can potentially make highly inaccurate predictions, whereas simpler more traditional approaches involve fewer assumptions and are less sensitive to these issues. We demonstrate how to develop flexible spatial models that can accommodate the complex processes driving animal distributions. Our findings highlight the importance of robust model checking protocols, and we illustrate how realistic assumptions can be incorporated into models using random fields.


| INTRODUC TI ON
Accurately and efficiently estimating the abundance and distribution of organisms is one of the most fundamental components of ecological research. For many species access to resources is becoming increasingly constrained either because contiguous habitats are becoming fragmented (Fischer & Lindenmayer, 2007;Løvschal et al., 2017), or the overall availability of resources is declining due to human use (Ogutu, 2016;Tilman, 2017). In addition, populations of many commercially harvested species are being exploited to the point of extirpation (Barnosky, 2011), and in some cases sudden outbreaks of disease cause massive declines (Milner-Gulland, 2015;Roelke-Parker et al., 1996). The combination of these pressures raises serious concerns about the long-term viability of many species and underscores the importance of getting rapid, accurate estimates of abundance. In most cases, practitioners and managers want to know not only the abundance of animals but also the drivers of long-term changes (such as the impacts of climate effects, over-exploitation or displacement) and also identify bottlenecks at specific life-history stages. This requires repeated estimates of abundance over time.
Ecological surveys are the primary tool for estimating abundance.
Although the design of a survey differs depending on the species (Hedges & O'Brien, 2012), in many cases the abundance is calculated by estimating the number of animals per unit area from a series of sampling points and inferring the density over the entire survey area (Norton-Griffiths, 1973). Samples are generally collected along a systematic pre-selected grid or transects that ensure maximal coverage and appropriate statistical power. Therefore, sampling requires adhering to a strict protocol; however, the practicalities of achieving this in the field can be very challenging (for example, transient obstacles such as low cloud-cover can force breaks in aerial transects and therefore sampling becomes discontinuous). Furthermore, observers may see groups of animals occurring off-transect which cannot be included in the enumeration, even if the observations are important. In many circumstances, a survey is constrained to a single species or a group of species that occur sympatrically resulting in multiple surveys required in areas that are rich in biodiversity.
The difficulties of estimating species abundance are compounded when organisms are highly aggregated. Many species occur in groups either because individuals distribute themselves based on the availability of resources in a patchy environment or because social interactions attract individuals to conspecifics (Krause & Ruxton, 2002;Westley et al., 2018). In addition, animals whose distributions are highly dynamic and subject to rapid fission-fusion processes (such as migratory species like wildebeest) pose a particular problem for aerial surveys because not only do they form fluid yet highly aggregated groups, the time window for conducting a population estimate is often very brief. This restricts the window of opportunity and limits the number of samples that can be feasibly collected to estimate the population size (i.e. the number of point observation, photos or transect). Furthermore, if the window is missed the survey cannot be repeated until the following year because animals have moved or dispersed, which places pressure on enumerators to be efficient and responsive.
The requirements of statistical models are therefore severalfold. Primarily, statistical methods must accurately infer population abundance from incomplete observations even when spatial distributions are highly hetereogeneous and subject to high degrees of autocorrelation. Where possible models must be sufficiently flexible to account for irregular survey design and the realities of the data collection procedure, thereby taking pressure off ecologists and wildlife managers in the field. Models must also be complex enough to make full use of the available data and maximize the accuracy of estimates while minimizing their uncertainty. Finally, models must be practical; fitting models to relatively large datasets should be achievable within reasonable time-scales and ideally implemented within well-supported open source software libraries and packages.
Surveys of animal aggregations often have strong spatial autocorrelation between observations. One approach to reduce the effect of spatial autocorrelation is to use the samples to calculate the mean density per transect or stratum and use this summary to calculate abundance, as in the Jolly II method (Jolly, 1969). In this case, all observations within a transect are assumed to be dependent (i.e. autocorrelated), but transects are assumed to be independent. The sample size is defined by number of transects (not the total number of observations across all transects), which can reduce the statistical power.
Recent developments in this area have focused on the use of latent Gaussian random fields (GRFs; Hristopulos, 2020;Rue & Held, 2005). Random fields are a major component of spatial statistical analysis and can be used to represent spatial and spatiotemporal correlation structures which remain unexplained by environmental variables. In this context, a likelihood is specified for an observed response variable and the random field defines the mean, or other parameter, of the likelihood function. A GRF is defined by a mean and a covariance function and has the property that any finite collection of points follows a multivariate Gaussian distribution with the specified mean and covariance. Hence, the covariance function determines the similarity between two locations of the field and is often assumed to be stationary and isotropic, meaning the similarity importance of robust model checking protocols, and we illustrate how realistic assumptions can be incorporated into models using random fields.

K E Y W O R D S
aerial surveys, animal abundance, multi-latent fields, spatial modelling between two locations depends only on their distance, and not their relative direction or absolute position. Latent Gaussian fields are explicitly or implicitly used in many spatial analysis methods, including INLA (Rue et al., 2009), kriging (Cressie, 1990), Gaussian process (GP) regression (Rasmussen & Williams, 2006), and spatial GAMs with smoothing splines (Wood, 2020). In all these implementations a key challenge is appropriately specifying the smoothness properties of the spatial field in a manner that accurately captures the autocorrelation of the data without overfitting to the noise.
Accounting for spatial structuring due to social behaviour, density dependence, or other effects unexplained by environmental variables has been found to improve certainty of abundance estimation. This is of particular importance in management and conservation of keystone species, and in estimating quotas for stock sustainability in fisheries (Calderon-Aguilera et al., 2019). GRFs have also been shown to be effective in dealing with data containing excess zeros when incorporated with zero-inflated or hurdle models (Welsh et al., 1996). This approach can be implemented using multiple random fields where one field defines the probability of the presence or absence, and the second field determines abundance conditional on the species being present (Barry & Welsh, 2002;Sadykova et al., 2017;Smith et al., 2019).
In this work, we explore a scenario in which spatial models are used to improve the accuracy of abundance estimates for an aggregated species. Our approach uses multiple random fields to estimate the abundance of wildebeest in the Serengeti National Park from aerial survey data. We specifically investigate the importance of model complexity and how it interacts with the spatial patterns that arise due to the ecological and behavioural dynamics of the wildebeest .
Wildebeest are a keystone species in the Serengeti. Each year, approximately 1.3 million wildebeest move through the range of the Greater Serengeti Ecosystem, a round-trip of >1,500 km (Hopcraft et al., 2015;Thirgood et al., 2004) that covers almost the entire 25,000 km 2 extent of the ecosystem. During the migration, herds move constantly and vary in size, from a few individuals to hundreds of thousands of animals (Hopcraft et al., 2014). Without the annual migratory cycle, much of the region's biodiversity would likely change (Dobson et al., 2010;Hopcraft et al., 2015), as the passage of the migrants affects virtually every other ecological interaction in the ecosystem from below ground nutrient cycling, to insects and avifauna abundance, to predator-prey interactions, to disease regulation .
Despite their large numbers, wildebeest operate at or close to their physiological limits. Historically, drastic population declines have been observed due to drought, access to forage and disease (Mduma et al., 1999). In adult females, the energetic demands of reproduction are continual and acute as the cycle of pregnancy and lactation encompasses a full calendar year; calving is highly synchronous and coincides with the period of peak primary production in February (Hopcraft et al., 2014). Approximately 400,000 calves are born each year of which about 50% survive to reproductive age (Hopcraft et al., 2015). The spatial diversity of the ecosystem means that at each point in the migration the quality and abundance of forage changes, as does the exposure to predators, as well as the animals' proximity to human infrastructure and poaching. As many as 100,000 wildebeest are illegally harvested every year (Rentsch & Packer, 2014), and there are upwards of 3,500 tourists per day entering the Serengeti. These anthropogenic inputs create additional hazards for animals to navigate as they strive to gain sufficient energy to survive and reproduce whilst avoiding predation. Any changes to the mass migration are likely to have severe environmental consequences as the wildebeest play a vital role in the ecosystem.
Estimating the abundance of wildebeest in the ecosystem is vital for the effective management of the park and is perhaps the most important metric of the ecosystem's health (Estes, 2014).
To estimate the abundance of the wildebeest population, transects are flown over the herds in March, April or May (Campbell & Borner, 1995;Norton-Griffiths, 1973)  To estimate population abundance from image counts, Jolly's method II (Jolly, 1969) is used. This provides an expected abundance calculated by multiplying the average within image density by the total survey area and a standard error calculated based on the standard deviation of the densities across each transect (the area of each image is calculated using the field of view of the lens and the altitude above-ground). It should be noted that the photographed area is typically less than 3% of the total survey area and typically between 5,000 and 10,000 wildebeest are manually identified out of a population total of around 1.3 million.
Here we compare the Jolly II method used to infer total wildebeest abundance with two hierarchical models that incorporate random spatial fields. We implement two spatial models; the first uses a single field to model the log intensity of a Poisson distribution for wildebeest abundance, while the second introduces a multi-latent model that makes use of a spatial mixture model with multiple random fields to account for the large regions of the survey area that contain zeros in the image counts. To evaluate the performance of these models we generate synthetic wildebeest distributions and vary the fraction of the landscape (denoted p) where wildebeest are present. Through our analysis of inferred wildebeest abundance we show that there is a non-monotonic relationship between model complexity and accuracy, as large errors are observed using a complex model when model assumptions are not met.

| Single-field model
Our first approach to incorporating a spatial field uses a log Cox GP.
The distribution of wildebeest is therefore, modelled as an inhomogeneous Poisson process with spatially varying intensity (x). A GP is used to model the log of the spatially varying intensity so that the full model for wildebeest density D w is defined as where GP defines a GP prior, is the mean of the log intensity field and K x, x ′ is a stationary covariance function that determines the similarity of two locations based on their distance.
Several methods are available for parameter inference and model fitting, including MCMC, variational inference and Integrated Nested Laplace Approximation (INLA; Rue et al., 2009). The appropriateness of different approaches will depend on the choice of covariance function and dataset size. Here we assume a Matérn covariance function that allows us to take advantage of the stochastic partial differential equation (SPDE) approach (Lindgren et al., 2011).
The solution to a given GRF with Matérn covariance is an SPDE, which in this approach is used to estimate a corresponding Gaussian Markov random field (GMRF). The GMRF can then be used to represent the GRF in computation, which due to its Markov property, considerably reduces computational cost, enabling complex spatial and spatiotemporal modelling within an efficient timeframe. Inference is implemented within inlabru (Bachl et al., 2019).

| Multi-latent model
A notable characteristic of the wildebeest distribution is that they are highly aggregated in certain regions of the survey area. Due to the survey design and the constraints of flying transects, this leads to large regions that contain no wildebeest. Counts of wildebeest with aerial photos, therefore, contain large numbers of zeros, which are well known to cause problems when fitting to Poisson distributions (Harrison, 2014).
A standard approach to dealing with zero-inflated data is to use a mixture model that combines a Bernoulli distribution with a Poisson or negative binomial distribution; in this approach the count is zero with a fixed probability, otherwise it is drawn from the discrete probability distribution that may be truncated to exclude zero (Wenger & Freeman, 2008). Although models of this type are better able to cope with zero inflated data, they do not account for the spatial autocorrelation in the presence of zeros. In the context of the wildebeest distribution, zero counts are highly correlated with one another and there are large clusters of dense and empty regions. It is unlikely these clusters are driven by environmental factors, instead they are driven by the social behaviour of the wildebeest and their aggregation into herds of various sizes.
To account for the spatial clustering of wildebeest we use a multilatent random field model (Saul et al., 2016). The model consists of two spatial fields, both modelled as latent GPs. The first field determines whether wildebeest are present in the local area or not, whereas the second field is used to model the log intensity of the wildebeest distribution, as in the single-field approach. Hence, the spatial distribution of wildebeest is determined by two independent random fields acting as a mixture as in a hurdle or zero-inflated model; one that determines presence or absence, and one that determines density given wildebeest are present. The model is then defined as follows: where, as in the single-field model, x is the spatial location, D w (x) is the density of wildebeest, which is modelled as an inhomogeneous Poisson process with intensity (x), and F(x) is a GP that is used to model the log of the spatially varying intensity. G(x) is a second, independent GP that is used to model the binary variable (x) which determines whether wildebeest are present. This variable is drawn from a Binomial distribution where G(x) determines the logit of the presence probability.
Given the additional complexity of the model, computational inference is more challenging. We again use an approximate method, this time using variational inference (Blei et al., 2017) to estimate the posterior distribution of the latent spatial fields. Variational inference follows a similar approach to the Laplace approximation used by INLA; however, it proceeds by firstly assuming a specific form for the posterior density and then minimizing the Kullback-Leibler divergence between the posited density and the target density.
Variational inference is a popular method within the machine learning community for fitting GP models (Hensman et al., 2013;Opper & Archambeau, 2009;Titsias, 2009) and has recently been applied to multi-latent GPs such as the model defined by Equation 2 (Saul et al., 2016). We follow this method (Torney, 2022) and use an implementation provided by the package GPFlow (Matthews et al., 2017).

| Synthetic data
To test the various inference schemes we generate synthetic wildebeest distributions and then simulate the survey by making (1) observations of the generated distribution within small areas along a sequence of transects. This process generates a set of counts corresponding to the aerial photographs collected during the real survey.
To generate the synthetic distributions that match the clustering of the real wildebeest, we generate two random fields. The first field is used to define a region where wildebeest are potentially present. We define this field to have a correlation length of 5 km and specify a fraction p of the survey area that may contain wildebeest.
This fraction is varied from p = 0.1 (i.e. 10% of the area may con-

| Empirical data
We further assess the multi-latent field approach using survey data A Cessna C182 aircraft was used to conduct 10.3 hr of photographic sampling flights along east-west transects that covered a straight-line distance of 2,040 km. Photographs were taken using a NIKON D800 through a 35-mm Nikor Lens. The camera was mounted in a port in the floor of the aircraft and an automatic trigger was used to collect images at 10-s intervals during each transect.
In total 1,584 georeferenced images were taken with a resolution of 7,360 × 4,912 pixels, which were later viewed and the number of wildebeest in each image was recorded.  The multi-latent approach displays a similar accuracy to the Jolly II method; however, it is less biased than Jolly II, which systematically appears to result in an undercount, and it is better able to quantify the uncertainty in the abundance. For all values of the fraction of the landscape potentially occupied by wildebeest, the true value was more likely to fall within the 95% credible intervals when compared to the 95% confidence interval of the Jolly method. The model was also able to accurately capture the latent presence-absence field used to generate the data as shown in Figure S4. It is likely that the uncertainty quantification of the multi-latent model could be further improved by careful selection of the covariance function used to encode spatial autocorrelation in each field. For the simulation study a Matérn 3/2 kernel was used for both fields but other functions that better approximate observed spatial patterns would likely result in greater accuracy.

| Empirical data
We next apply the multi-latent model to the 2015 wildebeest image counts. The spatial distribution of image counts is shown in Figure 3a. In this year, wildebeest were distributed across the landscape into several clusters with large regions of zero counts. After optimizing the model, the learned spatial fields define the probability of wildebeest being present, and their log-density given wildebeest are present. Figure 3b shows the presence probability for the 2015 data, whereas Figure 3c displays the log-intensity.
The optimized model parameters for the lengthscale of the presence-absence field and the log-intensity field are shown in Table 1 along with the population estimates for the year obtained using the multi-latent model and, for comparison, Jolly II. The posterior distribution of total wildebeest abundance is show in Figure 4 with lines indicating the posterior median and the upper and lower 95% credible intervals. Again, the Jolly II outputs are shown for comparison. Based on the results from the simulation study, we expect that the true abundance is underestimated by Jolly II and the credible intervals obtained from the multi-latent model are a more accurate reflection of the uncertainty in the estimates.

| DISCUSS ION
Our a hierarchical approach to spatial autocorrelation in which concentrated aggregations can occur within groups that are distinct from each other). This results in a comparable level of accuracy in abundance estimates with the Jolly II method. Furthermore, the multilatent model suffers from less bias than Jolly II and provides a more accurate quantification of uncertainty. Therefore, the inclusion of complexity in a model does not necessarily equate to improved performance, but increasingly complex models that make reasonable assumptions can provide advantages over simpler, more traditional approaches.
An extension to multi-latent field framework is to incorporate the effects of environmental features (see Figure S1). The use of environmental information allows the model to infer the relationship between known covariates and the species distribution and then use this relationship to make more informed predictions about density in unsampled regions. This approach is straightforward to implement via a modification to the GRF but could not be introduced into Jolly II due to its use of transect as the unit of observation. distribution, enabling the collection of data on multiple species in a single survey. In combination with deep-learning approaches to object recognition in images (Torney et al., 2019), this could be an extremely powerful approach, and potentially represents a significant advance in improving the efficiency of resource-use in counting animals for abundance estimation.
Our results demonstrate the traditional method for inferring wildebeest abundance, Jolly II, is a reliable approach that is relatively straightforward to implement. Given a fixed survey design that adheres to the requirements of the Jolly II method there is no clear advantage to implementing the more complex multi-latent field approach. However, the more complex model is both flexible and extensible. Potential extensions include the use of environmental covariates as described above, as well as explicitly modelling the temporal dynamics of the wildebeest herds (Blangiardo et al., 2013). The greater flexibility of the model translates to greater flexibility in the field when conducting the survey, as well as the prospect of including multiple species in the analysis and incorporating data from multiple sources (Graves et al., 2022;Robinson et al., 2021). Values shown indicate the expected log density of wildebeest. Areas where wildebeest are absent with probability >0.9 (as determined by (b)) are omitted.  We are grateful to two reviewers whose helpful comments and suggestions improved the manuscript.

CO N FLI C T O F I NTE R E S T
We declare that we have no competing interests.

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.13941.