Spatially explicit models of density improve estimates of Eastern Bering Sea beluga (Delphinapterus leucas) abundance and distribution from line-transect surveys

Design-based estimator

Although our primary focus is on developing density surface models for EBS belugas, previous abundance estimates for this stock were generated using a design-based estimator. A basic Horvitz-Thompson-like line-transect estimator of animal density is (Buckland et al., 2001; Burt et al., 2014): (1)${\displaystyle \hat{D}=\frac{1}{a}\sum _{j=1}^{n_g}\frac{S_j}{\hat{p}\left({\mathbf{z}}_j;\hat{\theta }\right)}}$ where

n_g

: number of groups detected;

S_j

: size of group indexed by j;

a

: area searched during line-transect survey, where a = 2Lw, L is the total length of transects surveyed, and w is the width of the strip searched on one side of the aircraft;

$\hat{p}\left({\mathbf{z}}_j;\hat{\theta }\right)$

: model-based estimate of the overall probability that an observer detects group j, given covariates z_j that affect detectability. This term accounts for all sources of perception and availability bias (Marsh & Sinclair, 1989; S4);

$\hat{\theta }$

: parameter estimates required to estimate detection probabilities.

To derive an estimate of the total number of animals in the study area ($\hat{N}$), we multiply the total study area size, A, by the density estimate from Eq. (1): (2)${\displaystyle \hat{N}=A\hat{D}.}$

Under this formulation, inference proceeds by first fitting detection function models to observed distances and other covariates to produce estimates of detection parameters (i.e., $\hat{\theta }$), before applying Eq. (1) in a second step. The abundance estimator in Eq. (2) is unbiased if certain assumptions about the survey design and realization hold (Buckland et al., 2001; Hedley & Bravington, 2014). Hence, this $\hat{N}$ is referred to as a design-based estimator.

In addition to allowing calculation of the design-based estimator, we also use these estimates of detection probability when fitting DSMs. In the following, we shorten notation such that $p_j=\hat{p}\left({\mathbf{z}}_j;\hat{\theta }\right)$ for sightings and $p_i=\hat{p}\left({\mathbf{z}}_i;\hat{\theta }\right)$ for segments. For more information about detection probability calculations, see S4.

DSMs: marginal likelihood

As with most DSM implementations, we construct a spatial model for counts of individuals, which in our case were summarized over 10-km transect segments (see Eastern Bering Sea beluga case study for further information). For all DSMs, we write a generic marginal likelihood of a parameter vector, ξ, given observed counts of individual animals, c, and other known variables, x, as (3)${\displaystyle \mathcal{L}\left(\xi ;\mathbf{c},\mathbf{x}\right)={\int }_{\eta }\left[\mathbf{c}|\xi ,\eta ,\mathbf{x}\right]\left[\eta |\mathbf{x},\xi \right]d\eta .}$

Here, [c|ξ, η, x] is the conditional probability density function of observed counts, given parameters, random effects (η), and known covariates. The counts represent the number of animals detected on rectangular transect segments (with width 2w, as in the design-based estimator). The component [η|x, ξ] represents the distribution of random effects. We use the integral to indicate that the random effects will be integrated out of the joint likelihood—in our case using the Laplace approximation available in TMB software (Kristensen et al., 2015). As usual in likelihood-based inferential statistics, the likelihood is viewed as a function of the unknown parameters, ξ.

In order to derive [c|ξ, η, x], we must first specify a suitable probability mass or density function. Although it is customary to specify probability mass functions for count data, initial exploration of Poisson and negative binomial distributions indicated considerable lack-of-fit when applied to our EBS beluga data set. Specifically, model diagnostic plots examining the relationship between the mean and variance in the residuals compared to the theoretical distribution (Ver Hoef & Boveng, 2007), and quantile residuals computed using a probability integral transform (PIT; Dunn & Smyth, 1996) and visualized using the R package DHARMa (Hartig, 2022), showed that Tweedie distributions (Jørgensen, 1987; Dunn & Smyth, 2005; Kendal, 2004) provided a better fit to the data. Therefore, we adopted a parameterization based on the Tweedie distribution.

