Comparing distribution of harbour porpoise using generalized additive models and hierarchical Bayesian models with integrated nested laplace approximation

Species Distribution Models (SDMs) are used regularly to develop management strategies, but many modelling methods ignore the spatial nature of data. To address this, we compared fine-scale spatial distribution predictions of harbour porpoise ( Phocoena phocoena ) using empirical aerial-video-survey data collected along the east coast of Scotland in August and September 2010 and 2014. Incorporating environmental covariates that cover habitat preferences and prey proxies, we used a traditional (and commonly implemented) Generalized Additive Model (GAM)


Introduction
Understanding environmental and anthropogenic drivers of species distributions is critical for identification of potential areas for protection, development of appropriate conservation management strategies, and mitigation of negative anthropogenic impacts. Collection of empirical environmental data, however, can be risky, costly, and time consuming; consequently, desk-based species distribution modelling provides a cost-effective, yet successful alternative method to address these factors (e.g. Bailey and Thompson, 2010;Norberg et al., 2019).
Many statistical techniques are used to model species distributions, but due to inherent complexity of analysing space, most do not fully address the spatial nature of data. Generalised Additive Models (GAMs; Hastie and Tibshirani, 1990) are one of the most commonly used and well-established methods of performing species distribution modelling (Becker et al., 2020;Redfern et al., 2006;Wood, 2006). Over the past decade, there has been rapid development of specific statistical methodologies aiming to address the spatial nature of empirical data (e.g. Cameletti et al., 2012;Diggle et al., 2013;Illian et al., 2008;Law et al., 2009;Sadykova et al., 2017); however, such methods are complex and require a high level of statistical expertise and computational power. Practitioners are, therefore, often faced with a choice of which model to apply to their data, which can have considerable impacts on outputs (Norberg et al., 2019).
The development of point process models in ecological contexts (Soriano-Redondo et al., 2019;Yuan et al., 2017a), by using exact locations of points in space within a GAM function are being explored by some studies (Fithian and Hastie, 2013;Miller et al., 2019;Renner et al., 2015). These use an inhomogeneous Poisson process, maximum entropy, or infinitely weighted logistic regression to model data in a point process framework thereby taking advantage of the spatial information; however, at the time of the production of this study, these techniques were still under development for generalist use, and were therefore not explored further here.
The aim of this study was to assess differences between use of traditional frequentist GAMs and explicit spatial-modelling approaches using Hierarchical Bayesian Models (HBM) fitted with Integrated Nested Laplace Approximation (INLA, Rue et al., 2009). To achieve this, we introduce the two methods (GAM and HBM-INLA) and identify three areas in which Bayesian developments offer potential benefits over traditional GAMs: (1) they account for spatial autocorrelation, (2) clustering of animals, and (3) model in continuous space. We then implement both methods using an example of real-world aerial-video survey data for harbour porpoise (Phocoena phocoena). Finally, we explore differences in resulting spatial predictions of porpoise distribution and conservation management implications.

Generalized additive models
Generalized Additive Models are popular for species distribution modelling because of their ability to encapsulate non-linear interactions between sightings and environmental covariates through use of smoothing functions (e.g. Booth et al., 2013;Embling et al., 2010;Williamson et al., 2016); however, their ability to capture spatial trends is limited without use of alternative, more complex additions (Fithian and Hastie, 2013;Miller et al., 2019;Renner et al., 2015;Scott-Hayward et al., 2014). Nonetheless, GAMs are relatively easy to implement by the user through various R libraries, such as mgcv (Wood, 2011), or similar. Counts y in a grid cell S i , i=1,…N are modelled (response variable), and the model usually assumes that counts in each grid cell follow a Poisson distribution. The expected value in a grid cell s i is modelled as a linear model of covariates x, i.e., where, g is a log link function, x(s i ) T is a matrix of N observations on each of the k explanatory variables, and β is a vector of regression coefficients. The GAM formulation in Eq. (1) does not address the problem of spatial autocorrelation, instead accounting for trends in data across larger geographical distances (Dormann et al., 2007). One common method for circumventing this problem is to include a bivariate smoothing term of latitude and longitude or a soap-film smoother (Wood et al., 2008) in the GAM itself. Soap-film x,y smoothers also include boundary information to prevent the model from predicting animals' presence in areas they are not capable of occupying (Wood et al., 2008), e.g. porpoises present on land. An x,y smoother accounts for similarity of the response variable in adjacent cells, reducing the effect (but not the amount) of spatial autocorrelation. The spatial autocorrelation problem is solved by using more complex methods, such as including autocorrelation structures in a generalized additive mixed effects model (Carvalho et al., 2020).

