Extending density surface models to include multiple and double-observer survey data

Differing density by platform/observation type

Including multiple platforms as in (2) assumes that any differences in observed counts are a result of detectability alone and that mean density does not differ between the two platforms. This may be unrealistic, especially when surveys were conducted from different platforms at different times (e.g., seasonal surveys) or when platforms operate simultaneously but record animals in differ behavioural states (e.g., SAS). A simple adaptation to model (2) would be to include a per-platform intercept, β_k:

(3) $$\mathbb{E}\left[n_{j_k}|\beta ,\lambda ,p_k({\hat{\theta }}_k;{\mathbf{z}}_{j_k})\right]=a_{j_k}{p}_k({\hat{\theta }}_k;{\mathbf{z}}_{j_k})\exp \left({\beta }_k+\sum _m{f}_m(x_{m,j_k})\right).$$

Model (3) assumes that the density only shifts the intercept via β_k and has no effect on the f_ms, which may be overly restrictive. We can extend our model using the hierarchical GAM framework (Pedersen et al., 2019) to allow the smooth parts of the model to vary by platform (factor-smooth interactions). We allow information to be shared between these smooths, such as how wiggly they are (shared smoothing parameters) or that they have similar shapes (smoothing towards a global term). To do so, we extend (2) as:

(4) $$\mathbb{E}\left[n_{j_k}|\beta ,\lambda ,p_k({\hat{\theta }}_k;{\mathbf{z}}_{j_k})\right]=a_{j_k}{p}_k({\hat{\theta }}_k;{\mathbf{z}}_{j_k})\exp \left({\beta }_0+\sum _{m_1}{f}_{m_1}(x_{m_1,j_k})+\sum _{m_2}{f}_{m_2}(x_{m_2,j_k},k)\right),$$where the m₁ smooths are as in (3) and m₂ are factor-smooth interactions. Using a spatial smooth (f(x, y)) as an example, we can then fit a spatial smooth for each platform: we have f(x, y, k) for k = 1,…, K. We may choose a subset of terms in the model that seem most likely to vary by platform or, include all terms as factor-smooths. Pedersen et al. (2019) enumerate all the models possible under this framework.

Incomplete detection at zero distance and availability

A fundamental assumption of distance sampling is that objects at zero distances are observed with certainty: that is, g(0) = 1 if g is the detection function (Buckland et al., 2001). In one approach to dealing with a violation of this assumption, two observers (or teams of observers) can be used to set-up a capture-recapture experiment where the probability of observing an animal at zero distance is estimated by considering one observer as setting-up trials for the other when animals are detected (mark-recapture distance sampling; MRDS). Using this approach we estimate g(0; z_g(0), θ_g(0)) where θ_g(0) are parameters specifically for the estimation of g(0) and z_g(0) are (optional) covariates. Burt et al. (2014) give details of models for g(0; θ_g(0)). Here we consider only the independent observer mode of MRDS where each observer’s detections are trials for the other.

Including MRDS models into the DSM framework is simply a case of using g(0; θ_g(0)) as an additional multiplier on the detection probability in the offset. We re-write $p(\hat{\theta };{\mathbf{z}}_i)$ in model (1) as the product $g(0;{\theta }_{g(0)})p(\hat{\theta };{\mathbf{z}}_i)$. Accordingly we can index the MRDS model with k and include it as one of our multiple surveys in models (2), (3) or (4).

As we will see below, we can also account for animals not being available to be detected (e.g., due to diving behaviour) in our models by using an estimate from other sources (such as tag data) and including this in the same way as g(0). For now we assume that the estimate of availability is independent of the surveys (i.e., data were collected at a different time, using different methods, perhaps in a different area) and therefore that we can add the squared coefficient of variation to the DSM’s to obtain a total uncertainty on our abundance estimates. Availability estimates can be specified at the segment level, so can vary between platform (see below).

Strip and plot sampling

If we assume that all objects within a given distance of the sampler are detected with certainty, we have strip or plot transects (as analogues to line and point transects). Strip transects are common when video/photo surveys are conducted from planes or drones, so detectability is not an issue or, when the truncation distance is sufficiently small that it is believed that all objects will be detected. In this case we simply replace the $p(\theta ;{\mathbf{z}}_i)$ in model (1) with one.

Variance estimation

Each fitted detection function included in the DSM has its own covariance matrix for that model’s parameters. If we assume independence between the detection functions, we can create a joint covariance matrix as a block diagonal matrix with each block representing one of the K platforms: V_θ = diag(V_θ1, V_θ2, …, V_θK), where V_θk is the covariance matrix for the parameters of the detection function for platform k. When a strip/plot transect is used there is no corresponding θ_k and therefore no uncertainty, so it does not appear in V_θ.