The Tweedie distribution provides increased flexibility compared to the Poisson and negative binomial distributions, allowing a diversity of shapes and accommodating zero-inflation. It is a specific case of an exponential dispersion model, with mean µ, and variance V(μ) = ϕμ^ρ (Dunn & Smyth, 2005). We specifically set the range of ρ to be on 1 < ρ < 2, a parameterization variously known as “compound Poisson”, “compound gamma”, or “Poisson-gamma” (Dunn & Smyth, 2005; Kendal, 2004). This distribution has support on the non-negative real line, although authors often use this distribution for non-negative integers (e.g., counts; Kendal, 2002; Miller et al., 2013; Sigourney et al., 2020), which is our approach in this paper. Kendal (2002) and Kendal (2004) discuss the relationship between the Tweedie distribution and Taylor’s power law in ecology, which explains clustered spatial distributions as manifestations of power function relationships between the variance and mean number of organisms in an area (Taylor, 1961).

For 1 < ρ < 2, the Tweedie distribution does not have a closed form, but can be evaluated numerically (e.g., using the ‘dtweedie’ function in the TMB library). Therefore, we symbolically write (4)${\displaystyle \left[\mathbf{c}|\xi ,\eta ,\mathbf{x}\right]\equiv \prod _i\text{Tweedie}\left(c_i;{\mu }_i,\varphi ,\rho \right).}$ That is, the joint likelihood of observed counts on transect segments indexed by i is a product of conditionally independent univariate Tweedie density functions, with mean μ_i and constant dispersion and power parameters, ϕ and ρ. The mean, μ_i, is a function of fixed and random effects, and “known” detection probability and survey coverage offsets, such that (5)${\displaystyle \mu =exp\left({\beta }_0+\delta +log\left(\mathbf{a}\right)+log\left(\mathbf{p}\right)\right).}$ Here, β₀ represents an intercept parameter (no other fixed effects were included in our models), δ is a vector of ‘realized’ random effects for transect counts, a is a vector of the area surveyed for each transect segment (a_i = 2L_iw), and p is a vector of the overall detection probability (including both availability and perception bias corrections) for each segment (p_i, see S4 Eq. 6).

Note that this parameterization requires that any covariates used to estimate detectability relate only to the transect segment; observation-specific covariates (e.g., color, group size) cannot be used in this parameterization. See Miller et al. (2013) for an alternative parameterization that allows observation-specific covariates by specifying the response variable to be an estimate of bias-corrected abundance.

The actual dimension of random effects often differs from the number of transect segments. Specifically, we model δ = Aη, where the matrix A has dimension (n_i, n_η), with n_η denoting the true number of random effects, and n_i the number of transect segments. Next, we elaborate on the random effects specifications.

Random effects specifications

We have yet to describe [η|x, ξ] in Eq. (3). This component defines the specification of spatially autocorrelated random effects. For a given survey period (year) the data likelihood Eq. (4) is the same for all the models that we considered, so the random effects specification is the only difference among the DSMs we developed.

For all models, random effects were assumed to be drawn from a multivariate normal distribuion with mean zero, and a spatially patterned covariance matrix, Σ: ${\displaystyle \eta \sim \text{Multivariate normal}\left(\mathbf{0},\Sigma \right).}$ Spatial autocorrelation is imparted by constraints on the (n_η × n_η) Σ matrix. In practice, we chose to work with a precision matrix Q = Σ⁻¹, which was often sparse, enabling greater computational efficiency.

We employed two related, but conceptually different types of models to specify Q: stochastic partial differential equations (SPDEs) to approximate Matérn geostatistical models (Lindgren, Rue & Lindström, 2011), and spline-based models. The latter are commonly used in generalized additive models (e.g., Wood, 2006), where smooth terms are often viewed as penalized fixed effects. However, it is also possible to conceptualize smooths as mean-zero random effects, with an associated precision matrix (Miller, Glennie & Seaton, 2020), which is the approach that we used in this paper. Further details on how Q is specified for individual models is provided below.

