1 Introduction
Line transect surveys are a standard method to survey both marine and terrestrial populations; distance sampling methods are then used to estimate abundance while accounting for detectability (Buckland et al., 2001). There is particular interest in spatiallyexplicit (modelbased) estimates of cetacean density for management and mitigation purposes (International Whaling Commission, 2012; Hammond et al., 2013; Roberts et al., 2016) where it is important that abundances (and their variances) are calculated over arbitrary areas. Spatial models for distance sampling data must take into account detectability data, usually via an offset in the model. Uncertainty about the (fixed) offset must be incorporated and the covariance between the detectability and spatial parts of the model needs to be estimated. This can be particularly problematic when sighting conditions change as a function of space, in this case the detection process and distribution processes are confounded. For example, during cetacean surveys weather changes in space as the survey progresses, so sightability and habitat covariates are both functions of space. In this paper we show how the detection function uncertainty can be accommodated painlessly within standard software.
There are many methods for estimating spatiallyexplicit abundance estimates from distance sampling data, here we consider density surface models (DSMs; Hedley and Buckland, 2004; Miller et al., 2013): twostage models that first use a detection function to model detectability, followed by a spatial model (in our case a generalized additive model; GAM, e.g. Wood, 2017) to capture the animals density as a function of environmental covariates such as sea surface temperature or bathymetry. Detectability is modelled using (for example) methods described in Marques and Buckland (2003), where covariates such as weather or observer are used to estimate the effect of observation conditions on detectability.
This twostage approach is convenient because it uses wellestablished software and diagnostic tests for each part of the problem. This does however complicate variance estimation, as uncertainty from the detection function must be propagated through to the spatial model. Joint fitting of the detection function and spatial model is not possible in most GAM packages. In this paper we fit a detection model, then fit the GAM including a quadratic approximation to the detection function loglikelihood centred on the point estimate which enters the linear predictor as a zeromean normal random effect. We therefore effectively (re)fit both models at once. This approach allows us to use welltested model checking and diagnostics used first for detection function fitting and then for GAMs. We also propose a diagnostic based on our variance propagation method that can be used to assess model specification issues.
Following a summary of DSMs and notation in Section 2, we present our new formulation in Section 3, including variance computation and diagnostics. Section 4 comments on problems with existing approaches to variance propagation in DSMs. In Section 5, we extend the formulation to cover DSMs where group size varies spatially and affects detectability (a common situation with whales and dolphins). Section 6 gives examples of the variance propagation method and the group size model. Some discussion is given in Section 7, including comparison with single stage methods such as Johnson et al. (2009) and Yuan et al. (2017).
2 Density surface models
We focus on modelling line transect distance sampling data: observers move along a set of survey lines, counting (groups of) animals, recording distances from the line to the observed groups (or their cues, such as blows for cetaceans), and the size of each detected group.
To fully describe the DSMs in this paper, we distinguish four different classes of variable.

Density covariates, , vary in space and potentially affect local animal abundance: e.g., latitude and depth; they are required for prediction and fitting, and are assumed known across the entire region of interest.

Effort covariate(s), (or ), affect detection probability: e.g., sea conditions measured on the Beaufort scale, observer name; assumed known along each transect, but not necessarily in unsurveyed areas.

