Developing and assessing a density surface model in a Bayesian hierarchical framework with a focus on uncertainty: insights from simulations and an application to fin whales (Balaenoptera physalus)

Data collection

Line transect data were collected as part of the Atlantic Marine Assessment Program for Protected Species (AMAPPS) conducted by the Northeast Fisheries Science Center (NEFSC) and the Southeast Fisheries Science Center (SEFSC). The study area ranges from Halifax, Nova Scotia, Canada to the southern tip of Florida, USA and from the coastline to slightly beyond the US exclusive economic zone covering approximately 1,193,320 km² (Fig. 1). A total of 16 shipboard and aerial surveys combined were conducted from July 2010 to August 2013 covering approximately 104,000 km of line transect survey effort (Table 1). Shipboard surveys were primarily conducted during summer months in offshore waters and aerial surveys were conducted throughout the year primarily in coastal waters. Each survey included two independent observer teams. For the purpose of this study, we only analyzed data collected from the NEFSC surveys and limited the southern extent of the survey area to Cape Hatteras, NC, USA.

We divided the study site into 10 × 10 km grid cells and into 8-day temporal time periods. For each spatial–temporal cell we calculated the amount of on-effort trackline, number of sightings and obtained the corresponding values of a suite of static physiographic variables and dynamic environmental variables (Table S1). Palka et al. (2017) provides more details on the methods to collect and process the line transect and environmental data.

Model overview

A general form of a DSM for a given unit of a study area can be written as

(1)${\displaystyle E\left(n_i\right)={\hat{p}}_i{A}_{i}exp\left(B_0+\sum _j{f}_j\left(X_{ij}\right)\right),}$ where $E\left(n_i\right)$ is the expected number of sightings in unit i, ${\hat{p}}_i$ represents the probability of detection within the search area of unit i, A_i is the amount of area (e.g., km²) searched in each unit i, $E\left(n_i\right)$ can represent either the expected number of individuals sighted in unit i (Miller et al., 2013) or the expected number of groups sighted in unit i (Moore & Barlow, 2011; Pardo et al., 2015), B ₀ represents an intercept term and f_j are smooth terms for the environmental covariates X_ij.

To fit this model to line transect data using the Two-Stage Method distance sampling (Buckland et al., 2001) is used to estimate ${\hat{p}}_i$ which is then multiplied by A_i and the product is used as an offset in a GAM typically using the software mgcv (Wood, 2011). Because ${\hat{p}}_i$ only includes probability of detection at the surface, it is common practice to further adjust this estimate of detection probability with an estimate of availability at the surface if one exists (Cañadas et al., 2018). If the model is formulated such that the GAM predicts the density of groups then the density of individuals can be calculated by multiplying by an estimate of average group size where group size is simply the number of individuals in a group (Pardo et al., 2015). However, if group size varies spatially a separate spatial model for group size may need to be considered (Redfern et al., 2008).

We take a hierarchical approach to modeling the spatial density of animals where we simultaneously estimate the components of Eq. 1 using MCMC methods (Gelman et al., 2004). We use number of groups sighted in each grid cell as the response variable. To model detection probability we include a detection function model based on distance sampling and an estimate of surface availability. For the density-habitat component we use a Bayesian GAM approach to fitting smooth functions (Wood, 2016). Finally, we also include a submodel for group size to estimate average group size. Below we outline the development of each subcomponent and its implementation in a Bayesian framework. We provide a list of all parameters and definitions in Table S2.

Detection function

To estimate surface detectability we used information from the double platform survey method. Information collected from this survey design allowed us to apply mark-recapture distance sampling (MRDS) methods (Laake & Borchers, 2004). To model the sightings data from the dual observers we adopt the formulation for point independence outlined by Laake & Borchers (2004). This estimator combines a mark-recapture analysis with conventional distance sampling to estimate detection probability such that detection on the trackline (i.e., g(0)) can be estimated directly, and therefore, is not assumed to be 1. The estimator is