We conducted all analyses in the R programming environment (version 4.3.2; R Core Team, 2023) using the TMB package (Kristensen et al., 2015) to formulate marginal log-likelihoods and generate DSM parameter estimates, the mgcv package (Wood, 2006) to set up spatial spline bases, and R-INLA (Rue, Martino & Chopin, 2009) to create Delaunay triangulation meshes for Stochastic Partial Differential Equation (SPDE) models. All data and code used in this paper have been uploaded to a public GitHub repository at https://github.com/meganferg/FergusonEtal_20250125_EBS_Beluga_DSM and are available on Zenodo at https://doi.org/10.5281/zenodo.16269135.

Prediction

For each model and year of analysis, we summed predictions of the number of belugas in each hexagonal grid cell h in our study area (see Eastern Bering Sea beluga case study for how these were defined), applying epsilon bias-correction to this total (see Correcting for detransformation bias). Specifically, we calculated (8)${\displaystyle {\hat{N}}_h=exp\left({\hat{\beta }}_0+{\hat{\delta }}_h+log\left(a_h\right)\right).}$ Note that a_h gives the area of ocean in hexagon h (i.e., omitting land). There is no offset for detection probability because we are interested in all belugas, not just those that are detectable and detected. The vector of “realized” random effects $\hat{\delta }$ is calculated as $\hat{\delta }={\mathbf{A}}^{pred}\hat{\eta }$ where $\hat{\eta }$ is the value of η that maximizes the joint likelihood conditional upon the maximum likelihood estimator (MLE) $\hat{\xi }$ for fixed effects (termed the empirical Bayes estimator for η). This predictor for ${\hat{N}}_h$ is called the “plug-in estimator” because it plugs in the empirical Bayes estimator as if it were fixed. For SPDE models, the (n_h, n_η) interpolation matrix, A^pred was constructed using the “inla.spde.make.A” function in R-INLA, using the centroid of each hexagon as the prediction location. For spline-based smoothers, we obtained A^pred using mgcv’s “predict()” function with “type=lpmatrix”, again using the centroids of each hexagon as prediction locations. An estimate of total abundance is then calculated as $\hat{N}={\sum }_h{\hat{N}}_h$. The epsilon bias-correction was applied to $\hat{N}$.

Model evaluation and final candidate model selection

To evaluate DSMs and select the final candidate DSMs for the ensemble model, we advocate using several criteria, including: examining conditional PIT residuals via the DHARMa package (Hartig, 2022); extrapolation metrics (defined below); and visual examination of maps showing the DSM predictions overlaid with the sightings used to build the models. We provide a detailed example in the Eastern Bering Sea beluga case study.

To identify models whose predictions might be unreliable due to extrapolation bias, we considered two types of ad hoc metrics. First, for each combination of model and cell (i.e., hexagon, h), we computed the following ratio: (9)${\displaystyle \frac{{\lambda }_{m,h}}{{\lambda }_{m,max}}.}$ Here λ_m,h is the predicted abundance from model m for unsampled location h, and λ_m,max is the maximum predicted abundance across all sampled cells. Second, for each model we counted the number of unsampled cells (i.e., cells that did not have line-transect survey effort) with predicted abundance exceeding the maximum predicted abundance in sampled cells. These procedures are motivated by a generalized version of Cook’s independent variable hull (Cook, 1979; Conn, Johnson & Boveng, 2015).

Uncertainty estimation

Ensemble model

Fitting multiple DSMs to sightings raises the question of which model, or collection of models, should be used to generate a final abundance estimate and density surface. The question is particularly important when different models produce markedly different estimates of abundance. We chose to base ultimate inference on an ensemble (Araújo & New, 2007), whereby estimates from different models are averaged to produce a final estimate. Specifically, we compute (13)${\displaystyle {\hat{N}}_{ens}=\sum _m{w}_m{\hat{N}}_m}$ where ${\hat{N}}_m$ is the MLE of abundance from TMB for each model m. The model weight is w_m, where ∑_mw_m = 1.0. The advantage of averaging models is that there is often a reduction in prediction error (Burnham & Anderson, 2002; Dormann et al., 2018).