Individual covariates, , that affect detection probability and are a persistent property of each group (independent of whether the group is observed or not) during its window of observability: e.g., size (number of animals), and perhaps behaviour.
is assumed known for each observed group (see Discussion). The random variable
varies from one group to the next, and its statistical distribution may vary spatially. may have a direct effect on abundance (via the mean group size), as well as on detection probability, in which case it is also necessary to estimate certain properties of such as its local mean. 
Observation variables, , which are random properties of one observation on one group: e.g., perpendicular distance between the group and the trackline. In certain settings, may contain other elements. For example, in a multiobserverplatform survey (e.g., MRDS; Buckland et al., 2004), might also include which of the active observers saw the group; in a cuebased setting, might include the bearing between sighting and trackline.
These classes are assumed to be mutually exclusive; overlaps can lead to fundamental problems for distance sampling which we do not address here (such as nonuniform animal distribution within the sample unit; e.g., Marques et al., 2012). The distinction between individual and effort covariates is often glossed over but they have rather different implications for abundance estimation (see below).
In the first stage of DSM, the detection function , which involves unknown parameters as well as and , describes the probability of making an observation at some distance . The parameters are usually estimated by maximizing this loglikelihood across observations :
(1) 
where is the transect containing sighting . Here is the overall detection probability for a group, defined by
where is the distribution function of . In standard distance sampling where consists only of perpendicular distance, is uniform between and some fixed truncation distance, beyond which observations are discarded. This formulation encompasses a wide range of models, including multiple covariate distance sampling (MCDS; Marques and Buckland, 2003) with and , multiobserver markrecapture distance sampling (MRDS; Borchers et al., 1998), and cuebased “hazard probability models” (Skaug and Schweder, 1999).
The second part of DSM models the local count of observations via a GAM to capture spatial variation in animal density. This allows us both to estimate abundance within any subregion of interest, and to compensate as far as possible for uneven survey coverage (whether by design, or by virtue of field logistics and weather conditions). Since transects are generally very long in comparison to their width and therefore contain a range of density and density covariate values, we divide transects into smaller segments, which are the sample units for GAM; environmental covariates are assumed not to change much within each segment. The relationship between counts, , per segment, , and density covariates, is modelled as an additive combination of smooth functions:
(2) 
where each segment is of area ), and follows some count distribution such as quasiPoisson, Tweedie or negative binomial and we assume a link. The are smooth functions, represented by a basis expansion (i.e., for some basis functions ) penalized by a sum of quadratic penalties;
is an intercept term, included in parameter vector
; is a vector of smoothing (hyper)parameters which control the wiggliness of the. We take a Bayesian interpretation of GAMs (Wood, 2017), in which controls the variance of a multivariate improper Gaussian prior (Marra and Wood, 2011):with scale parameter , smoothing parameters and penalty matrices ( indicates pseudoinverse). We estimate via REML (Wood, 2011)
, an empirical Bayes procedure. Fully Bayesian approaches, placing hyperpriors on
are also possible. Ignoring group size for now we write (we return to group size in Section 5). We may think of the product as the effective area of segment , analogous to the effective strip (half) width (Buckland et al., 2001).We are interested in the uncertainty of a predicted abundance estimate, . We assume below that we have created some prediction grid with all density covariates available for each cell in the grid. Abundance is predicted for each cell, and summed for an overall abundance, , over some region of interest which may not be the entire surveyed area. Although does not appear explicitly in the prediction:
the GAM offsets clearly do affect , so it is important to account somehow for detection probability uncertainty. (2) assumes the offset is fixed, so extra steps are required.
3 Variance propagation for Density Surface Models
Let be the true probability of detection in segment and omit , therefore assuming that there are no individuallevel covariates (e.g., that group size is always 1) for now (see Section 5). If is the true (unknown) value of , and is its MLE, we use the shorthand and when the dependence is clear. The expected number of encounters in segment is where is the underlying density.
Given , we can rewrite (2) on the link scale as:
(3) 
is the (known) ^{th} row of the design matrix, i.e., the values of the basis functions in segment , so and is an offset. The complication is that we only have an estimate of . To tackle this, we first rewrite the linear predictor as
and then take a Taylor series expansion of about :
(4) 
By defining the vectors and , we can rewrite 3 as
(5) 
Asymptotic arguments suggest that the term can be neglected for reasonable sample sizes. Since estimation of spatial models requires a large amount of data, the requirements for the quadratic approximation are probably met in most cases (if the spatial model is reasonable).
We have approximately that where the covariance matrix is calculated as the negative inverse Hessian of 1. In other words, the “posterior distribution” of from fitting the detection function now becomes a prior distribution for . To first order, then plays the same structural role in (5) as the basis coefficients . The design matrix for ( in (5)) is obtained by differentiating the logdetection probabilities, with respect to at . Derivatives do not need to be accurate, and need only be computed once, so simple 3point numerical differentiation is perfectly adequate. should be readily available (via the Hessian) from distance sampling stage regardless of the complexity of the model. This method can be applied automatically to almost any distance sampling setup. Simultaneous estimates and can be obtained from standard GAM fitting software. Posterior inferences about (therefore and abundance) automatically propagate the uncertainty from fitting the detection function.
The only technical difference from fitting a standard GAM, is that is usually unknown and has to be estimated (i.e., the prior on has known covariance, but unknown scale), whereas the prior on is completely determined from the detection function fitting (i.e., in effect , where is the scale parameter). This setup cannot be specified directly in the R package mgcv because of implementation details (at least up to version 1.8; it may be possible within other GAM implementations), unless is fixed (as opposed to estimated). With persegment count data, a fixed (of 1) would correspond to a purely Poisson response, which might be appropriate if maximum count per segment is only 1 or 2. In our experience, better fits can usually be obtained by estimating and using a Tweedie response distribution. In order to implement (5) for a general response distribution using mgcv (Wood, 2011), we therefore use a onedimensional search over to maximize the marginal REML. At each iteration, given the working value , we refit the GAM fixing and . Speed can be improved by reusing some of the setup computations (design matrices, etc) at each iteration.
3.1 Diagnostics and interpretation of
Posterior modes are available after fitting, and may help with diagnosing inconsistencies between the detection function and spatial model. For example, if weather is systematically worse in some parts of the survey region, then both and will contribute to the expected pattern of sightings, and the two sets of parameters will be confounded. Fortunately, because both and are penalized (i.e., partly constrained by their prior distributions), the spatial model has limited ability to affect , but it is nevertheless worth checking whether there is any such tension between the two parts of the model. This can be checked in a couple of ways. The first is to compare the inferred spatial distribution and abundance from fitting (3) with the “naive” estimates where detection uncertainty is ignored and the offset is treated is exact. The second is to check whether the detection probabilities (by covariate level) would be substantially changed by fitting the spatial model; in other words, whether
is close enough to zero given its prior distribution, or, perhaps more usefully, whether the overall detectability by covariate level has changed. Since the fitted spatial model still includes the information from the first stage, any shift of more than about 1 standard deviations (based on the covariance from the detection function stage) might merit investigation.
As a general diagnostic tool for density surface models, we have found it useful to compare total observed and expected numbers of sightings, grouped by detection covariates (e.g., Beaufort). This can be helpful in diagnosing detection function problems, e.g., failure of certain detectability at zero distance under poor weather conditions, as well as failures of the spatial model (e.g., an abrupt change in density). In addition one could also use standard detection function model checking (e.g., quantilequantile plots) with the adjusted parameters.
3.2 Calculating
Given a model that includes detection function uncertainty, we now discuss two methods to calculate the variance of .
Posterior simulation. The posterior for given data and smoothing parameters , are distributed as ; where is the covariance matrix for the GAM coefficients Wood (2017, Sections 5.8 & 6.9.3). We form the prediction matrix, , which maps model coefficients to values of the linear predictor for the prediction data so (Wood, 2017, Section 6.10). The following algorithm can be used:

For :

Simulate from , to obtain .

Calculate predicted abundance, (where is a row vector of areas for the prediction cells).


Calculate the empirical variance or percentiles of the s.
In practice does not have to be particularly large. Marra et al. (2011) achieve reasonable results with .
Analytically. One can avoid the computational burden of posterior simulation (which can be large for complex models) by calculating:
(6) 
(the delta method; Seber 1982, Section 1.3.3) where derivatives are evaluated at the estimated values of the model parameters.
4 Previous methods for estimating uncertainty in Density Surface Models
Several approaches have been suggested to combine detection function and spatial model predicted abundance uncertainties; we review them briefly here. We need to estimate the following:
where here is a random variable for the probability of detection and the subscripts indicate the expectation/variance taken over that variable. The first part of this can be derived from GAM theory and is detailed below, the second is more tricky.
Assuming independence. When is the same for all observations, then , so and are independent. The total variance of the abundance estimate can be calculated by combining the GAM variance estimate with the variance of the probability of detection summing the squared coefficients of variation () (Goodman, 1960). Hence
When there are not covariates in the detection function we calculate:
(7) 
In the case where detectability is a function of covariates, it is hard to justify the use of the CV decomposition as there are correlations between the spatial distribution and the covariates affect detectability. Program Distance (Thomas et al., 2010) uses the HorvitzThompsoncorrected counts per segment as the response (thus removing the detectability from the right hand side of (2)). Variance is then calculated by taking the probability of detection averaged over the observations by first calculating the HorvitzThompson estimate of the abundance in the covered area (, where is group size of the observation and is the probability of detecting that group) then using that to calculate the implied average detectability, had the analysis not contained covariates (, where is the number of observed groups). The numerical derivatives of with respect to can then be used in (7) to derive a variance for this probability of detection, averaged over the observations.
We do not recommend this approach, for several reasons. Transforming the response through multiplication by a random variable breaks the meanvariance and independence assumptions of the GAM, so that the computed is invalid (unless there are no covariates in the detection function). Additionally, there is no coherent way to generalize the formula to smallarea predictions — the effort covariates within a small area will not have the same range as those in the larger survey area (i.e., weather conditions will not be homogenous throughout the survey area). Hence, the uncertainty that applies to the overall is usually not the appropriate uncertainty to apply to a small area where observing conditions may be untypical. While some ad hoc correction could perhaps be contrived, it is not obvious how to do so reliably in general.
The bootstrap. Bootstraps are sometimes seen as an attractive alternative to deal with all aspects of variance in DSMs, bypassing awkward analytical calculations. There are many different ways to bootstrap a DSM; Hedley and Buckland (2004) describe two possible implementations (one parametric, one nonparametric), which are not easy to choose between and which do not necessarily give similar answers. Ignoring computational time issues the first practical difficulty in setting up a “good” nonparametric bootstrap for a DSM is sampling units, independent of the fitted model: in situations where there is more than one survey occasion, it is unclear exactly how one decides on the units used for resampling.
The second, more substantial, issue is the fundamental statistical problem with combining smoothers (which essentially come from a Bayesian/random effects paradigm) with bootstraps (a deeply frequentist concept). The problem does not seem to be wellknown in the statistical ecology literature, so we give an explanation here. The basic problem is that (most) bootstraps use only the posterior modes of random effects, thus omitting a key part of the posterior uncertainty. To see this, consider a simple “spatial model” where the region is divided into blocks, each with its own independent random effect, and a bootstrap that generates new data at each original observation/transect, either parametrically or nonparametrically. If one of the blocks is unsampled in the original data, it will be unsampled in every realization too, and the “spatial model” simply sets the point estimate of that random effect to zero in every bootstrap realization; hence a bootstrap will ascribe zero uncertainty about the density in that block. The correct inference would of course be for the random effect to retain its prior variance. This phenomenon has been wellknown in statistics since at least Laird and Louis (1987) (see also the discussants), who coined the term “naïve bootstrap” for such procedures that ignore the point estimate shrinkage inevitable in mixed or random effect models (fixed effect models are not susceptible in the same way). They proposed some parametric modifications (“type II” and “type III” bootstraps) that are more sensible in the limited situations they consider. However, the underlying theory is complex (Carlin and Gelfand, 1991; Carlin and Louis, 2008) and it is far from clear whether simple yet reliable bootstraps can be devised for complicated multistage random effect situations like DSMs.
The GAM formulation, whereby each random effect spreads its influence across a wide area and the posterior distributions of the are not independent, disguises the naïve bootstrap issue but does not fully remove it, especially in line transect settings. Figure 1 shows a simple unidimensional Poisson process, sampled at either end but not in the middle (rug plot). Bootstrap replicates (shown in light grey, of which there are 500) largely fail to capture our uncertainty in the unsampled middle area. The analytical estimate (dark grey band) illustrates how little we know about the unsampled area.
The above does not imply that simple or indeed complicated bootstraps will never give reliable results in line transect density surface modelling; given plenty of observations and good, even coverage, many approaches to inference will give similar and good results. However, it is sometimes not obvious whether this holds for a specific dataset, nor what to do bootstrapwise if not. Instead, the (empirical) Bayes framework of GAMs offers a coherent and generalpurpose way to capture uncertainty which we explore next. See Wood (2017, Section 6.10.3) for further comment.
5 A new model for group size
Our variance propagation method addresses when detectability depends only on effort covariates, not for individual covariates such as group size. Incorporating individual covariates in the detection function is not problematic but it is not obvious how to allow for these different detection probabilities in the GAM. Further, it is not obvious how to combine predictions of different group sizes since average group sizes may vary spatially (Ferguson et al., 2006).
One approach is to use HorvitzThompson estimates of abundance per segment as the response in the model, correcting for detectability in the response rather than incorporating it via the offset of the model (Hedley and Buckland, 2004), but this does not allow variance propagation. One could fit a separate spatial model to subsets of the data for each group size, though it seems inefficient to not share information between subsets of the data. To account for group size variation in space and in detectability we propose a new model that allows us to propagate variance in the manner described above.
We form categories of group sizes, denoted , where groups within each category have similar detectabilities and fit a detection function incorporating these group size categories. We then fit a GAM to an fold replicate of the dataset, with the response in the replicate of the segment being , the number of groups in category that were seen in that segment. (The total number of observations is unchanged; each observation is allocated to just one of the “replicates”). Group size category (as a factor) is included as an explanatory variable, and smooths are modified to allow similar variations in density of groups with different sizes. There are no extra assumptions in this formulation from the model in Section 3, except to assume that the numbers of groups of different size categories in a given segment are independent, given the underlying density (which is allowed to vary with group size).
Factorsmooth interactions
We extend (2) to include multiple smooths of space which correspond to different categorizations of group size, so our model is:
for where is the number of observed groups in group class in segment and is the spatial smooth (where is a spatial coordinate) for group size class . Smoothers like are referred to as factorsmooth interactions (Wood, 2017; Pedersen et al., 2018). are any other smooths (of covariates ). For clarity we make the dependence on group size class explicit: , i.e., the probability of detection given segmentlevel detection covariates and group size . In practice we duplicate the full segment data as many times as there are group size categories (so our model now has (number of segments) rows of data) with a response of the number of groups in that category that occur in that segment.
There are a number of different possible forms for . These vary in two main ways: (1) do levels share a smoothing parameter, or have separate ones? (2) do smooths tend toward a “global” smooth that dictates a general spatial effect. Here we adopt the “fs” basis in mgcv which can be thought of as a nonlinear version of a random slopes model: smooths are generated for each factor level with smooths defined as deviations from a reference level. Each of the levels has the same smoothing parameter; this is appealing as we might expect that the spatial smooths for each group size are similar but there might be some process that generates larger groups in certain places (e.g., large prey aggregations attracting large groups of animals). This approach is easily extended to nonspatial smoothers (e.g., depth).
Abundance and uncertainty estimation with group size smooths
Abundance is estimated by summing over the predictions for each group size category () and weighting them by the corresponding mean group size ():
We can find (where is the mean group size) from the variance propagation procedure above, but we need , which we can obtain from the law of total variance:
(8) 
where is the variance within category , estimated by its empirical variance. The effect of on should be small (because categories are narrow, and mean must lie within category), and also should not vary much spatially, so further spatial adjustment to that variance component is not required. Estimates are of course correlated as they are fitted within a single GAM, so is a covariance matrix.
6 Examples
We reanalyse line transect data for two species from the SCANSII surveys in European Atlantic waters. For survey details and a serious analysis see Hammond et al. (2013). Here we assume certain detection on the trackline and ignore island/coastline effects in the spatial model. We computed a diagnostic for both models comparing observed versus expected counts at the observed levels of observationlevel covariates (in addition to usual distance sampling/GAM diagnostics). Models were fitted using the dsm package in R, using the dsm_varprop function to calculate the variance of .
6.1 Minke whales
Minke whale (Balaenoptera acutorostrata) shipboard survey data (from the North Sea; strata U and V in Hammond et al. 2013), included Beaufort as a covariate in a hazardrate detection function (recorded at each segment, binned as [0, 1], (1, 2], (2, 4], only including segments where sea state ). The following spatial model was then fitted:
where , constrained to have a minimum power parameter of 1.2 since values in the range 1–1.2 make little practical difference to variances but can lead to numerical problems. are projected coordinates. Distances were truncated at 870m for detection function fitting, as in Hammond et al. (2013); we ignore group size as most group sizes were predominantly 1 (total 68 observations, with 3 size 2 and 1 size 5, all others singletons). Table 1 shows the number of predicted minke whales per kilometre of transect from the model given Beaufort — minke density is correlated with weather conditions, so we need to account for this when calculating variance.
Table 2 compares estimated detectability at different Beaufort levels from the detection function with those “corrected” during the fitting of the variance propagation model; we see there is little difference in the probability of detection after refitting. 0.2465 using variance propagation method versus 0.2643 when using (7). Our estimates are lower, indicating there is some negative covariance between density and detectability. increases in the most severe sea state from the middle level (last row of Table 2), so this feature is perhaps driving covariance between the spatial model and detection function.
Beaufort  Predicted abundance per km  Segments  