Hierarchical Bayesian modelling using INLA
Recent developments in spatial statistical methodology have focused on producing a flexible, computationally efficient approach that is increasingly accessible to end users. Here, we used a continuously indexed Gauss-Markov random field based on the Stochastic Partial Differential Equation, SPDE (Lindgren et al., 2011). This provides an approximation to a continuous spatial field that allows users to work efficiently with complex spatially explicit models (Blangiardo et al., 2013). In addition, models were fitted based on INLA (Rue et al., 2009), a recently developed computationally efficient method for hierarchical Bayesian inference that is particularly suitable for complex spatial models, such as those based on the SPDE approach. Combining the two methods (SPDE and INLA) provides a suite of flexible models that allow for realistic spatial modelling within feasible computation timescales.
As in the GAM model, counts y in a grid cell s i , i=1,…,N are modelled, but a Bayesian approach is used and spatial autocorrelation is modelled by a random field ξ, i.e., where, g is a log link function, and β is a vector of regression coefficients. In this HBM formula, the residual spatial autocorrelation is modelled explicitly by including the ξ term, which specifies a spatial or spatiotemporal random field (Rue and Held, 2005). This framework allows for smooth and nonlinear effects of covariates, time trends, seasonal effects, random intercepts and slopes, and spatio-temporal random effects; therefore, this class of model is highly flexible and can accommodate a wide range of paradigms.

Discrete vs. continuous space
Spatial data are usually gridded or segmented when using GAMs, which causes spatial information to be lost; however, spatially explicit approaches such as SPDE/INLA can model in continuous space. Since this implies much less information loss, SPDE/INLA approaches were expected to be more appropriate and capable of reflecting small-scale spatial behaviour that cannot be captured by gridding or segmenting sightings. Consequently, in this paper we made two comparisons. We initially used both approaches, GAM and HBM, to model data on the same spatial grid to establish how the approaches differed when the actual models were similar, with the expectation of similar results. We then moved towards modelling in continuous space using the SPDE/ INLA approach that also allowed us to incorporate information on the exact location of transects, as well as the spatial structure formed by the locations of individual animals, i.e. their clustering (Yuan et al., 2017b).

Spatial clustering of animals
Animals do not always move independently. At times, some species can form large groups in small areas for several reasons, such as social cohesion or specific attraction to the location as a result of physical and/ or environmental covariates. Plants and animals are often clustered spatially according to a Poisson process, which inappropriately assumes complete spatial randomness between clusters. To account for this, a Cox process (Møller and Waagepetersen, 2007;Waagepetersen and Schweder, 2006) is a flexible type of spatial point procedure which improves accuracy by modelling spatial patterns relative to observed or unobserved spatial trends . In this way, the spatial pattern can be modelled along with traits of individual observations such as markings, scars, size, or colour and environmental covariates . Since clustering violates Poisson process assumptions, this obliges the SPDE to capture low-range peaks in detection, rather than a large-scale distribution trend. An alternative (and often implemented) approach is to use the negative binomial distribution to account for overdispersion which means that the SPDE does not have to accommodate spikes to capture the large-scale trend; however, while often used to avoid violation of model assumptions, it does not explicitly model overdispersion. Moreover, there is no obvious way of extending this to continuous space, so analysis is restricted to using gridded bins.
When using hierarchical Bayesian models, localised clustering can be differentiated from large-scale patterns using the Log-Gaussian Cox Process (LGCP) method which models exact locations of points (Diggle et al., 2013). The model assumes point locations are scattered independently in space as a Poisson process (top level of the hierarchy), given a model of the point density (formally the "intensity" of the processthe middle level of the hierarchy) and a set of priors for all model parameters (bottom level of the hierarchy). The log of intensity Λ in a location s ∈ R is modelled as a linear model of covariates x and a random field ξ that is again approximated by an SPDE:

Case study
Harbour porpoise is the most numerous cetacean species in European Atlantic Shelf waters (Hammond et al., 2021) which is also where most anthropogenic activities occur. Consequently, they are threatened by, inter alia, unsustainable incidental catches in fishing gear and prey removal, various forms of chemical and underwater noise pollution, and other anthropogenic pressures (e.g. Culik, 2004;Hammond et al., 2021;MacLeod et al., 2007). Knowledge of porpoise distribution in relation to environmental drivers is therefore critical for effective implementation of management and mitigation measures.
Depth, slope, bottom temperature, and sediment have been significant in previous models of porpoise distribution (e.g. Booth et al., 2013;Brookes et al., 2013;Sadykova et al., 2017;Williamson et al., 2016). Other variables included chlorophyll-a, net primary productivity, and fronts and tidal mixing (e.g. current speed, vertical shear, and potential energy anomaly to quantify stratification strength), which cause prey aggregation, and have been correlated previously with the distribution of fish and marine predators (e.g. Embling et al., 2012;Hofmann and Powell, 1998;Olson and Backus, 1985;Scott et al., 2010).
Previous distribution modelling in the Moray Firth, north-east Scotland (part of the current study area; Fig. 1) found a relatively high density of animals in the Smith Bank, a shallow (30-40 m) sand bank in the middle of the Firth (Brookes et al., 2013;Williamson et al., 2016). This is likely because the Smith Bank provides suitable habitat for sandeels (Ammodytidae) and whiting, Merlangius merlangus (Hopkins, 2011), two of the most common prey species for porpoise in this area (Santos and Pierce, 2003). Knowledge of porpoise distribution along the East coast of Scotland, however, is lacking. Consequently, the present analysis combines previously collected aerial-video-survey data with a newer dataset, which allowed us to investigate distribution along the entire east coast of Scotland.
Aerial video surveys use video cameras to record strip transects below the survey aircraft and are thus different to visual surveys by human observers (Buckland et al., 2012). Moreover, video surveys provide a permanent record of sightings, which can be beneficial for clarifying mis-or ambiguous identifications, anomalous results, improving data transparency, and audit; indeed, some regulatory authorities now require such a permanent record (BSH, 2013). Estimates of relative density of porpoise modelled from video-survey data have been found previously to correlate strongly with relative density estimated from visual surveys (Williamson et al., 2016).

Aims
The purpose of this study was to compare two modelling frameworks to investigate differences between results obtained from a spatialmodelling technique (HBM using INLA) and a GAM which is used commonly by ecologists. We illustrated differences between these methods in how they model spatial autocorrelation, groups of individuals, and discrete vs. continuous space. We performed this using a case study of aerial-video transect-survey data to analyse harbour porpoise distribution along the East Coast of Scotland and discussed the potential management implications of results.

Case study data collection
Aerial digital video surveys covering a total of 5762 km were performed by HiDef Aerial Surveying Ltd. in August and September 2010 and 2014 along the east coast of Scotland ( Fig. 1 and Table 1). Video data were processed by qualified, trained, and experienced observers, who extracted all non-avian objects for identification by specialists at WWT Consulting Ltd., resulting in a total of 303 individual porpoises observed.