There are different approaches for setting the model weights (Dormann et al., 2018). For instance, a common approach is to use Akaike’s information criterion (AIC) associated with fitted models to calculate weights (Burnham & Anderson, 2002). However, calculation of AIC weights relies on the complexity of a model, often computed as the effective degrees of freedom from a generalization of the hat-matrix, and this is difficult to compute in a hierarchical model using maximum-likelihood methods. Instead, we used equal model weights, which have been shown to perform well in prediction of species distributions (Dormann et al., 2018). This procedure has the added advantage that a single model with an extremely high or low abundance estimate will not dominate inference.

The variance of model-averaged predictions was calculated using the standard unconditional variance estimator (i.e., Burnham & Anderson, 2002): (14)${\displaystyle \widehat{Var}\left({\hat{N}}_{ens}\right)={\left[\sum _m{w}_m\sqrt{Var\left({\hat{N}}_m\right)+{\left({\hat{N}}_m-{\hat{N}}_{ens}\right)}^2}\right]}^2.}$

Eastern Bering Sea beluga case study

The EBS beluga survey design and data collection methods are similar to those previously described in Ferguson, Conn & Thorson (2024); the analytical methods and results reflect improvements to Ferguson, Conn & Thorson (2024). The data collection methods and sighting and effort summaries are presented in S2 for the aerial line-transect surveys and S3 for the aerial imagery. All aerial surveys were approved by the Alaska Fisheries Science Center/Northwest Fisheries Science Center Institutional Animal Care and Use Committee (IACUC Number NW/AK2016-6). The analytical methods used to estimate detection probabilities and corresponding results are presented in S4. See Frost, Lowry & Nelson (1985) and Frost & Lowry (1995) for details about the VHF telemetry data and analyses. There were no estimates of uncertainty for availability probability (Ferguson et al., 2023); therefore, this parameter was included as a known constant in the offset for the DSM Eq. (5).

To derive spatially-explicit estimates of EBS beluga abundance, we constructed density surface models separately for each year, 2017 and 2022. The 2017 study area boundary corresponded to the area covered by the geographic strata in Lowry et al. (2017) (Fig. 1). The ABWC advocated for a larger survey area in 2022. They noted that Indigenous knowledge has confirmed that the southern extent of the EBS beluga stock’s distribution during early summer extends farther south than the historical survey boundaries. Therefore, the 2022 survey area extended south of  Hooper Bay (Fig. 1). The study area boundary used to estimate abundance in 2022 excluded Scammon Bay because that bay was not included in the 2022 survey design, there was no survey effort inside the bay (i.e., east of the barrier islands), and we have insufficient information about beluga habitat to make inference about beluga density in Scammon Bay based on density in the area surveyed outside Scammon Bay.

DSMs were constructed using aerial line-transect sighting and effort summaries for 10-km segments of transect effort. This segment length is approximately the distance between adjacent transects (9.3 km). The segments were created by sequentially slicing transect effort conducted in Beaufort Sea State ≤4, beginning with the start of each transect. End segments <10 km were added to adjacent segments so that all segments used in the analysis were ≥10 km. Predictions from the DSM were based on a hexagonal grid with cell midpoints located 10 km apart. All geospatial data were projected into an equidistant conic projection (false easting: 0.0; false northing: 0.0; central meridian: −164.0°; latitude of origin: 63.5°; standard parallels: 62.5°, 64.5°; WGS84 datum; linear unit: kilometer), and this projection was used when calculating cell areas and distances for the spatial correlation functions or splines.

The DSMs required segment-specific estimates of detection probability, p_i Eq. (5). The best-fitting MCDS detection function for EBS belugas included covariates for Beaufort Sea State (integer-valued) and turbidity (binary) (S4). To build the DSMs, effort data for these variables were summarized by segment. The segment-specific Beaufort Sea State variable was calculated as the average value of integer-valued Beaufort Sea State for all records that were located on the segment; all records were weighted equally. The segment-specific turbidity variable was calculated by assigning the binary turbidity variable an integer value (no = 0; yes = 1), computing the average of the integer-valued turbidity values for all records located on the segment, and rounding the result. For example, if segment i comprised three data records with turbidity “yes”, “yes”, and “no”, the average of their integer-valued analogs would be 1 + 1 + 0 = 0.67, which rounds to 1, so the segment would be designated as turbidity = “yes”.