If we are willing to assume independence between the detection function and spatial model, we could apply the delta method (Seber, 1987) to combine the detection function and GAM variances. The delta method is unappealing if there are implicitly spatial covariates in the detection model (e.g., wind speed varies in space) as there will be non-zero covariance between those model components. Instead we can adapt the variance propagation approach of Bravington, Miller & Hedley (2021) to include multiple detection functions. This approach uses a quadratic approximation to the detection function in a refit of the spatial model to adjust the detectability estimates in light of the additional spatial information from the GAM covariates. The quadratic adjustment takes the form of a random effect with mean zero and covariance matrix V_θ, so the detection function uncertainty is included in the final variance estimate. In this way V_θ can be plugged-in to the variance propagation method of Bravington, Miller & Hedley, 2021 and a posterior covariance matrix for all model components (V_{β, θ}) can then be estimated. Note that the assumption of independence is made a priori but the variance propagation procedure can estimate the off-diagonal elements of V_{β, θ}. In the next section we show how this uncertainty about the model parameters can be applied to the uncertainty in abundance estimates.

Estimating abundance and its variance

DSMs are most often used to provide estimates of abundance over the study area (or some subset(s) of it), which are calculated as sums of predictions over grids (Miller et al., 2013). We can express this in matrix form as $\hat{N}=\mathbf{a}\exp {\mathbf{X}}_p{\beta }_{}$, where a is a row vector of grid cell areas, X_p is the matrix that maps model coefficients to predictions (a design matrix for the predictions) and $\hat{\beta }$ are the estimated GAM parameters.

To calculate the uncertainty in estimates of $\hat{N}$, we can use posterior sampling (sometimes referred to as parametric bootstrapping), using the procedure outlined in (Wood, 2017, Section 7.2.6). We use the posterior distribution of β to generate possible abundance estimates, then calculate appropriate summary statistics. The following algorithm can be used:

1. Calculate the matrix X_p using the model and prediction grid.

2. For b ∈ {1, …, B}:

1. Generate new model parameters ${\beta }_b^{\ast }$.

2. Calculate the new linear predictor ${\eta }^{\ast }={\mathbf{X}}_p{\beta }_b^{\ast }$ for each prediction cell.

3. Calculate predictions on the response scale ${\hat{N}}_b^{\ast }=\mathbf{a}\exp {\mathbf{X}}_p{\beta }_b^{\ast }$.

4. Store ${\hat{N}}_b^{\ast }$ for this iteration.

3. Calculate empirical variance and mean of the B stored ${\hat{N}}_b^{\ast }$ values.

The posterior distribution of β is approximately multivariate normal with mean β and variance V_{β, θ}. To sample from this distribution we can directly sample using a multivariate normal random number generator, but may obtain better results by avoiding this assumption and using the Metropolis–Hastings algorithm to sample from the model posterior (such a sampler is provided as part of the package mgcv, see Supplementary Material B for an example of its use).

This method is extremely general. For example, we can extend this algorithm to calculate per-cell estimates by calculating per-cell abundances at 2.(c) above. A more complex example is when we include β_k in the model as in (3). In this case we are assuming that although the study area is the same, the densities are different for the different platforms. When making predictions we predict for each platform (k = 1, …, K) then sum these per-platform predictions.

Multiple surveys with uncertain detection on the trackline-fin whales

These data consist of observations of fin whales (Balaenoptera physalus) as part of NOAA’s Atlantic Marine Assessment Program for Protected Species. Data were collected during two distance sampling surveys: one shipboard (requiring adjustment for g(0)) and one aerial (requiring an adjustment for availability and g(0)). The left panel of Fig. 1 shows effort and detections. Details on field methods are available in Sigourney et al. (2020). Here we reproduce the analysis from that article using the DSM approach.

Shipboard observations (truncated at 6 km) were conducted in independent observer mode mark-recapture distance sampling with the model for g(0) only including distance as a covariate; a hazard-rate detection function was used, with Beaufort wind speed and a subjective measure of the effect of weather conditions on detectability as covariates. Observations from the aerial survey (truncated at 900 m) were modelled using a hazard-rate detection function with Beaufort wind speed as a covariate; a fixed g(0) and availability adjustment was applied based on results in Palka et al. (2017).

We specified the following DSM:

(5) $$\mathbb{E}(n_i)=a_i{p}_{i}g{(0)}_i{u}_i\exp \left[f(\text{DIST125})+f(\text{DEPTH})+f(\text{DIST2SHORE})+f(\text{SST})\right]$$where for segment i, n_i is the count (assumed to be Tweedie distributed), a_i is the area, u_i is the availability (one for shipboard segments, 0.37 for aerial segments), p_i is the probability of detection and g(0)_i is the correction for detection at zero distance (0.67 for aerial segments and based on the MRDS model for shipboard segments). f indicates a smooth constructed using thin plate regression splines with shrinkage (Marra & Wood, 2011), each with a maximum basis size of five. The covariates used were DIST125, distance to the 125 m isobath; DEPTH, depth; DIST2SHORE, distance to shore; SST, sea surface temperature. We used the variance propagation method of Bravington, Miller & Hedley, 2021 to propagate uncertainty from the two detection functions. Model checking gave reasonable results (see Supplementary Material A).

Abundance was calculated on a grid over shelf and shelf slope waters of the North East United States and South East Canada from Delaware to Nova Scotia. To estimate uncertainty in abundance, we can first take a Metropolis–Hastings sample from the posterior of our model (post variance propagation). The resulting samples capture uncertainty in the spatial model, aerial detection function, shipboard detection function and shipboard g(0). Assuming independence between the aerial estimates of g(0) and availability, we can use the fact that squared CVs add to combine uncertainty for these estimates in a way which is comparable to the estimate given in Sigourney et al. (2020).

Multiple platforms on one survey-fulmars

RRS Discovery conducted a survey for seabirds in the mid-Atlantic in June 2017 as part of cruise DY080. Figure 2 show the transects and observations (and Fig. S1 shows the study area in context of the North Atlantic). A modified version of the Eastern Canada Seabirds At Sea (ECSAS) protocol (Gjerdrum, Fifield & Wilhelm, 2012) was used while on effort, comprising a line transect survey for birds on the water and a strip transect survey for birds in flight. A single observer, located on one side of the bridge (varied according to conditions to avoid, e.g., glare) searched for all birds flying or sitting on the water within a 300 m wide strip to one side of the transect line, while a second person recorded bird sightings and sighting conditions. Birds were detected using the naked eye and identified using 10 × 40 or 8 × 40 binoculars. Birds on the water were detected in one of four distance bins: [0–50 m], (50–100 m], (100–200 m] and (200–300 m]. A “snapshot” method was used to record birds first detected in flight, flagging records of these birds if they were within a 300 m × 300 m box at the moment of the snapshot (Tasker et al., 1984). We selected observations of fulmars (Fulmarus glacialis) for analysis. Due to their distinctive plumage and flight behaviour, fulmars are generally easy to distinguish from other species at sea, but confusion can occur with large immature gulls. However, the latter are virtually absent from the study area during summer (Wakefield et al., 2021), so misidentification of fulmars was assumed to be negligible.

We therefore have two platforms. For fulmars on the water, we fitted a half-normal detection function, with precipitation (factor; yes/no) and visibility (in km) as covariates. For flying birds, we simply assume that detection is perfect out to 300 m (the width of the strip transect).

To explore the formulations defined above, we fitted models based on (2)–(4). Explanatory covariates were limited to spatial smooths of projected location (x, y) with variations on how the differing animal behaviour (“platform” in our terminology) was accounted for:

• (A) Based on (2), with linear predictor β₀ + f(x, y), the bivariate smooth of space used thin plate regression splines with shrinkage (Marra & Wood, 2011).

• (B) Based on (3), with linear predictor β_k + f(x, y) where β_k is an intercept depending on the behaviour. The spatial smooth was as above.

• (C) Based on (4), with linear predictor β₀ + f(x, y, k) where we have a factor-smooth interaction for the spatial effect. The factor-smooth model creates a smooth for each platform as deviations from a reference level.

All models were constructed so that the maximum basis size was 100. As above, the variance for the detection function was propagated using the method of Bravington, Miller & Hedley, 2021. A complete analysis of these data is available in Wakefield et al. (2021).

We compared our models using the following techniques:

1. Residual model checking procedures outlined in e.g., (Wood, 2017, Chapter 7) can be used to assess the models and potentially remove from consideration models that have violated assumptions.

2. AIC scores can be used to compare models, as we might usually do for GAMs. Since here we use the same detection model in each case we need not include the detection function AIC in our considerations (but if different detection functions were used we would).

3. Goodness-of-fit can be assessed in terms of the number of observed vs. predicted animals swimming or flying. Aggregating at the level of some binned covariate is important as the smoother will (by its nature) tend to predict small values where exact zeros were observed. We can aggregate by any covariate (whether or not it is included in the model), here we use the platform covariate to assess our fitted models since we are trying to address differences in platform.

4. For models that include additional terms to account for platform (models B and C), plotting the difference surface between per-platform predictions spatially. We expect that if that extra complexity is not needed we would not see a difference in predicted densities.