(2)${\displaystyle {\hat{\theta }}_i=\hat{g}\left(0,{\mathit{z}}_{\mathit{i}}\right)\frac{{\int }_0^{W}g\left(y,{\mathit{z}}_{\mathit{i}}\right)dy}{W},}$ where $\hat{g}\left(0,{\mathit{z}}_{\mathit{i}}\right)$ represents the estimate of detection probability on the trackline and is estimated from the double observer data; $g\left(y,{\mathit{z}}_{\mathit{i}}\right)$ represents the detection function at distance y and is estimated from the distance data; W is the truncation distance and z is an g × k matrix of detection covariates that influence surface detectability where g is the total number of grid cells and k is the total number of detection covariates such as Beaufort sea state, glare, etc. To model $g\left(y,{\mathit{z}}_{\mathit{i}}\right)$, we considered half-normal and hazard rate detection functions. To model $g\left(0,{\mathit{z}}_{\mathit{i}}\right)$, we adopted the approach outlined by Laake & Borchers (2004). Specifically, we modeled the binary outcome of whether or not an observer successfully detected an animal group that was present at distance y as the outcome of a Bernoulli trial. Further details of analyzing the double platform line transect data are provided in Appendix S1.

For the aerial surveys, the secondary team was positioned toward the back of the plane but had an obstructed view of the trackline complicating a direct implementation of the MRDS approach. Therefore, we estimated an average g(0) for the aerial surveys independently where we treated the front team as a single platform and estimated g(0) using a trial configuration (i.e., using detections by the rear observers as “trials” for the front observers). The resulting estimate and coefficients of variation (CVs) for g(0) was 0.67 (CV = 0.16). We used this estimate to develop an informative prior in the Bayesian Method. Information on estimating g(0) and applying it to the aerial data is provided in Appendix S2.

For each survey, we determined the best detection function through a stand-alone MRDS analysis using the mrds package in R. A hazard rate likelihood provided the best fit to both the aerial survey data and the shipboard data. Because sample sizes were low for sightings that were positively identified as fin whales, we pooled these sighting with ambiguous sightings that were either a sei whale (Balaenoptera borealis) or a fin whale (see Table S3). We compared models using AIC and the top model structure for each survey was used in the Bayesian Method (see Table S3).

Surface availability

Because most marine animals spend some amount of time below the surface there is a need to also correct for surface availability (a) (Laake et al., 1997; Forcada et al., 2004). We developed an informative prior distribution for this parameter using the estimate and corresponding standard error of surface availability for fin whales reported in Palka et al. (2017). Their method for estimating a was based on the probability of an animal being at the surface and available for detection during a survey, and took into consideration the species diving and aggregation behaviors, in addition to the amount of time the observer had to analyze any spot of water from each of the survey platforms. This correction tended to be larger for aerial surveys than for shipboard surveys as the window of observation is considerably smaller for aerial surveys. The estimate for fin whales for aerial surveys was 0.37 (CV = 0.34) and assumed to be 1 for shipboard surveys (see Appendix S2 for more details).

Our final survey-specific correction for detection probability in each grid cell i is:

(3)${\displaystyle p_i^S={\Theta }_i^S}$ for the shipboard surveys (S) and

(4)${\displaystyle p_i^A={\Theta }_i^{A}a}$ for the aerial surveys (A). To investigate possible correlations between the density-habitat model and detection probability model, we also performed a stand-alone analysis of the distance data with the Bayesian Method by omitting the density-habitat model in a manner analogous to the stand-alone MRDS analysis.

Density-habitat function

We take a GAM approach to parameterize the density-habitat function. The general basis expansion for a single smooth term can be written as

(5)${\displaystyle f\left(x\right)=\sum _{k=1}^K{\beta }_k{b}_k\left(x\right),}$ where b_k(x) are basis functions and β_k are parameters to be estimated. The basis size K is usually chosen by the user to be large enough to allow an appropriate amount of flexibility in $f\left(x\right)$. To avoid overfitting quadratic penalty terms are included which take the form

${\displaystyle \sum _k{\gamma }_k{\beta }^T{S}_k\beta ,}$ where S_k are penalty matrices and γ_k are smoothing parameters to be estimated.

A multivariate normal distribution can be constructed for the GAM parameters as ${\displaystyle \beta \sim MVN\left(0,\sum _k{\gamma }_k{S}_k\right)}$ where ∑_kγ_kS_k represents a precision matrix (instead of the more common variance–covariance matrix) and the penalty terms are given a vague, gamma prior such as

${\displaystyle {\gamma }_k\sim Gamma\left(0.05,0.005\right),}$

The terms can be estimated efficiently using Gibbs sampling with conjugate priors.

To calculate the precision matrices we used the jagam function in the R package mgcv (Wood, 2016). This function allows the user to specify a number of different smooths (cubic splines, tensor products, etc.) and provides the basic code and input of a JAGS model. In addition, it centers the smooths to facilitate faster convergence.