For 2017 and 2022, 13 and 17 DSMs, respectively, were constructed and examined. For each year, we allowed at most one Matérn and one SPDE Matérn model with barriers to be included in the ensemble. Overviews of key aspects of each model are provided in Tables S6.1 and S6.2. To ensure that all models for a given year were allowed the same flexibility to construct the spatial surface, we specified the DSMs for each year using approximately the same number of random effects. For each year, the number of random effects for all DSMs was based on the number of random effects (i.e., vertices) in that year’s SPDE Matérn model, which was determined using the guidelines provided in Belmont (2022). The knot locations for the soap film smoothers corresponded to the locations of the vertices for the corresponding SPDE Matérn model that were within the soap boundary. This resulted in fewer total knots for the soap film smoother than the corresponding SPDE Matérn model; however, both models had the same number of knots inside the soap boundary, and the soap film smoother was specified to have the same number of random effects as the corresponding SPDE Matérn model. The total number of random effects for the bivariate tensor product smoother corresponds to the product of each of the basis dimensions used to build the smoother, and it was not possible to exactly match the number of random effects in the corresponding SPDE Matérn model.

Identical DSMs were fit in mgcv and TMB, with the exception of the barrier SPDE models, for which mgcv functions defining this type of basis were not available, so they were constructed only in TMB. We constructed identical models in both software platforms for two reasons: (1) to apply the methods presented in Miller, Glennie & Seaton (2020) to the EBS beluga data and confirm that nearly identical results could be derived using mgcv and TMB; and (2) to evaluate whether the existing “ds_varprop” function in the dsm package (Bravington, Miller & Hedley, 2021; Miller et al., 2022) would be an alternative to Eq. (10) and the methods presented in S5 for propagating uncertainty from the MCDS detection function model into the overall estimate of uncertainty in the abundance estimate. However, as of dsm version 2.3.3, the “dsm_varprop” function could not propagate errors through SPDE models (M Ferguson with D Miller, pers. comm., 2023).

Because the variance in the MLE of abundance from TMB for each candidate DSM, $Var\left({\hat{N}}_m\right)$, did not include uncertainty from the estimate of transect detection probability in the beluga case study (i.e., ${\hat{p}}_{MR}\left(0,{\mathbf{z}}_j;{\hat{\theta }}_{MR}\right)$ in S4), we applied the delta method to add this component of uncertainty to $\widehat{Var}\left({\hat{N}}_{ens}\right)$, resulting in ${\widehat{Var}}_{tot}\left({\hat{N}}_{ens}\right)$ (see Eq. 4 of S5).

For comparison with the 2017 EBS beluga DSM results, we examined the conventional abundance estimate derived using the post-stratified design-based abundance estimator of Ferguson et al. (2023). Their variance estimator had three components: (1) variation from uncertainty in estimating the MCDS parameters; (2) variation from uncertainty in estimating transect detection probability; and (3) variation in abundance due to random sample selection. Two conventional abundance estimates for 2022 were derived using the post-stratified design-based estimator. One conventional estimate for 2022 used the same strata as the conventional abundance estimate for 2017 (i.e., the Lowry et al. (2017) strata; Ferguson et al. (2023)). The other conventional estimate for 2022 used geographic strata covering the full extent of the 2022 study area (Fig. 1). The modified geographic strata for 2022 incorporated the latitudinal boundaries from Lowry et al. (2017) (Figs. 1, S2.2). Unlike the Lowry et al. (2017) strata, the modified 2022 strata extended to the western border of the 2022 study area and two new strata were created to encompass the northernmost and southernmost portions of the 2022 study area. Both conventional abundance estimates for 2022 used the methods of Ferguson et al. (2023), with a MCDS detection function model based on the pooled data from 2017 and 2022 (discussed briefly above and in S4).

The detection function model used in the post-stratified design-based estimator for 2017 was based on only a single year of data (Ferguson et al., 2023), whereas the detection function model used in our model-based estimators for 2017 and 2022 was based on data from both years pooled (S4). However, the CV of the former detection function model was 0.043 and the CV of the latter was 0.037, only trivially smaller; therefore, we do not believe that this difference in detection function models affected our overall comparison of the precision in the different abundance estimators.