[0,1]  0.026  0.508  255 
(1,2]  0.008  0.216  316 
(2,4]  0.007  0.247  386 
Beaufort  

[0,1]  0.508  0.115  0.482 
(1,2]  0.216  0.078  0.203 
(2,4]  0.247  0.088  0.304 
, its standard error and the “corrected” detectability from the variance propagation model denoted
for the minke whale data. Beaufort is treated as a factor in the analysis, hence detectability is nonmonotonic with increasing sea state.[0,1]  (1,2]  (2,4]  

Observed  41.00  14.00  20.00 
Expected  42.58  15.28  14.59 
6.2 Harbour porpoise
Harbour porpoise (Phocoena phocoena) aerial survey data from SCANSII in the Irish Sea, coastal Irish waters and Western coastal Scotland (strata O, R and N, respectively) was used to illustrate the group size model. Aerial surveys were preferred as error in group size is likely less of an issue than in shipboard surveys for harbour porpoise (Phil Hammond, Debi Palka, personal communication, November 2017). Three group size bins were formed: size 1 (131 observations), 2 (35 observations) and 35 (14 observations). We fitted the following model:
where is the detectability (estimated from a hazardrate detection function) in segment for group size class (), when the Beaufort level was (same as for minke whales). where the power parameter was constrained to be . Each
had a maximum basis size of 20 (total basis size of 60); the fitted model had a total effective degrees of freedom of 20.47 for
. Table 4 shows observed vs expected counts by Beaufort — there is some misfit at the highest state (perhaps due to our assumption certain detection at zero distance; Hammond et al. estimated this as 0.45 for Beaufort 01 and 0.31 in Beaufort 23). GAM checking showed reasonable fit to the data (low deviance explained of 26.4% is not uncommon with DSMs). Plots of the pergroup size bin predictions and the combined prediction are given in Figure 2. Using (7), the CV of abundance was estimated to be 0.0236, when our new variance propagation method was used the CV was estimated as 0.0965. So variance calculated assuming independence (via (7)) underestimates uncertainty in the case were group size (and therefore detectability) vary in space. Only a small piece of code implementing (8) was required in addition to the dsm_varprop function.[0,1]  (1,2]  (2,4]  