Likelihood

We developed a likelihood by assuming count data followed a Tweedie distribution which has been shown to provide a good fit to cetacean data (Miller et al., 2013; Roberts et al., 2016). The Tweedie distribution is a three parameter family of distributions that can take the form of more commonly used distributions such as the normal, Poisson and gamma. If the power parameter (ρ) is in the range 1<ρ<2 than the distribution can also be referred to as the compound Poisson-gamma (CPG). Because the Tweedie random variables are a sum of G gamma variables (M) where G is Poisson distributed (Jorgensen, 1987), it can be expressed in terms of a Poisson and a gamma distribution such that ${\displaystyle G_i\sim Poisson\left({\lambda }_{\mathrm{p}}\right)}$ ${\displaystyle M_{ij}\sim Gamma\left(\alpha ,\beta \right)}$ where

(6)${\displaystyle {\lambda }_g=\frac{\alpha }{\beta }}$ and

(7)${\displaystyle Y_i=\left\{\begin{array}{cc}\sum _{i=1}^{G_i}{M}_{ij}\hfill & G_i>0\hfill \\ 0\hfill & G_i=0\hfill \end{array}\right.}$ where Y_i represents the response variable and the expectation of Y_i is then E(Y) = λ_p λ_g. Lauderdale (2012) shows that under a specific parameterization, the coefficients of the regression model can be estimated by estimating both the Poisson and gamma components separately. Specifically, this parameterization can be written as

(8)${\displaystyle {\lambda }_P=e^{\frac{\mathit{X}\left(\beta -\varphi \right)}{2}}}$

(9)${\displaystyle {\lambda }_g=e^{\frac{\mathit{X}\left(\beta +\varphi \right)}{2}}}$

where X is a design matrix, β is a vector of regression coefficients and ϕ is a vector of coefficients that control the extent to which the regression coefficients vary between the Poisson component and gamma component of the compound distribution (Lauderdale, 2012).

Group size

To model group size we use a zero truncated Poisson such that the group size of each sighting is modeled as ${\displaystyle \left(s_t-1\right)\sim Poisson\left({\lambda }_s\right),}$ where s_t is the t_h observation of group size and λ_s+ 1 represents the average group size. This approach assumes that group size is unrelated to detection probability. This assumption is supported by our analysis of fin whale sightings data which did not indicate a strong influence of group size on detection probability. Moore & Barlow (2011) used a similar approach to model group size in their analysis of fin whale data where they also concluded that group size was unrelated to detection probability.

Density and abundance estimation

The density of individuals within a grid cell is the product of group size, the density of groups in grid cell i and the area of grid cell i such that

(10)${\displaystyle N_i=\left({\lambda }_s+1\right)D_i{A}_i^G}$

where D_i represents the density of groups in grid cell i from the GAM, λ_s + 1 represents the average group size (s), $A_i^G$ is the total area of grid cell i and N_i is the predicted number of individuals in grid cell i. The total abundance of individuals within the study area is the sum N_i over all grid cells within the study area.

Model fitting

We fitted the Bayesian Method outlined above using MCMC sampling implemented with the JAGS software (Plummer, 2003). We used vague prior distributions for all parameters with the exception of g(0) for the aerial surveys and $\hat{a}$ where we used estimates and associated CVs to develop informative beta prior distributions. We included a burnin of 20,000 samples and two chains of 50,000 with a thinning rate of 50. We assessed convergence by examining trace plots and calculating Gelman–Rubin diagnostics (see Table S4 for results).

Simulation study

To evaluate the Bayesian Method and compare its performance to the Two-Stage Method we used a simulation study. We simulated spatial variation in abundance among 1,000 hypothetical grid cells. We next simulated line transect sampling in a subset of the grid cells. We fitted both the Bayesian Method and the Two-Stage Method to each of 500 simulated datasets to estimate abundance and compare results to the true abundance used in the simulation. We use these results to quantify bias and precision of each method and evaluate statistical interval coverage. We provide a more detailed summary of the simulation study in Appendix S3.

Application to fin whales

To further evaluate our modeling approach and investigate sources of uncertainty we analyzed a four year dataset of fin whales sightings collected during the AMAPPS surveys with both the Bayesian Method and the Two-Stage Method. For the Two-Stage Method we used the formulation of Miller et al. (2013) where we used estimates from the stand-alone MRDS analysis as an offset for $\widehat{p_i}$ and A_i is the same as described above for the Bayesian Method. For the density-habitat model, we fitted GAMs in R using the mgcv package version 1.8-28 (Wood, 2014). We used thin plate regression splines and restricted maximum likelihood (REML) to estimate parameters. We set the basis size (K) to 5 for all smoothing functions which was large enough to allow enough flexibility in fitting the smooth terms. To account for overdispersion, a Tweedie distribution was assumed. We used a combination of Akaike’s Information Criterion (AIC) and deviance explained to determine the best set of covariates and the best structure of the smooth terms. Because our interest was in keeping as many components of the model structure consistent across methods for the purpose of comparison, we did not perform model selection with the Bayesian Method and instead opted to use the same model formulation for both the detection function and the density-habitat model in the Bayesian Method that we used in the Two-Stage Method.

To compare the Bayesian Method to a scenario where there is no variance propagation, we intentionally did not propagate uncertainty from the first step when applying the Two-Stage Method such that the precision of final estimates only reflected the density-habitat model uncertainty. To further investigate how much each component of a DSM contributes to uncertainty of the final estimate of abundance we ran several parallel analyses with the Bayesian Method where in each analysis we fixed one component from the first step to its point estimate such that the uncertainty in that component was not propagated. For the detection function components, we used point estimates from the stand alone MRDS analysis. Analogous to the Two-Stage Method, we also ran a model in the Bayesian framework where we fixed all components to their point estimates such that only the density-habitat model uncertainty was represented. For each analysis, we calculated an abundance estimate and CV and compared changes in CV to the full Bayesian model.

Simulation study

Results from the simulation study showed that the Bayesian Method was able to achieve close to the nominal 95% rate of coverage with low bias. In addition, results were comparable to the Two-Stage Method. Estimated coverage for both methods was 94.9%. The average CV of abundance estimates was 0.16 for the Bayesian Method which was similar to the average CV of 0.14 for the Two-Stage Method although the distribution of CVs were more positively skewed for the Bayesian Method (Fig. S1). Bias was low (<1.5%) for both methods although marginally more negative in the Bayesian Method (Fig. S1).

Application to fin whales

A comparison of the resulting detection functions between the stand-alone MRDS and the Bayesian Method showed overall detection probabilities were similar although higher in the Bayesian Method (Table 2). For the shipboard surveys, estimates of detectability from the distance sampling component was approximately 13% higher than the estimates from the stand-alone MRDS whereas estimates of g(0) were the same. For the aerial surveys estimate of detectability from the distance sampling component was approximately 17% higher than the stand-alone MRDS analysis. Estimates of detection probabilities from the stand-alone Bayesian analysis were closer to the stand-alone MRDS analysis then estimates from the full Bayesian Method for the aerial survey whereas estimates for the shipboard survey were unchanged (Table 2).

The posterior estimate of mean group size was 1.4 (CV = 0.14). Most observed group sizes were fewer than 2 animals with approximately 4% greater than 3 animals (Fig. S2).

The top model for the density-habitat function included distance to shore, depth, sea surface temperature and distance to 125 m isobath as covariates. Results from fitting the Bayesian Method to the observed sightings data showed good agreement between predicted and observed number of groups per grid cell although there was some tendency of the model to under predict as the number of sightings increased (Fig. S3). In comparison to the Two-Stage Method, density estimates for the grid cells during summer were similar between the two methods (Fig. 2). This similarity in model performance resulted in predicted density distribution patterns that were largely similar among methods (Figs. 3A and 3B). The pattern in uncertainty was also similar for the two methods while the Bayesian Method produced noticeably higher CVs overall (Figs. 3C and 3D).

The analysis of the variance components suggested that variance associated with detectability had a significant influence on the precision of the final abundance estimate whereas variance in mean group size made a negligible contribution (Table 3). The CV for the abundance estimate decreased by 46.9% and 40.6% when we fixed detection probabilities for the aerial survey and shipboard survey to their point estimates, respectively (Table 3). The effect of fixing only surface availability to its point estimates also reduced the CV substantially (Table 3). The CV for the abundance estimate from the full Bayesian Method was greater than 2x the CV from the Two-Stage Method. The abundance estimate and CV for the Bayesian Method were similar to the Two-Stage Method when all components were fixed in the Bayesian Method (Table 3). Abundance estimates of fin whales in the study area ranged from 4012 to 4551 depending on which components were fixed to their point estimates and which method was used (Table 3).