Data preparation for gridded analysis
Once porpoise observations were extracted, a 5 × 5 km grid, matching previous marine mammal work (Jones et al., 2015), was created which covered the entire survey area. Transect lines were intersected to the grid and the total length of transects and number of porpoises sighted within each grid cell was calculated.
Depth on a raster grid of approximately 180 × 180 m and polygons of sediment type at a 1:250,000 scale were provided by SeaZone Solutions Ltd. (2005a;2005b). Sediment types were classified using a Folk triangle (Folk, 1954), and expressed as the proportion of sediment that was sand or gravelly sand within each grid cell, based on previous studies of porpoise habitat association in this area (Brookes et al., 2013;Williamson et al., 2016). Seabed slope was calculated in degrees using the Slope tool in ESRI ArcGIS 10.2.1 (ESRI, Redlands, California, USA).
Modelled environmental data were provided from runs of the NEMO-ERSEM 3D coupled hydrodynamic-ecosystem model (Nucleus for European Modelling of the Ocean; Edwards et al., 2012;Madec, 2008;O'Dea et al., 2012). Parameters available included: bottom temperature (⁰C), potential energy anomaly (J/m 3 ; which is the energy required to mix the water column completely), depth-averaged current speed (m/s), vertical velocity shear (m/s), depth-averaged vertical velocity (m/day), depth-integrated net primary production (tons/day), and maximum chlorophyll-a (mgC/m 3 ). These data were available for the "summer season" (July, August, September, and October), as climatological means across 25 years (1989-2014) at a 7 × 7 km spatial scale. An estimate of stratification was calculated using the log 10 of depth (H) divided by current speed (u) cubed (log 10 (H/u 3 ); Simpson et al. (1981)). In addition, satellite data on a 1 × 1 km spatial scale were obtained from NEODAAS (NERC Earth Observation Data Acquisition and Analysis Service) and processed to derive ocean-front metrics (Miller et al., 2015): distance to the closest major ocean front (km), and front sidewhether a location is on the cold (negative values) or warm (positive values) side of the closest major front, or directly on that front (zero). These satellite data were provided as weekly means during the survey period; however, due to cloud cover obscuring parts of the study area each week, the mean in August and September of 2010 and 2014 were calculated and joined to the survey data from the corresponding year. Where surveys overlapped between the two years, the mean of August and September in both years was used.
All environmental variables were up or down sampled (depending on spatial resolution) and joined to the 5 × 5 km grid for inclusion in the models. Variables were also standardised by subtracting the mean value. The above variables were selected not necessarily because they were thought to influence porpoise distribution directly, but often because they are thought to be drivers of the distribution of porpoise prey (e.g. Scott et al., 2010). After inspection of histograms of the data, cells with a vertical shear greater than 0.15 m/s and a front side of < − 2 were removed, as data outside these ranges were sparse.

Data preparation for LGCP analysis
The same 5 × 5 km grid of environmental data was used for this analysis; however, number of sightings in each grid cell, and length of survey effort in each cell were not joined to this grid. Instead, the exact location of sightings, as well as transect start and end points were used directly.

GAM analysis
The first model was created using GAM. Only cells in which surveys were performed (e.g. survey effort >0 km) were included during model development. A 10 × 10 grid of points spaced evenly throughout the study region in which points that did not overlap with land were used for creation of the soap-film smoother (resulting in 27 internal knots; Fig. S1). Smoothers with more (up to 30 × 30), and fewer (5 × 5) knots were also tested, but not found to influence results. Models were tested with the maximum likelihood method using families such as negative binomial, Poisson, and zero-inflated Poisson. The negative binomial likelihood was selected for analysis after inspection of summary plots of models. In every model, an offset of the log of the area surveyed was included. The area surveyed was calculated by multiplying the effort (length surveyed in each cell in km) by the strip width (width of sea surface in view of the camera in km).
All environmental variables were tested for collinearity with each other, and those with a correlation of 0.5 or greater were not included in the same models. After initial investigation, relationships between most covariates and sightings were found to be linear; therefore, for simplicity all covariates were treated as linear in the model except for the soap-film smoother. Backwards selection was used to select the best model, in which a full model was fitted which included all environmental covariates except the soap-film smoother (e.g. without the influence of space), and then each environmental covariate was removed sequentially and the model re-run. This process was repeated consecutively using the model which had the lowest Akaike Information Criterion (AIC; Akaike, 1973) until there was no-longer a decrease in AIC when further covariates were removed. An Auto Correlation Function (ACF) plot was used to check for temporal autocorrelation, which was not found to exist; therefore, use of mixed models was not investigated. The selected model was run with the soap-film smoother and then used to predict relative density of porpoise in all cells within 10 km of the survey transects (Fig. 2). A variogram was created to investigate if spatial autocorrelation was evident in the model (Diggle and Ribeiro, 2007), which was not found to be an issue (Fig. S5). All analysis was performed in R version 3.2.4 and GAMs were modelled using the mgcv package (Wood, 2011(Wood, , 2017. R code can be seen in Appendix 1 and is available on GitHub.

Hierarchical Bayesian analysis using INLA
Two models were created using hierarchical Bayesian methods fitted using INLA. Both models were fitted using the inlabru version 2.2.4.9000 and INLA 20.03.17 packages in R version 3.4.1 (Bachl et al., 2019;R Core Team, 2017;Rue et al., 2009). The first model used the same gridded data that were used in the GAM analysis with the addition of the SPDE to account for spatial autocorrelation (Lindgren et al., 2011). The SPDE approach uses a Matérn covariance function to approximate the Gaussian field using a flexible stochastic model which provides a continuous relationship between the response and explanatory variables throughout the study area. Similar to the knots in the soap-film smoother in GAM, the SPDE approach requires a triangulation mesh of the modelled area. The mesh provides a lower bound on the spatial resolution for analysis; therefore, a mesh should be developed which is fine enough so that no further changes in the results can be observed when a finer mesh is used (Lindgren et al., 2011). The triangulation mesh for the SPDE was created bounded by the coastline and the boundary of the study region and allowed to place vertices randomly as needed (Fig. S2). Default uninformative prior specification (Gaussian with mean of 0 and variance of 100) was used for all variables. The same effort offset that was used in the GAM was included here (log of transect length * strip width) and again, environmental covariates were treated as linear. A negative binomial likelihood was selected to address clustering of the animals and to be able to compare more closely with the GAM results. Model selection (excluding the SPDE) was performed in the same way as described for the GAM; however, a Widely Applicable Information Criterion (WAIC) -also known as Watanabe AIC (Watanabe, 2010) -was used instead of AIC. This selected model was run once more with the SPDE included and was then used to predict relative density throughout the study area (Fig. 2).
The second, continuous space model created using HBM used the point process likelihood (Simpson et al., 2016) which has been adapted to line transect surveys (Yuan et al., 2017b). This approach represents the locations of animal sightings as realisations of LGCP observed at the area covered by the transects. As opposed to the gridded approach, where data are aggregated, the LGCP allows the locations of observations and transects to be used directly. Taking the actual transect locations into account allows us to distinguish between areas where we did not find any animals because no survey was performed in that location, compared to those areas where we did not find anything, but they were also surveyed. This is especially useful when transects are not placed randomly throughout a study area. Information about the relative and absolute animal location is thereby retained, which allows for a higher precision in estimating the relationship between covariates and density, as well as the residual spatial correlation structure. The LGCP model was created using the same set of covariates and the SPDE described above. Similarly, backwards selection using WAIC was used based on the environmental covariates without the influence of space and once the best environmental variables were identified the SPDE was included, and relative density was predicted throughout the study area. R code used to create the HBM using INLA and the LGCP can be seen in Appendices 2 and 3.
Both HBMs were also predicted onto the same grid used in the GAM models to facilitate comparison between methods. Spearman's correlation coefficients were calculated between GAM and HBM results. Moreover, for visual comparison, results of each model were divided by their maximum value (to standardise) and HBM results were then subtracted from GAM results to investigate if GAM or HBMs predicted higher relative densities (Fig. 4).

Results
The negative binomial likelihood was selected for GAM analysis after inspection of summary plots of the models (Fig. S3). A negative binomial likelihood was also selected for the HBMs using INLA on gridded data for comparison to the GAM. The best models selected using each of the three methods included many of the same covariates (Table 2); however, potential energy anomaly, net primary productivity, vertical shear, and front distance were not retained in any of the best models. This best GAM model had a soap-film smooth term that included latitude, longitude, and the boundaries of the study region to account for spatial autocorrelation. Spatial correlation was not included as a random effect in this model because a variogram of model results did not show any relationship with distance (Fig. S5). This model predicted the highest relative density of porpoise on the Smith Bank, and an area of high density along the southern coast of the Moray Firth. This model explained 33.2% of deviance ( Fig. 2 and S4).
HBM models used the SPDE approach to account for spatial autocorrelation. When comparing the HBM of gridded data with the GAM, the same areas of highest density were apparent with an additional third area of high density along the east coast of Scotland (Fig. 2). In addition, much finer-scale clustering was apparent in the HBM results. The predicted distribution of porpoise from the GAM showed broad areas of high density. This contrasted with the HBM on a grid, particularly in the southern half of the study area where the area of high density was much more restricted than in the GAM. Comparing the HBM using LGCP (Fig. 2) with the GAM, we again saw much finer localised clustering,  with only the Smith Bank highlighted as a high-density area. This suggested that the model was identifying areas where animals were either found in groups, or they used very fine-scale areas (1-2 km), which were not apparent in the gridded GAM analysis. Maximum coefficients of variation for the models fitted using HBMs were much lower than those created using GAM (GAM CV = 33, HBMgridded CV = 8, HBM -LGCP CV = 11, Fig. 3), which showed higher confidence in the models. The influence of environmental covariates on models can be seen in Figs S4, S7, and S8. Effects of environmental covariates were similar between the three models, with the sign of each variable being the same even though the exact numbers varied ( Table 2).
The GAM also had quite large areas with high Coefficient of Variation (CV), whereas both HBMs had a maximum CV less than half that of the GAM, but also the highest areas of CV in the HBMs were very localised and only at the edges of the study area in areas with very few data (Fig. 3). Correlation with GAM results was higher for the HBM of gridded data (Spearman's ρ = 0.76) than for the HBM using LGCP (Spearman's ρ = 0.68). The GAM consistently predicted higher relative porpoise density just south of the Smith Bank and both HBMs predicted higher density on the Smith Bank and along the east coast (Fig. 4).
Scatterplots were created to investigate how well predicted relative densities matched observed data (number of sightings divided by the survey effort; Fig. S6). Correlation between these ranged between 0.09 and 0.51.

Discussion
All three models had similar predicted distributions of porpoise (Fig. 2), with important differences. Firstly, compared to GAMS, HBMs (particularly when using LGCP), revealed more fine-scale areas of porpoise higher relative density, with an additional fine-scale (<5 km) area, that was not identified in GAMs. Secondly, compared to GAMS, HBM LGCP predicted a more realistic patchier porpoise distribution. This trend was not apparent in the other two models, which used the negative binomial distribution and gridded data, which smoothed out the highest-density areas. Finally, the best models from each of these three techniques included tidal mixing and chlorophyll-a, while two models also included depth, slope, and vertical shear ( Table 2), suggesting that these may be key drivers of porpoise distribution (Brookes et al., 2013;Embling et al., 2012;Sadykova et al., 2017;Scott et al., 2010;Williamson et al., 2016).
In terms of relating these model outputs to real porpoise ecology, animals in the field tended to be observed singly or in small groups of 2-3 animals. Clustering observed in models can likely be explained in terms of animal use of fine-scale environmental features that influenced prey availability (Table 2; e.g. Scott et al., 2010;Williamson et al., 2017), wider social factorssuch as cooperative hunting behaviour (Ortiz et al., 2021) or simply porpoise group spacing on survey sampling days. In terms of oceanography, thermal fronts are fine-scale, transient features which aggregate fish (Hofmann and Powell, 1998;Olson and Backus, 1985), and consequently have been related significantly to porpoise and other cetacean detections (e.g. Mendes et al., 2002;Philpott, 2013). Finally, the scale at which modelling was performed may have Influenced the impact/importance of variables between models (Mod et al., 2020).

Variability
Compared to GAM, variability for both HBMs was less than half, showing twice the level of certainty in the results. Some values of high CV in the GAM were in areas where there were survey data, whereas in the HBMs, areas of highest CV were restricted to the perimeter of the study area, where there was no survey effort. The GAM predicted much larger areas (25+ km) of high density than HBMs because the GAM smoothed out some spatial trends by gridding the data and using the negative binomial distribution, whereas the HBMs took the spatiality explicit autocorrelation into account.

Discrete vs. continuous space
The continuous modelled data using the SPDE/INLA, offered clear advantages because it eliminated information loss associated with data gridding, as performed in traditional GAMs. This is expected to be similar when segmenting data and can result in more accurate predictions of distribution. Both the HBM of gridded data (analysed in continuous space), and the HBM using LGCP (which used exact locations of sightings and transects) have smaller areas of high density (5 km) with a patchier distribution. L.D. Williamson et al. In general, the continuous representation and triangulation using the mesh provided a more flexible way of approximating a spatial area. This becomes relevant both when the observation area is complex (with a lack of data, etc.), and if the study area represents a large area of the earth's surface and a fixed grid distorts special relationships. In this case, models can be defined on the sphere without the need to project into two-dimensional space and the resulting distortion, which is inevitable when working with a grid.

Spatial clustering
Results of the HBM using LGCP showed a slightly different pattern to either of the other models fitted to the gridded data. This LGCP model showed more localised areas of high density, with overall low density everywhere else. Again, this was because animals were observed in clusters. Using the negative binomial likelihood (as in both the GAM and HBM using gridded data) is an arbitrary fix for this grouping problem, which smooths data within the grid, but there is not a straightforward interpretation of how the negative binomial likelihood reflects overdispersion. Consequently, if appropriate for the data and species being studied, using Bayesian modelling with the LGCP provides the best advantage because it does not smooth out patchy areas; however, porpoise is a highly mobile species, and performing the same survey on a different day would likely give different results, especially because the animals' distribution is related strongly to environmental covariates, which can change on timescales of minutes. With repeated surveys, the LGCP approach offers potential to differentiate between fine-scale clustering caused by grouping of animals, and large-scale clustering caused by the underlying distribution pattern, but it was not possible to investigate this here due to a lack of data from repeat surveys. Nonetheless, the exercise based on survey data available, produced a much clearer and more flexible approach than using the negative binomial likelihood, but was not as straightforward to implement when using GAMs. Without several repeated surveys, it is unclear if these very finescale trends reflect the true distribution of porpoise. A single survey of this type might be appropriate for sessile organisms, but this appears to be modelling at a finer resolution than is realistic for data from a single survey of mobile organisms.
Relatively low correlation was observed between predicted relative densities and observed data (Fig. S6). A potential reason for this arises from a spatial mismatch between where surveys or sightings fell within the 5 × 5 km grid cell and the centre of the cell where predictions were calculated. Also, models considered information from neighbouring cells when calculating predictions, which could be increasing/ decreasing predictions and would not be reflected by survey data within a particular cell. In terms of relating this to empirical observations of porpoise spatial clustering of animals in the field, the situation is more complex. For example, in energetically active tidal flows, there are linkages between porpoise presence and small-scale heterogeneity amongst different environmental covariates (e.g., tidal phase, time of day) and significant spatiotemporal variability in site use at scales of hundreds of metres and hours (Benjamins et al., 2017).

Future developments
Inlabru (inlabru.org; Bachl et al., 2019), the R package used here to implement HBM using INLA (Yuan et al., 2017b), is still in development, but aims to make these techniques more accessible to users and expand versatility. Inlabru makes it easier to use data from transects, as well as actual locations of sightings using the LGCP method. Moreover, because inlabru was motivated originally by the need to incorporate distance sampling with HBM, it is possible to implement information relevant to sighting conditions and a detection function into modelling. This was not used here, as data analysed were strip transects in which the detection function was assumed to be uniform throughout the strip width; however, this could be a benefit for line-transect surveys (e.g. boat-based surveys for marine mammals or birds), where sighting conditions (e.g. seastate, visibility, glare, etc.) influence observations strongly.
Given that models fitted with inlabru are Bayesian, the user can include prior information about parameters. Choosing appropriate priors is still an area of ongoing research and the Penalised Complexity (PC) prior approach discussed in Simpson et al. (2017) and Sørbye et al. (2019) aims at making prior choice transparent and user-driven based on the distance and strength of spatial correlation. Default uninformative priors available in inlabru (and R-INLA) were used here, but different approaches to prior choice should be assessed in the future.

Practicalities of use
The increased computing time, power, statistical complexity, and high end-user capability of fitting models using HBM with INLA, combined with the fact that they are still under development by statisticians, has meant that their uptake within the ecological community has been slow. Moreover, these newer methods have limited user-friendly documentation available (especially regarding use of the inlabru package), which may further hinder implementation (Norberg et al., 2019). Nonetheless, compared to GAMS, models incorporating INLA methodology are more flexible and can assess additional complex relationships. Moreover, HBMs are capable of modelling complex spatio-temporal data and patterns, as well as modelling interspecific interactions and species Fig. 4. Comparison of GAM and HBM results. Positive (red) indicates GAM predicted higher relative density, negative (blue) indicates GAM predicted lower density and zero is the same density. GAM vs. HBM of gridded data is on the left and GAM vs. HBM using LGCP is on the right. characteristics (e.g. joint species distribution modelling: Illian et al., 2012;Sadykova et al., 2017) which is not straight-forward in GAMs.

Management implications
In terms of translating differences between GAM vs. HBM-INLA into applied, real-world scenarios and case studies, this work has demonstrated that the loss of spatial detail in a traditional GAM, could have important implications for development of protected areas, or mitigation during fine-scale offshore operations such as pile driving, drilling, dredging, etc. For example, any elevated fine-scale clustering of porpoise at scales <1 km due to, inter alia localised bathymetric features (e.g. Smith Bank), adjacent to low-use areas, could be overlooked with a GAM, which might misinform any fine-scale designation of protected areas. Moreover, relevant to this study, for two offshore windfarms in the Moray Firth, GAMs were used to assess porpoise distribution at a minimum 4 × 4 km spatial resolution (Brookes et al., 2013;Williamson et al., 2016), which provided a quick overview of distribution or coarse estimates of relative density (e.g. for initial impact assessment or analysis at the scale of national jurisdictions); however, if fine-scale trends in animal distribution at the sub-km resolution was required to develop a plan of where to place the individual Moray Firth turbines, then LGCP analysis could have disentangled relationships between chance and feature-related clustering.
This study has demonstrated that modelling gridded data using a GAM and the negative binomial distribution highlighted two areas of high density within the Moray Firth, Scotland. Modelling the same gridded data using HBM-INLA (also with a negative binomial distribution) highlighted the same two areas, with an additional third area along the east coast that was not apparent when using GAMs. Using HBM-INLA, this time with the LGCP, (which takes spatial clustering of animals into account), showed that the Smith Bank had the area of highest relative density. Consequently, we conclude that there are three parts of the study area which host higher relative densities of harbour porpoise (one of which is missed by GAMs).

Author contributions
LW, KB, and BS conceived the study; KB collected data; PM provided satellite data; LW, JBI, and ML performed modelling; LW, BS, and VLGT led the writing of the manuscript. All authors contributed critically to drafts and gave final approval for publication.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.