Observed  98  36.00  30.00 
Expected  96  34.32  37.14 
7 Discussion
We have shown how to (i) propagate uncertainty from detectability models to the spatial models for a particular class of detection function (i.e., those without individuallevel covariates) and (ii)
include group size as a covariate in the detection function while still being able to propagate uncertainty and address spatial variation in group size. These methods are implemented in the dsm package for R but can be implemented in any standard GAM fitting software. We also have reiterated the issues with using bootstraps for estimating uncertainty in random effects models, a model class that includes DSMs. It is straightforward to apply our group size approach to covariates which affect detectability and vary in space, but do not directly affect abundance; e.g., behaviour: if feeding groups are less conspicuous (by being underwater) than resting groups (at the surface) of similar size, and feeding is more prevalent in some parts of the surveyed area. A major advantage of this approach over simple (or complex) stratification schemes is that we are now sharing information between the levels of our categorized variable, this makes the results less sensitive to over specifying the number of categories, as the model will shrink back towards the simpler model in the absence of strongly informative data.
We have assumed that all variables are measured without much error. Measurement error for individuallevel covariates such as group size can be a serious problem in distance sampling (Hodgson et al., 2017) — distance between observer and group can affect not just detectability, but also the extent of group size error. If group size varies spatially, it is hard to see how to separate the spatial modelling stage from the distance sampling stage. A full discussion is beyond the scope of this paper, but we suspect that speciallydesigned observation protocols and bespoke analyses may be the only way to tackle such thorny cases.
Allinone fitting of both detection and spatial models is also possible (e.g., Johnson et al., 2009; Yuan et al., 2017). If models are specified correctly the allinone approach could be slightly more efficient, but adding the full detection function likelihood to the spatial model likely would not make much difference in addition to our quadratic approximation used here. Nevertheless, our own preference is for the twostage approach, mainly because in our experience the careful fitting of detection functions is a complicated business which can require substantial model exploration and as few as possible “distractions” (such as simultaneously worrying about the spatial model). The twostage process allows any form of detection function to be used, without having to make deep modifications to software. In summary, if one knew one had the correct model to begin with, onestage fitting would be more efficient, but this is never the case in practice. It is valuable to check for any tension or confounding between the detection function and density surface parts of the model, which can occur if there are largescale variations in sighting conditions across the survey region (Section 3.1), this is readily diagnosed in a twostage model. Although this does not appear to lead to problems in the datasets we have analysed with the software described in this paper, we have come across it in other variants of line transectbased spatial models with different datasets. It may not be so easy to detect partial confounding when using allinone frameworks.
Finally, we note that the approach outlined here (of using a firststage estimate as a prior for a secondestimate, and propagating variance appropriately) is quite general and is comparable to standard sequential Bayesian approaches to socalled “integrated data models”. The firststage model need not be a detection function (or such a disparate model) but instead could be from another GAM (or other latent Gaussian model). Again, this allows us to ensure that firststage models are correct before moving to more complex modelling. Modelling need not only be twostage and could extend to multistage models.
Supplementary Materials
R code implementing the group size model in Section 5 along with data and code for the minke whale and harbour porpoise analyses in Section 6 are available at https://github.com/dill/varpropsuppmaterials.
Acknowledgements
The authors thank Natalie Kelly, Jason Roberts, Eric Rexstad, Phil Hammond, Steve Buckland and Len Thomas for useful discussions, Devin Johnson for the suggestion of this as a general method, and Simon Wood for continued development of mgcv and GAM theory. Data are from the SCANSII project supported by the EU LIFE Nature programme (project LIFE04NAT/GB/000245) and governments of range states: Belgium, Denmark, France, Germany, Ireland, Netherlands, Norway, Poland, Portugal, Spain, Sweden and UK. This work was funded by OPNAV N45 and the SURTASS LFA Settlement Agreement, and being managed by the U.S. Navy’s Living Marine Resources program under Contract No. N3943017C1982, US Navy, Chief of Naval Operations (Code N45), grant number N002441010057 and the International Whaling Commission.
References
 Borchers et al. [1998] D. L. Borchers, W. Zucchini, and R. M. Fewster. MarkRecapture Models for Line Transect Surveys. Biometrics, 54(4):1207, 1998.
 Buckland et al. [2001] S. T. Buckland, D. R. Anderson, K. P. Burnham, D. L. Borchers, and L. Thomas. Introduction to Distance Sampling. Estimating Abundance of Biological Populations. Oxford University Press, Oxford, UK, 2001.
 Buckland et al. [2004] S. T. Buckland, D. R. Anderson, K. P. Burnham, J. L. Laake, D. L. Borchers, and L. Thomas. Advanced Distance Sampling. Estimating abundance of biological populations. Oxford University Press, Oxford, UK, Aug. 2004.

Carlin and Gelfand [1991]
B. P. Carlin and A. E. Gelfand.
A sample reuse method for accurate parametric empirical Bayes confidence intervals.
Journal of the Royal Statistical Society: Series B, 1991.  Carlin and Louis [2008] B. P. Carlin and T. A. Louis. Bayesian Methods for Data Analysis, Third Edition. CRC Press, June 2008.
 Ferguson et al. [2006] M. C. Ferguson, J. Barlow, P. Fiedler, S. B. Reilly, and T. Gerrodette. Spatial models of delphinid (family Delphinidae) encounter rate and group size in the eastern tropical Pacific Ocean. Ecological Modelling, 193(34):645–662, Mar. 2006.
 Goodman [1960] L. A. Goodman. On the Exact Variance of Products. Journal of the American Statistical Association, 55(292):708, 1960.
 Hammond et al. [2013] P. S. Hammond, K. Macleod, P. Berggren, D. L. Borchers, L. Burt, A. Cañadas, et al. Cetacean abundance and distribution in European Atlantic shelf waters to inform conservation and management. Biological Conservation, 164(C):107–122, Aug. 2013.
 Hedley and Buckland [2004] S. L. Hedley and S. T. Buckland. Spatial models for line transect sampling. Journal of Agricultural, Biological, and Environmental Statistics, 9(2):181–199, June 2004.
 Hodgson et al. [2017] A. Hodgson, D. Peel, and N. Kelly. Unmanned aerial vehicles for surveying marine fauna: assessing detection probability. Ecological applications : a publication of the Ecological Society of America, 27(4):1253–1267, 2017.
 International Whaling Commission [2012] International Whaling Commission. Requirements and Guidelines for Conducting Surveys and Analysing Data within the Revised Management Scheme. Journal of Cetacean Research and Management, 13(SUPPL):509–517, 2012.
 Johnson et al. [2009] D. S. Johnson, J. L. Laake, and J. M. Ver Hoef. A ModelBased Approach for Making Ecological Inference from Distance Sampling Data. Biometrics, 66(1):310–318, 2009.
 Laird and Louis [1987] N. M. Laird and T. A. Louis. Empirical Bayes Confidence Intervals Based on Bootstrap Samples. Journal of the American Statistical Association, 82(399):739–750, Sept. 1987.
 Marques and Buckland [2003] F. F. C. Marques and S. T. Buckland. Incorporating covariates into standard line transect analyses. Biometrics, 59(4):924–935, Dec. 2003.
 Marques et al. [2012] T. A. Marques, S. T. Buckland, R. Bispo, and B. Howland. Accounting for animal density gradients using independent information in distance sampling surveys. Statistical Methods & Applications, 22(1):67–80, Nov. 2012.
 Marra and Wood [2011] G. Marra and S. N. Wood. Practical variable selection for generalized additive models. Computational Statistics and Data Analysis, 55(7):2372–2387, July 2011.
 Marra et al. [2011] G. Marra, D. L. Miller, and L. Zanin. Modelling the spatiotemporal distribution of the incidence of resident foreign population. Statistica Neerlandica, 66(2):133–160, Nov. 2011.
 Miller et al. [2013] D. L. Miller, M. L. Burt, E. A. Rexstad, and L. Thomas. Spatial models for distance sampling data: recent developments and future directions. Methods in Ecology and Evolution, 4(11):1001–1010, Aug. 2013.
 Pedersen et al. [2018] E. J. Pedersen, D. L. Miller, G. L. Simpson, and N. Ross. Hierarchical generalized additive models: an introduction with mgcv. PeerJ Preprints, 6:e27320v1, Nov. 2018. ISSN 21679843. doi: 10.7287/peerj.preprints.27320v1. URL https://doi.org/10.7287/peerj.preprints.27320v1.
 Roberts et al. [2016] J. J. Roberts, B. D. Best, L. Mannocci, E. Fujioka, P. N. Halpin, D. L. Palka, L. P. Garrison, K. D. Mullin, T. V. N. Cole, C. B. Khan, W. A. McLellan, D. A. Pabst, and G. G. Lockhart. Habitatbased cetacean density models for the U.S. Atlantic and Gulf of Mexico. Scientific Reports, pages 1–12, Mar. 2016.
 Seber [1982] G. A. F. Seber. The Estimation of Animal Abundance and Related Parameters. Macmillan, New York, 1982.
 Skaug and Schweder [1999] H. J. Skaug and T. Schweder. Hazard models for line transect surveys with independent observers. Biometrics, 55(1):29–36, Mar. 1999.
 Thomas et al. [2010] L. Thomas, S. T. Buckland, E. A. Rexstad, J. L. Laake, S. Strindberg, S. L. Hedley, et al. Distance software: design and analysis of distance sampling surveys for estimating population size. Journal of Applied Ecology, 47(1):5–14, Feb. 2010.
 Wood [2011] S. N. Wood. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(1):3–36, 2011.
 Wood [2017] S. N. Wood. Generalized Additive Models: An Introduction with R, Second Edition. Chapman & Hall/CRC Texts in Statistical Science. CRC Press, 2017.
 Yuan et al. [2017] Y. Yuan, F. E. Bachl, F. Lindgren, D. L. Borchers, I. J. B., S. T. Buckland, et al. Point process models for spatiotemporal distance sampling data from a largescale survey of blue whales. The Annals of Applied Statistics, 11(4):2270–2297, 2017.
Comments
There are no comments yet.