Estimating seal pup production in the Greenland Sea using Bayesian hierarchical modeling
EESTIMATING SEAL PUP PRODUCTION IN THE GREENLAND SEAUSING BAYESIAN HIERARCHICAL MODELING
MARTIN JULLUM , THORDIS THORARINSDOTTIR , AND FABIAN E. BACHL Abstract.
The Greenland Sea is an important breeding ground for harp and hooded seals.Estimates of the annual seal pup production are critical factors in the abundance estimationneeded for management of the species. These estimates are usually based on counts fromaerial photographic surveys. However, only a minor part of the whelping region can be pho-tographed, due to its large extent. To estimate the total seal pup production, we propose aBayesian hierarchical modeling approach motivated by viewing the seal pup appearances as arealization of a log-Gaussian Cox process using covariate information from satellite imageryas a proxy for ice thickness. For inference, we utilize the stochastic partial differential equa-tion (SPDE) module of the integrated nested Laplace approximation (INLA) framework. Ina case study using survey data from 2012, we compare our results with existing methodologyin a comprehensive cross-validation study. The results of the study indicate that our methodimproves local estimation performance, and that the increased prediction uncertainty of ourmethod is required to obtain calibrated count predictions. This suggests that the samplingdensity of the survey design may not be sufficient to obtain reliable estimates of the seal pupproduction. Introduction
Three stocks of harp seals (Pagophilus groenlandicus) and two (possibly three) stocks ofhooded seals (Cystophora cristata) inhabit the North Atlantic Ocean where they have beenharvested for centuries (Sergeant, 1974, 1991; Kovacs & Lavigne, 1986). Monitoring the abun-dance of seals is vital for controlling the biodiversity in the region. State-of-the-art seal pop-ulation models are dynamically built based on historical catch data (Øig˚ard et al., 2014a,b).The main ingredient in these models is the total pup production in a given year which needsto be quantified based on on-site observational data since other quantification methods basedon catch-at-age and mark-recapture data etc. are considered unreliable (ICES, 2014). Thewhelping regions in the North Atlantic typically cover several thousand square kilometers sothat the total pup production needs to be estimated based on observations from a minor partof the region. Both the estimated total pup production and the associated uncertainty are thenused as input in dynamic population models (Øig˚ard et al., 2010).The observational data consists of seal pup counts obtained by manual counting on pho-tographs. The photographs stem from an aerial photographic survey conducted by flying alongtransects sparsely covering the whelping region. The survey methodology is discussed in moredetail in Section 2. The traditional method for estimating the total pup production based onsuch count data is that of Kingsley et al. (1985) which assumes a homogeneous dispersion ofseals across the entire whelping region. Salberg et al. (2009) propose a generalized additivemodeling (GAM) approach (Hastie & Tibshirani, 1990), assuming the counts follow a negativebinomial distribution and taking the spatial location of the counts into account. For data thatis close to homogeneous, the negative binomial GAM approach and the Kingsley method yieldsimilar estimates. However, the Kingsley method may possess a positive bias when the spatialdistribution of the pups is clustered (Salberg et al., 2008; Øig˚ard et al., 2010). Additionally,the GAM method produces much smaller uncertainty bounds than the homogeneous Kingsleyapproach.
Norwegian Computing Center, Oslo, Norway. School of Mathematics, University of Edinburgh, Edinburgh, Scotland.
E-mail addresses : [email protected], [email protected], [email protected] . a r X i v : . [ s t a t . A P ] J u l STIMATING SEAL PUP PRODUCTION 2
In this paper, we propose a new method for estimating the total seal pup production. Weview the seal pup appearances as a spatial point process (Møller & Waagepetersen, 2003)and model the point pattern of the seal pups as a log-Gaussian Cox process (LGCP; Mølleret al., 1998) with a spatial latent field which also allows additional covariate information to beaccounted for. In a Bayesian formulation with priors on the model parameters, the seal pupproduction estimate is represented by the posterior predictive distribution found by integratingthe posterior distribution over the spatial domain of the whelping region, instead of a singlepoint estimate accompanied with a variance estimate. This Bayesian hierarchical model can befitted by utilizing the stochastic partial differential equation (SPDE) approach of the integratednested Laplace approximation method (INLA; Lindgren et al., 2011; Rue et al., 2009). Thefinal posterior predictive distribution can subsequently be computed from this fitted model by asampling approach. Although more traditional Markov Chain Monte Carlo (MCMC) methodsin theory could be used to arrive at the same posterior, application of INLA allows results tobe produced magnitudes faster, at a negligible cost in terms of accuracy.To illustrate and test this methodology, we use seal pup photo counts from an aerial pho-tographic survey in the Greenland Sea in March 2012 with two different types of seals, harpand hooded seals. The data set contains the spatial location of each photo and the correspond-ing pup count. To be more informative about the non-observed areas, we include covariateinformation extracted from satellite imagery captured on the very same date as the aerialphotographic survey was conducted, to act as a proxy for ice thickness. This is importantas the seal pups can only be observed on ice, with non-observable pups accounted for withinthe dynamic population model (Øig˚ard et al., 2010). Compared to the other procedures ourmethod gives larger uncertainties, especially for the harp seals. To validate these differences,we compare our proposed method with a number of reference methods in two cross-validationexperiments, one where random sets of photos are removed and one where whole transects areremoved from the data set prior to inference. Performance assessment based on proper scoringrules suggests our method performs best on a local level, and comparable on a more globalscale. Further calibration assessment suggests that the larger uncertainty in our method isindeed more realistic.The rest of the paper is organized as follows: Section 2 describes the survey method usedto gather seal pup observational data, in addition to specific details related to our particularseal pup data set. The satellite imagery and the covariate information extracted therefrom arealso discussed. Relevant background related to point processes and aggregated point patternsare given in Section 3. Section 4 describes the details of our suggested modeling approach andthree references methods that we compare our method against. The validation schemes usedto verify and compare the different approaches are also described. The results are presentedin Section 5, including specifics of the model fitted with our procedure, and the validation andcomparison results. The final Section 6 contains concluding remarks and pointers to futurework. 2.
Data
In this section we describe the survey method and additional modeling information we haveobtained through satellite imagery.2.1.
Survey method.
Before conducting the aerial photographic survey with the purpose ofmonitoring the seal pup production, the marine researchers typically perform a helicopter re-connaissance survey. This is done in order to locate the patches where the seals whelp forlimiting the survey area for the more expensive airborne photographic survey. The actual pho-tographic survey is conducted by flying a survey aircraft equipped with advanced photographicequipment and GPS along a number of transects at a fixed distance that sparsely cover thesurvey area.In this particular survey in March 2012, the airplane flew at an altitude of about 330m, andtook a total of 2792 photos along 27 parallel transects, approximately 3Nm ( ≈ × STIMATING SEAL PUP PRODUCTION 3 two southernmost transects, which were flown at an altitude of 250m, with the photos covering170 × ≈ . Seal pup counts.
Following the airborne survey, experienced marine researchers manu-ally count the number of seal pups of each species in each photo . Quality checks with multipleexaminations are performed to limit the measurement error introduced in this step (Øig˚ardet al., 2014a,b). The seal pup count data set used on the subsequent analysis contains thecoordinates and extent of each photo, in addition to the number of seal pups of each speciesobserved. The data are plotted in Figure 1 along with the transect locations and the extentof the whelping region. As seen from the figure, there tends to be more seal pups clusteredtowards the middle eastern boundary and southern corner of the whelping region. Compara-tively more harp seal pups are observed than hooded seal pups and the spatial distribution ofthe former seems less homogeneous. −80−40040 −25 0 25 50 x (km) y ( k m ) Lines
TransectWhelping region
Seal type
HarpsHooded
Seal count
Figure 1.
Harp and hooded seal pup count data from the Greenland Sea 2012survey, with east-west direction on the x -axis and north-south direction on the y -axis.2.3. Satellite imagery.
For whelping, the seals require large ice floes with access to the oceanfor the adult seals to access food. Information regarding which areas are covered by ice floesand which areas are merely open water, is thus potentially highly relevant when estimatingthe seal pup production. In an attempt to account for this, we have collected high resolutionsatellite imagery (Modis) from the whelping region captured on the same day as the airbornephotographic survey was conducted. From this satellite imagery we have extracted a variablewhich acts as a proxy for the ice thickness. This density variable is displayed in Figure 2.Comparing the satellite data to the seal pup counts in Figure 1, we see that the seal pup The exact position of the seals are not recorded.
STIMATING SEAL PUP PRODUCTION 4 counts appear to be higher in the areas with high ice density, than in areas with lower icedensity. −100−50050 −30 0 30 x (km) y ( k m ) Density
Whelping region
Figure 2.
Satellite data showing ice density used as a covariate in the modelfitting process. 3.
Background
In the present paper, we assume that the seal pup appearances can be thought of as arealization of a point process. The aim of the analysis is to estimate the total number of sealpups, N (Ω), in the whelping region Ω ⊂ R (indicated in black in Figures 1 and 2).3.1. Point processes.
A spatial stochastic point process is a mathematical description of therandom process through which points or locations are distributed in space. The collection ofsuch observations is called a point pattern. In informal mathematical terms, a spatial pointprocess Y is a random collection of points in a bounded observation region Ω ⊂ R where boththe number of points and their locations are random. The most fundamental type of pointprocess is the Poisson point process . It may be specified by a deterministic intensity function λ : Ω (cid:55)→ [0 , ∞ ) which determines the density of points in any location in Ω. The number ofpoints, N ( B ), of any Borel set B ⊆ Ω is Poisson distributed with mean µ ( B ) = (cid:82) B λ ( s ) d s ,i.e. N ( B ) ∼ Po( µ ( B )). Further, N ( B ) is independent of N ( B ∗ ) for any other non-overlappingBorel set B ∗ ⊆ Ω with N ( B ∪ B ∗ ) = N ( B ) + N ( B ∗ ).The Cox process (also known as the doubly stochastic Poisson point process) introducedby Cox (1955) is a generalization of the Poisson point process where the intensity function λ is in itself stochastic. A popular special case of this hierarchical model is the log-GaussianCox process (LGCP) (Møller et al., 1998), where the intensity is assumed to be log-Gaussian;i.e. there is an underlying latent Gaussian random field Z , and given λ = exp( Z ), Y is a Poissonpoint process with intensity λ . This type of model is quite flexible and useful for modeling agreat variety of natural processes, in particular when only a single realization of the processis available (see e.g. Diggle, 2013; Illian et al., 2008; Møller & Waagepetersen, 2003). TheLGCP model can describe various types of clustering in the point pattern through the positivesemi-definite correlation function of the underlying Gaussian process (Møller et al., 1998).Wolpert & Ickstadt (1998) consider a Cox process where the stochastic intensity function isgiven by λ = Γ for a latent gamma random field Γ. As before, the conditional distribution of N ( B ) given λ is a Poisson distribution. For this model, the marginal distribution of N ( B ) is STIMATING SEAL PUP PRODUCTION 5 available in closed form and given by a negative binomial distribution (Mat´ern, 1971; Diggle &Milne, 1983).3.2.
Models for aggregated point patterns.
As mentioned in Section 2.2, the exact po-sitions of the counted seal pups are not available in the survey data. Instead, the seal pupcounts are provided as aggregated counts per photo. Thus, our data is a partly observed pointpattern aggregated to counts on an irregular lattice, as opposed to the actual point pattern.Methods that fit doubly stochastic Poisson process models to point pattern data, such as thoseproposed by Wolpert & Ickstadt (1998), Simpson et al. (2016) and Yuan et al. (2017), are thusnot directly applicable. Rather, our setting can be viewed as an extension of the approachconsidered in Rue et al. (2009) where a fully observed point pattern is approximated by thecorresponding counts on a regular lattice.To this end, let N ( A i ) for i = 1 , . . . , n denote the number of seal pups in each of the n = 2792photos with domains A i ⊂ Ω and (potentially varying) areas | A i | for i = 1 , . . . , n . Per propertiesof Poisson driven point processes, conditional on the intensity λ , being either deterministic(Poisson process) or random (Cox process), we have that N ( A i ) is Poisson distributed withparameter µ = (cid:82) A i λ ( s ) d s , i.e. N ( A i ) | λ ∼ Po (cid:16) µ = (cid:90) A i λ ( s ) d s (cid:17) , i = 1 , . . . , n. (1)Working under a Bayesian paradigm, we obtain an estimate of the total number of seal pupsin Ω given by the posterior predictive distribution of N (Ω), p (cid:0) N (Ω) | N ( A ) , . . . , N ( A n ) (cid:1) = (cid:90) p (cid:0) N (Ω) , λ | N ( A ) , . . . , N ( A n ) (cid:1) d λ = (cid:90) p (cid:0) N (Ω) | λ (cid:1) p (cid:0) λ | N ( A ) , . . . , N ( A n ) (cid:1) d λ, (2)where N (Ω) | λ ∼ Po( µ = (cid:82) Ω λ ( s ) d s ). The estimate for N (Ω) is thus given by a mixture ofPoisson distributions rather than a single Poisson distribution under both deterministic and therandom models for the intensity function λ . As the data are overdispersive, this may improvethe fit of the Poisson model.The negative binomial distribution is a common choice for modeling random counts in ap-plications with clustering. Due to the overdispersion in the data, Salberg et al. (2009) modelthe seal pup counts using a generalized additive model (GAM) based on a negative binomiallikelihood. Diggle & Milne (1983) show that only two point process models yield a negativebinomial marginal count distribution. One such model is the so-called compound Poisson pro-cess where each point of an underlying Poisson process is replaced by a random number ofcoincident points where the numbers are independent and identically distributed according toa logarithmic distribution. The second example is the Cox process with an intensity functiongiven by a gamma random field considered by Wolpert & Ickstadt (1998). Similarly, if a Poissonprocess has a constant intensity and the inference is performed using a conjugate gamma priordistribution for λ , the resulting posterior predictive distribution in (2) is a negative binomialdistribution, see Section 4.3 below. 4. Methods
We consider four different approaches to modeling the seal pup counts, see the summaryin Table 1. Our new proposed method is based on the LGCP framework and accounts forpotential clustering, or overdispersion, in the data through the underlying latent Gaussianrandom field. This approach is compared to the following reference approaches: A GAMapproach previously proposed by Salberg et al. (2009) and Øig˚ard et al. (2010) under both aPoisson and a negative binomial likelihood, and a homogeneous Poisson model. The methodsand the associated inference approaches are described below.
STIMATING SEAL PUP PRODUCTION 6
Table 1.
Summary of fitted models and inference approaches. Here, s is aspatial index, λ , α and β are coefficients, x is a vector of covariates, f denotesa latent spatial random effect, and S a deterministic spatial smoothing compo-nent.Name Model Intensity EstimationLGCP Poisson exp( α + β (cid:62) x ( s ) + f ( s )) SPDE-INLAGAM Po Poisson exp( α + β (cid:62) x ( s ) + S ( s )) GAMGAM NB Negative Binomial exp( α + β (cid:62) x ( s ) + S ( s )) GAMHom Po Poisson λ Bayesian4.1.
Stochastic intensity model.
Our proposed approach (LGCP) is based on a Poissonpoint process with a stochastic intensity function (i.e. a Cox process). The stochastic intensityfunction takes the form(3) λ ( s ) = exp( Z f ( s )) , where Z f is a Gaussian random field, or an LGCP with a conditional distribution for theaggregated counts on the form of (1). The continuous Gaussian random field takes the form Z f ( s ) = α + β (cid:62) x ( s ) + f ( s ) , (4)where α is an intercept term, β are regression coefficients for x ( s ) = ( q ( s ) , s , s , (cid:112) s + s ) (cid:62) ,with q ( s ) containing the ice density variable from the satellite imagery, while s , s , and (cid:112) s + s model linear spatial effects for s = ( s , s ) (cid:62) ∈ Ω. Finally, f ( s ) is a (non-linear)continuous Gaussian random field meant to model spatial dependence not captured by thecovariates in x . Specifically, f ( s ) is given the Mat`ern covariance function of the formCov( f ( s ) , f ( t )) = σ ν − Γ( ν ) ( κ (cid:107) s − t (cid:107) ) ν K ν ( κ (cid:107) s − t (cid:107) ) , (5)for s , t ∈ Ω where ν > K ν is the modified Bessel function ofthe second kind, κ > σ is the marginal variance. Further, foridentifiability of the intercept term in (4), we restrict f ( s ) to integrate to zero over the modelingregion.Note that the model defined by (1), (3) and (4) is equivalent to the Poisson log-normal modelof Christensen & Waagepetersen (2002) which is a special case of the spatial generalized linearmixed model framework proposed by Diggle et al. (1998). Christensen & Waagepetersen (2002)employ a Markov chain Monte Carlo (MCMC) algorithm for model inference using a slightlysimpler covariance structure than the Mat´ern covariance function in (5) to analyze a data setwhere the observation sets { A i } ni =1 are given by discs of fixed radius. Similar methods basedon Poisson kriging are discussed in e.g. Bellier et al. (2010) and De Oliveira (2014) where moretraditional geostatistical inference methods are employed.4.1.1. Inference.
To obtain an approximation to the posterior predictive distribution in (2), weapply the integrated nested Laplace approximation (INLA) of Rue et al. (2009) that allows forcomputationally feasible approximate Bayesian inference with discrete space latent Gaussianmodels, see Appendix A for a brief description. The stochastic partial differential equation(SPDE) approach of Lindgren et al. (2011) extends the INLA framework to also handle modelswith continuous latent fields as in (4). The SPDE approach is based on transforming thecontinuous latent field to a certain Gaussian Markov random field (GMRF), formulated throughthe solution of a SPDE. The key point is to approximate the continuous field Z ( s ) by a field Z GMRF ( s ) living on a triangular mesh. For a triangular mesh with m triangle vertices, we write Z GMRF ( s ) = m (cid:88) j =1 z j φ j ( s ) , (6) STIMATING SEAL PUP PRODUCTION 7 where z = ( z , . . . , z m ) (cid:62) is a multivariate Gaussian random vector and { φ j ( s ) } mj =1 is a setof deterministic linearly independent basis functions which are piecewise linear between thevertices and chosen such that φ j ( s ) is 1 at vertex j , and 0 at all other vertices. A consequenceof the representation in (6) is that Z GMRF ( s ) is fully determined by z . That is, Z GMRF ( s )takes the value z j at vertex j while its values inside the triangles are determined by linearinterpolation.Assume now that Z ( s ) is equipped with the Mat`ern covariance function in (5). As this typeof field is a solution to a certain SPDE, the precision matrix Q of z takes an analytical formwhich can be approximated by a sparse matrix ˜ Q . Since Z GMRF ( s ) is completely determinedby z , this allows continuous field computations to be carried out approximately using the INLAimplementation. Note that following certain guidelines for constructing the triangular mesh,the resulting approximation error is typically small (Lindgren et al., 2011; Simpson et al., 2012).For a complete introduction and review of the INLA framework, including the SPDE approach,see e.g. Blangiardo & Cameletti (2015) and Rue et al. (2016).4.1.2. Model specification and fitting.
The model parameters to be estimated are the regressionparameters ( α, β ) in (4) and the hyperparameters ( ν, κ, σ ) of the Mat´ern covariance function in(5). To improve identifiability, we fix the Mat`ern smoothing parameter at ν = 2, as is commonwhen applying the INLA framework (e.g. Blangiardo & Cameletti, 2015). The other hyperpa-rameters related to the latent field are equipped with mesh dependent default priors specifiedby the INLA software. For our mesh these are θ ∼ N (1 . ,
10) and θ ∼ N ( − . , θ = log( τ ) and θ = log( κ ), with σ = 1 / (4 πκ τ ). The intercept term α is assignedthe improper prior N (0 , ∞ ) while β ∼ N (0 , I ), where I denotes the 4 × p ( Z | N ( A ) , . . . , N ( A n )) , see Blangiardo & Cameletti (2015, Ch. 8.2). This allows us to use a Monte Carlo approximationfor the integral in (2): (cid:90) p ( N (Ω) | Z ) p ( Z | N ( A ) , . . . , N ( A n )) d Z ≈ K K (cid:88) k =1 p ( N (Ω) | ˜ Z k ) , STIMATING SEAL PUP PRODUCTION 8
Figure 3.
The triangular mesh used in the SPDE-INLA analysis for the LGCPmodel. The bottom right corner shows a zoomed in version of the mesh for theblue area above.where ˜ Z k is the k -th sample of the posterior latent field. We use K = 10 000 in this and similarMonte Carlo integrations throughout the paper. Further, by the point process properties, p ( N (Ω) | Z = ˜ Z k ) ∼ Po (cid:16) µ = (cid:90) Ω exp( ˜ Z k ( s )) d s (cid:17) . The integral (cid:82) Ω exp( ˜ Z k ( s )) d s can be solved by e.g. a simple Riemann midpoint rule. This isachieved by dividing Ω into J ≈
40 000 rectangles B , . . . , B J centered in s , . . . , s J , similarin size to the individual photos { A i } ni =1 , and using (cid:82) Ω exp( ˜ Z k ( s )) d s ≈ (cid:80) Jj =1 exp( ˜ Z k ( s j )) | B j | where the values ˜ Z k ( s j ) are derived using (6). To find the J sets, the full modeling regionis covered by a regular grid and the cells within the whelping region are selected. The finalapproximation to the posterior predictive distribution is thus p (cid:0) N (Ω) | N ( A ) , . . . , N ( A n ) (cid:1) ≈ K K (cid:88) k =1 Po (cid:16) µ = J (cid:88) j =1 exp( ˜ Z k ( s j )) | B j | (cid:17) , (7)i.e. a Poisson mixture distribution.4.2. Inhomogeneous intensity model.
An alternative to the stochastic intensity functionin (4) is a deterministic inhomogeneous intensity. With only one realization of the point pat-tern, it is generally not possible to distinguish between a doubly stochastic Poisson process andan inhomogenous Poisson process (Møller & Waagepetersen, 2003). However, the underlyingassumptions regarding the data structure (the LGCP e.g. assumes clustering) and the compo-nents of the intensity model vary somewhat so that the resulting predictive distributions maydiffer.The models we use here are motivated by the work of Salberg et al. (2009) who modelthe seal pup production by a generalized additive model (GAM) based on a negative binomiallikelihood, except that they do not include covariate information. Apart from adding covariates,the below approach follows that of Salberg et al. (2009) and Øig˚ard et al. (2010). Salberg et al.(2009) argue that a negative binomial likelihood should be used rather than a Poisson likelihooddue to overdispersion in that data. However, as the inhomogeneous intensity and the Poissonmixture in the posterior predictive distribution may improve the fit of the Poisson model, wefind it natural to include the inhomogenous Poisson model in our list of models.
STIMATING SEAL PUP PRODUCTION 9
Poisson model.
The inhomogeneous Poisson model takes the conditional Poisson formof (1) with intensity function λ ( s ) = exp( η ( s )) given by η ( s ) = α + β (cid:62) x ( s ) + S ( s ) , (8)where S ( · ) is a spatial smoothing component given by a thin-plate smoothing regression spline(Wood, 2003), i.e. a smooth, nonlinear deterministic spatial effect. As noted by Wood (2017,Ch. 5.8) there is a certain duality between the smoothing component S ( · ) in (8) and theGaussian random field f ( · ) in (4) in that S ( · ) can be interpreted as the mean of a latentrandom field. An assessment of uncertainty in the estimation of S ( · ) then corresponds to anassessment of the uncertainty in the mean of the random field.To fit the model in (8), we follow Salberg et al. (2009), and rely on the gam function in the R -package mgcv (Wood, 2017). This function fits the model N ( A i ) ∼ Po (cid:0) µ = | A i | exp( η ( s i )) (cid:1) , for i = 1 , . . . , n, with the | A i | values as fixed offsets with overlapping cubic regressions on a set of artificial knotsin space, and A i centered in s i . A generalized cross-validation (GCV) criterion is used to selectthe right amount of smoothing.The method essentially returns maximum likelihood (ML) estimates for the regression pa-rameters α, β and parameters ν associated with the spatial smoothing component S ( · ), θ =( α, β , ν ). To account for the uncertainty involved in the Poisson model we rely on the asymp-totic distribution of θ , N (cid:0) θ , Cov( (cid:98) θ ) (cid:1) , (9)where Cov( (cid:98) θ ) is the asymptotic covariance of (cid:98) θ estimated by gam . We sample from (9) andapproximate the integral over Ω by a simple Riemann midpoint rule using the J rectangles B , . . . , B J centered in s , . . . , s J with B ∪ . . . ∪ B J ≈ Ω, as in (7). This gives the approximation p GAM Po (cid:0) N (Ω) | N ( A ) , . . . , N ( A n ) (cid:1) ≈ K K (cid:88) k =1 Po (cid:16) µ = J (cid:88) j =1 | B j | exp( η ˜ θ k ( s j )) (cid:17) , (10)where η ˜ θ k ( s j ), is the value of (8) using θ -parameters corresponding to the k -th sample from(9).4.2.2. Negative binomial model.
An alternative to the inhomogeneous Poisson model is to con-sider the negative binomial model N ( A i ) | λ ∼ NegBin (cid:16) µ = (cid:90) A i λ ( s ) d s , τ (cid:17) , (11)with λ ( s ) = exp( η ( s )) where η ( s ) is given in (8), and shape parameter τ . This model is fit inthe same manner as the Poisson model above using N ( A i ) | λ ∼ NegBin (cid:0) µ = | A i | exp( η ( s i )) , τ (cid:1) , for i = 1 , . . . , n. In contrast to the Poisson distribution, the negative binomial distribution is not closed underaddition for different mean values. Thus, to arrive at a full posterior predictive distribution forthe total seal pup count N (Ω) under the model in (11), we employ a sampling procedure thatrelies on an underlying conditional independence condition stating that conditional on µ and τ ,the point counts in disjoint sets are independent. Again, using the J disjoint sets B , . . . , B J ,we sample counts N k ( B j ) | λ ∼ NegBin (cid:16) µ = | B j | exp( η ˜ θ k ( s j )) , τ = (cid:98) τ j (cid:17) , for j = 1 , . . . , J ; k = 1 , . . . , K, where s j is the center point of B j , η ˜ θ k ( s j ) is sampled from (9) as in (10), and (cid:98) τ j is the estimatedshape parameter for B j . The posterior predictive distribution for N (Ω) is then given by the STIMATING SEAL PUP PRODUCTION 10 empirical distribution function of { N k (Ω) } Kk =1 where N k (Ω) = J (cid:88) j =1 N k ( B j ) , for k = 1 , . . . , K. Both the sampling from the asymptotic normal distribution in (9) and the conditional inde-pendence assumption employed above, are also used by Salberg et al. (2009) to quantify theuncertainty around the total seal pup production.4.3.
Homogeneous intensity model.
A simple reference model is a homogeneous Poissonmodel with a constant intensity. That is, we set λ ( s ) ≡ λ for a fixed scalar λ and assume N ( A i ) | λ ∼ Po( µ = | A i | λ ) , for i = 1 , . . . , n. For a Bayesian inference, we equip λ with a non-informative conjugate gamma prior, λ ∼ Γ( a = 10 , b = 10) , for both seal types, where a and b are the shape and rate parameters, respectively. Thisresults in the following posterior distribution, λ | N ( A ) , . . . , N ( A n ) ∼ Γ (cid:16) a = a + n (cid:88) i =1 N ( A i ) , b = b + n (cid:88) i =1 | A i | (cid:17) . Moreover, the posterior predictive distribution in (2) takes the form of a negative binomialdistribution and is available in a closed form p Hom Po (cid:0) N (Ω) | N ( A ) , . . . , N ( A n ) (cid:1) = (cid:90) Po( | Ω | λ )Γ( a, b ) d λ = NegBin( µ = | Ω | a/b, τ = a ) . (12)where a and b are specified in the posterior distribution for λ above.The motivation for the homogeneous Poisson model is not only that it is a simple specialcase of our proposed LGCP model, but also that it is similar to the traditional approach toestimate seal pup production based on aerial photographic transect surveys, often referred toas Kingsley’s method (Kingsley et al., 1985). Kingsley’s method is fundamentally simple: Foreach transect T , . . . , T covering the space A T k = (cid:83) { A i ⊂ T k } A i , compute the seal pup count N T k = (cid:80) A i ∈ T k N ( A i ). Then, the estimate of the total number of seal pups is (cid:98) N (Ω) = | A Ω | (cid:80) k =1 | A T k | (cid:88) k =1 N T k . (13)Kingsley et al. (1985) also provide an estimate of the variance related to the seal pup productionestimate, based on serial differences between the transects. Salberg et al. (2008) later provideda modification to this variance estimate, which we have used here. Since this method onlyprovides a point and a variance estimate, it is difficult to properly compare it against theremaining Bayesian procedures. We will therefore not perform validation tests, as described inSection 4.4, for this method.4.4. Verification.
We compare the various modeling approaches using a cross-validation schemewhere we rely on two procedures for subsetting the data. The first procedure is a standard10-fold cross-validation setup, where we randomly remove 10% of the photos each time, suchthat each photo is removed exactly once. In the second procedure we remove all photos inone full transect at a time, such that each transect is removed exactly once, leaving us with27 different subsets. For both procedures, we fit the competing models for every subset andcompute posterior predictive distributions for every photo that is removed along with posteriorpredictive distributions for the sum of the removed photos (corresponding to the full transectfor the latter procedure).We compare the predictive performance of the various modeling approaches using two per-formance measures: The logarithmic score (Good, 1952) and the continuous ranked probability
STIMATING SEAL PUP PRODUCTION 11 . . . . . Value P r obab ili t y Value P r obab ili t y S c o r e Figure 4.
Illustration of the logScore (solid lines) and CRPS (dashed lines)functions for different values of y true , with two different predictive distributions(indicated in gray), Po( µ = 3) (left plot) and Po( µ = 20) (right plot).score (CRPS) (Matheson & Winkler, 1976), which both are proper scoring rules that assessfull predictive distributions (Gneiting & Raftery, 2007). Denoting a generic posterior predic-tive distribution by g ( x ), the corresponding cumulative distribution function by G ( x ), and theobserved count by y true , the two performance measures takes the formlogScore( g, y true ) = − log( g ( y true )) , (14) CRPS( G, y true ) = (cid:90) ∞−∞ ( G ( x ) − { x>y true } ( x )) d x, where {·} ( x ) denotes the indicator function. For both measures, smaller values reflect a bettermodel.Both the logarithmic score and the CRPS are optimized in expectation when the true datadistribution is issued as the forecast (Gneiting & Raftery, 2007). However, in real-life situations,average scores are often associated with high uncertainty (Thorarinsdottir & Schuhen, 2018).We thus apply two scores that penalize prediction errors in slightly different manners, seeFigure 4. In particular, the logarithmic score is more sensitive to outliers, cf. the left plotin Figure 4. For every cross-validation schemes, we present mean scores, where the mean iscomputed over the different folds. To indicate the variation in the scores, we compute 90%bootstrapped confidence intervals (CI) of those means by repeated re-sampling (10 000 times,with replacement) of the fold scores.We further assess the calibration, or the prediction uncertainty, of the methods by assessingthe coverage of the posterior predictive distributions. That is, we check how often y true lieswithin different credibility intervals, compared to their intended coverage – small uncertainty isof no value if it is not reflecting the true uncertainty of the model. For a calibrated forecast, anevent predicted with probability p should be realized with the same frequency in the observeddata. In order to mimic quantification of the prediction uncertainty for the complete whelpingregion as closely as possible, we perform this exercise on the transect level.5. Results
Here, we present the results obtained when applying the various models/estimation methodsin Table 1 to per-photo count data from the 2012 survey of the Greenland sea whelping region.We model hooded and harp seals separately as their occurrences are expected to be indepen-dent conditional on the covariate information. Specifically, we compare our LGCP approachestimated using SPDE-INLA with the GAM-based procedure, both with a negative-binomialdistribution for the counts (GAM NB) and with the simpler Poisson distribution (GAM Po).
STIMATING SEAL PUP PRODUCTION 12 −80−40040 −25 0 25 50 x (km) y ( k m ) −2.50.02.5 mean −80−40040 −25 0 25 50 x (km) y ( k m ) sd Figure 5.
The mean and standard deviation of the random field Z f ( s ) in (4),fitted for hooded seals with our LGCP approach. Table 2.
Summary table for the posterior predictive distributions of the totalcount of hooded seal pups in the whelping region.mean median mode IQR 0.025-quantile 0.975-quantileLGCP 11649 11503 11472 1699 9472 14741GAM NB 11178 11157 11093 807 10075 12395GAM Po 11296 11292 11093 572 10467 12147Hom Po 11166 11161 11172 555 10372 11987As a baseline model we use a homogeneous Poisson model with no covariates, spatial term, orother random effects (Hom Po).5.1.
Hooded seals.
Within the flight transect sparsely covering the whelping region, a totalof 777 hooded seal pups were counted. The blue dots in Figure 1 show how these are spreadon the 2792 photos, with between 0 and 12 pups per photo.Figure 5 shows the mean and standard deviation of the random field fitted using our LGCPprocedure outlined in Section 4.1. As seen, the latent field captures the high intensity of hoodedseal pups in the middle-eastern part of the whelping region. This area has a medium range icethickness, cf. Figure 2. There is also an increased seal pup intensity further south, in particularcloser to the open water. Generally speaking, as seen from the standard deviation plot, theuncertainty is rather large where the intensity is low, while it is smaller where the intensityis high. This means that apart from the north, one is fairly certain that there are some sealpups in areas where seal pups are observed nearby, while there could very well exist seal pupsin locations where none are observed at the neighboring observation sites. As a result of veryfew seals being observed in the north, the mean intensity here is so low that it is unlikely thata significant amount of seals have settled there.The range of the latent field, defined as the distance at which the spatial correlation isapproximately 0.1, has a posterior mean of 3.63 km, or about 2/3 of the distance between twotransects. This means that in an area lying between two transects, the latent field is essentiallydetermined by the two neighboring transects. Further, the fitted model gives the followingposterior means for the intercept ( α ) and the fixed effects ( β ): mean α = − . , mean β ,q =9 . , mean β ,s ,s ,s = (0 . , − . , − . (cid:98) α = − . , (cid:98) β q = 9 .
59, and (cid:98) β s ,s ,s = (0 . , − . , . STIMATING SEAL PUP PRODUCTION 13 l Hooded pup count P r obab ili t y Model
LGCPGAM NBGAM PoHom Po l Kingsley
Figure 6.
The posterior predictive distributions for the total hooded pupcounts in the whelping region for the five competing models. For Kingsley’smethod, we show the point estimate +/- 2 times the estimated standard devia-tion, corresponding to an approximate 95% confidence interval under a normaldistribution assumption.Figure 6 shows the posterior predictive distribution for the total pup count in the whelpingregion using our LGCP procedure, along with the corresponding results for the two GAM-basedprocedures (negative binomial and Poisson response) and the homogeneous Poisson model. Asimple summary of Kingsley’s method is also given for reference. Table 2 summarizes thepredictive distributions. While most of the mass coincides for the four different approaches,the GAM NB method yields the lowest estimate and our LGCP method the highest estimate.The Hom Po approach and the GAM Po methods have the lowest prediction uncertainty. Here,our LGCP method has an interquartile range approximately twice as large as the GAM NBapproach, and three times larger than GAM Po and Hom Po approaches. For all the methods,the mean predictor is slightly higher than the median predictor, indicating a small degree ofskewness with a heavier upper tail. This effect is most pronounced for the LGCP approach.Table 3 shows the results from the validation scheme applied to the four methods we comparehere, as outlined in Section 4.4. At the photo level, we issue a prediction for the pup countper photo, for either 10% of the photos or all photos in a single transect at a time. Here, ourLGCP method yields very good results, in particular for the random 10-fold cross-validationwhere observations are generally available in the neighborhood of the prediction locations.It is significantly the best method for this setting under the CRPS, defined as having non-overlapping 90% CI strictly below all others, and almost significant in terms of the logScore.When leaving out a full transect at a time, our LGCP method and GAM NB method performvery similar, and somewhat better than the others.At the aggregate/transect level we issue a joint prediction for the total pup count per 10%of the photos or per transect, respectively. Here, the baseline homogeneous Poisson modelperforms very well with the random leave-out scheme, indicating that the data may be close tohomogeneous across the photos. However, this does not transfer to the second set-up where weleave out one transect at a time in which case our LGCP method and the GAM NB methodagain perform well. As seen in Figure 1, the data do not appear homogeneous across thetransects. Note that these reported scores are averages over relatively few predictions, only 10distinct predictions for the random 10-fold cross-validation study and 27 distinct predictionsfor the leave-out-transect setup. The score values are thus associated with a large degree of
STIMATING SEAL PUP PRODUCTION 14
Table 3.
Validation results on photo level (one prediction per photo) and tran-sect level (one prediction per transect), respectively. Lower and upper bounds of90% bootstrapped confidence intervals for the scores are shown in parenthesis.Cells shown in italics are the best (smallest) per column. Those which are sig-nificantly smaller than the others (defined as having non-overlapping confidenceintervals) are also in bold.
HOODED SEALS: PHOTO LEVEL
Random 10-fold CV Leave-out full transectCRPS logScore CRPS logScoreLGCP
GAM Po 0.22 (0.20, 0.24) 0.54 (0.51, 0.58) 0.24 (0.22, 0.26) 0.58 (0.54, 0.62)Hom Po 0.26 (0.24, 0.28) 0.77 (0.71, 0.84) 0.26 (0.24, 0.29) 0.78 (0.72, 0.85)
HOODED SEALS: AGGREGATE/TRANSECT LEVEL
Random 10-fold CV Leave-out full transectCRPS logScore CRPS logScoreLGCP 5.43 (4.04, 6.99) 3.68 (3.51, 3.86) 9.91 (5.99, 14.80)
GAM NB 5.93 (4.95, 7.00) 3.79 (3.68, 3.91)
Transect number P r ed i c t i on Model
LGCPGAM NBTruth
Figure 7.
The posterior predictive distributions per transect for the LGCP andGAM NB methods plotted against the true count for hooded seals. For eachmethod, the solid line shows the median, the light colored box shows the 50% CI,while the transparent box shows the 90% CI. The y-axis has a log ( x + 2)-scaleto better show differences.uncertainty. The uncertainty is particularly high for the homogeneous Poisson model that, dueto the assumption of homogeneity, is especially sensitive to inhomogeneities across the differentcross-validation sets.Based on the results in Table 3, our LGCP method performs somewhat similar to the GAMNB method and these are both clearly superior to the other two alternatives; the two methodscan only be distinguished in terms of CRPS on the photo level. Despite this, their resulting STIMATING SEAL PUP PRODUCTION 15 −80−40040 −25 0 25 50 x (km) y ( k m ) −4048 mean −80−40040 −25 0 25 50 x (km) y ( k m ) sd Figure 8.
The mean and standard deviation of the random field Z f ( s ) in (1)fitted for harp seals with our LGCP approach.posterior predictive distributions are quite different. To get further insight into this phenom-enon, Figure 7 shows the posterior predictive distributions for the two methods per transectin the leave-out-transect setup, plotted against the true transect counts. As seen from thefigure, the GAM NB method seems to have too narrow credibility intervals, while the LGCPapproach appears more calibrated. In fact, the 90% CI covers the true count in 26 / ≈ / ≈
67% of the transects for the GAMNB approach. The 50% CI is covered in respectively 16 / ≈
59% and 11 / ≈
41% ofthe transects for the LGCP and GAM NB approaches. Note that when looking at the MAE(mean absolute error) of the median count estimate for the two methods, the GAM NB methodachieves MAE of 12.4, while the LGCP method is slightly worse with MAE of 13.8. Thus, theGAM NB method seems to do somewhat better as a point estimator, while the LGCP is betterat estimating the uncertainty.5.2.
Harp seals.
A total of 6034 harp seal pups where observed on the 2792 photos from theaerial photographic survey. As illustrated by the red dots in Figure 1, there are much largerpacks of harp seal pups than hooded seal pups, indicating a higher degree of inhomogeneity.Here, the pup count per photo ranges from 0 to 160.Figure 8 shows the mean and standard deviation of the fitted random field Z f ( s ) using ourLGCP procedure. Compared to the latent field for the hooded seal pups in Figure 5, the meanfield here has a much higher degree of spatial variation with higher and steeper peaks. However,the locations where the seal pups mainly appear are similar to those for the hooded seal pups,except for some additional colonies in the north and north-west of the region. Otherwise, theproperties of the two fields are fairly similar. The range of the latent field has a posterior meanof 2.89 km, a slightly smaller value than for the hooded seals which corresponds to roughlyhalf the distance between two transects. The fitted model gives the following posterior meansfor intercept ( α ) and fixed effects ( β ): mean α = − . , mean β ,q = 14 . , mean β ,s ,s ,s =(0 . , . , − . (cid:98) α = 2 . , (cid:98) β q = 19 .
05, and (cid:98) β s ,s ,s =(0 . , . , − . STIMATING SEAL PUP PRODUCTION 16 l Harp pup count P r obab ili t y Model
LGCPGAM NBGAM PoHom Po l Kingsley
Figure 9.
The posterior predictive distributions for the total harp pup countsin the whelping region for the five different approaches. The x-axis is plotted onlog-scale. For Kingsley’s method, we show the point estimate +/- 2 times theestimated standard deviation, corresponding to an approximate 95% confidenceinterval under a normal distribution assumption.
Table 4.
Summary table for the posterior predictive distributions for the totalnumber of harp seal pups in the whelping region.mean median mode IQR 0.025-quantile 0.975-quantileLGCP 147919 127965 110996 72347 69267 357185GAM NB 98617 98035 91876 12895 81023 119349GAM Po 84852 84852 84910 1536 82681 87094Hom Po 85751 85747 85719 1539 83529 88004Po). A simple summary of Kingsley’s method is also given for reference. Table 4 summarizesthese distributions. The LGCP prediction of the total harp seal pup count is highly uncertain,essentially saying that the total number of seal pups could very well be above 250 000, butalso less than 100 000. In contrast, the GAM NB method’s upper tail ends at about 130 000seal pups, while the GAM Po and Hom Po methods agree that there are between 80 000 and90 000 harp seal pups within the whelping region. Thus, there are significant differences bothbetween the centrality and width of the different methods’ posterior predictive distributions.The predictive distributions for GAM Po and Hom Po are essentially symmetric while the onefor GAM NB is minimally skewed with a heavier upper tail. The LGCP approach, however,yields a severely skewed predictive distribution with a predictive mean that is 15% larger thanthe predictive median due to the large uncertainty in the estimation of the random field Z f ,cf. Figure 8.To further investigate the differences between the four predictions, we have plotted theposterior median intensity fields for the four methods in Figure 10, cf. the expressions for theintensity fields in the third column in Table 1. The regression coefficient for the ice densitycovariate q shown in Figure 2 is somewhat higher for the GAM approaches than for the LGCPapproach (posterior mean of 19.05 for GAM NB compared to 14.70 for LGCP), resulting infairly smooth median intensity fields that reflect the spatial structure of the ice density field.The median LGCP intensity field, however, has stronger inhomogeneities across space andappears more strongly influenced by the data density shown in Figure 1 with noticeable peaksin locations with higher observed seal pup density. Overall, however, the median intensity STIMATING SEAL PUP PRODUCTION 17
GAM Po Hom PoLGCP GAM NB −25 0 25 50 −25 0 25 50−80−40040−80−40040 x (km) y ( k m ) −4048 Median intensity(log scale)
Figure 10.
The estimated posterior median intensity field (plotted on a loga-rithmic scale) for all four methods fitted for harp seals.is higher for the GAM methods than for the LGCP. The reason that LGCP still produceshigher predictions (cf. Table 4) is partly due to the high peaks, but mainly due to the highdegree of uncertainty. The posterior of Z ( s ) is more or less symmetric with a high degree ofuncertainty, such that large values are sampled quite frequently. These large values boosts theintensity λ ( s ) = exp( Z ( s )) considerably, resulting in large predictions. The lower uncertaintyof the GAM methods does not have the same effect. Further, note that although the two GAMmodels have very similar posterior median intensities, their posterior predictive distributionsare quite different, cf. Figure 9 and Table 4.As for the hooded seals, we take a closer look at the LGCP and GAM NB methods to betterunderstand how well their different estimates of the prediction uncertainty match the actualuncertainty, see Figure 11. As expected, the GAM NB method has a much narrower credibilityintervals which too often fail to cover the true seal pup count in the transect, while our LGCPmethod shows much better calibration. Out of the 27 transects, the 90% credibility intervalsfor the LGCP and GAM NB approaches covers the true count in respectively 24 / ≈ / ≈
67% of the transects. The corresponding coverage for the 50% interval are14 / ≈
52% and 11 / ≈ STIMATING SEAL PUP PRODUCTION 18
Transect number P r ed i c t i on Model
LGCPGAM NBTruth
Figure 11.
The posterior predictive distributions per transect for the LGCPand GAM NB methods plotted against the true count for harp seal pups. Foreach method, the solid line shows the median, the light colored box showsthe 50% CI, while the transparent box shows the 90% CI. The y-axis takesa log ( x + 2)-scale to better show differences.hooded seals, the posterior median of the GAM NB method does better as a point estimatorin terms of MAE (130.2) than our LGCP method (179.1).Table 5 shows the results from the validation procedure for the harp seals, yielding similarmodel rankings as for the hooded seals. On photo level, the LGCP approach gives a significantlybetter CRPS and logScore under random 10-fold cross-validation. Leaving out full transectsgives no significantly best method although the GAM Po method tends to generally do well here.However, as before, these average scores are associated with a very high degree of uncertaintyso that a ranking of the methods based on these results is not advisable.6. Conclusions and discussion
We have presented a point process based approach to estimate seal pup production basedon observational data from an aerial photographic survey. Using the SPDE-INLA framework,we fit a Bayesian hierarchical model with Poisson counts following a log-Gaussian Cox process(LGCP) model formulation. As an additional contribution to seal pup production estimation,we adopt the use of satellite imagery as covariates in the modeling process, to act as a proxyfor ice thickness. The approach is applied to 2012 survey data from the Greenland Sea, withboth harp and hooded seal pup counts, and compared to several reference methods that canbe associated with non-homogeneous or homogeneous point process formulations rather thanthe doubly stochastic setting of the Cox process.The competing methods are compared in two cross-validation studies. The proposed LGCPapproach generally performs best locally, while no method stands out as the best on a moreregional scale. However, this lack of discrimination in the comparison at the regional scale isnot surprising given the relatively small size of our data set resulting in large uncertainties inthe scores, see e.g. the discussion and examples in Thorarinsdottir & Schuhen (2018). Themost distinguishing characters of the LGCP method are higher count predictions and a largeprediction uncertainty compared to the other methods. Our analysis suggests that the wideuncertainty bounds are indeed necessary to issue calibrated predictions. This further suggests
STIMATING SEAL PUP PRODUCTION 19
Table 5.
Validation results on photo level (one prediction per photo) and tran-sect level (one prediction per transect), respectively. Lower and upper bounds of90% bootstrapped confidence intervals for the scores are shown in parenthesis.Cells shown in italics are the best (smallest) per column. Those which are sig-nificantly smaller than the others (defined as having non-overlapping confidenceintervals) are also in bold.
HARP SEALS: PHOTO LEVEL
Random 10-fold CV Leave-out full transectCRPS logScore CRPS logScoreLGCP
GAM Po 2.32 (2.10, 2.55) 2.09 (2.00, 2.17) 2.46 (2.22, 2.71) 2.17 (2.08, 2.26)Hom Po 2.64 (2.40, 2.90) 3.47 (3.28, 3.67) 2.66 (2.42, 2.92) 3.49 (3.30, 3.69)
HARP SEALS: AGGREGATE/TRANSECT LEVEL
Random 10-fold CV Leave-out full transectCRPS logScore CRPS logScoreLGCP 95.98 (51.20, 148.78) 6.57 (5.93, 7.33) 152.95 (111.93, 198.59)
GAM NB 88.33 (70.94, 106.87)
STIMATING SEAL PUP PRODUCTION 20
For the LGCP inference, we have applied the SPDE-INLA approach directly as implementedin the R package INLA (Lindgren et al., 2011; Rue et al., 2009). Alternatively, the more recentpackage inlabru (Bachl et al., 2019) is based on the SPDE-INLA software to model pointpattern data from surveys with varying detection probabilities over the sampled area. Oursetting is slightly different, in that the detection probability is assumed to be equal to oneover the entire sampled area and we only have sampled counts per each photo rather than theprecise locations of the seal pups. While this somewhat simpler setting can also be analyzedwith inlabru , we have chosen to implement out own version which allows for a slightly higherflexibility in the generation of the mesh. Besides small differences in the manner in which edgeeffects are treated when approximating the integral over the latent field, our implementationshould give comparable results to inlabru .In our model specification, we only consider linear effects of the satellite covariate and linearspatial effects. We investigated including non-linear effects and squared terms both for theLGCP and the deterministic intensity model formulations. However, as this did not improvethe performance of the models, such terms were not included in the final model specifications.In the present work, harp and hooded seals have been modeled separately, based on a singlesurvey. As an alternative, one may consider building a joint model for harp and hoodedseals, for instance by using a common spatial field, in addition to seal specific ones (seee.g. Waagepetersen et al. (2016)). Further, due to drifting ice and moving seals, the spatiallocations of the seal pups cannot be directly compared from one survey to the next. However,most of the seals tend to stay in more or less the same packs from one year to another. Suchinformation could potentially be utilized to construct informative priors, which may reduce themodeling uncertainty. It would be interesting to see investigations on such attempts at borrowstrength, either from previous surveys or between seal types. Such investigations are, however,out of scope of the present paper.
Acknowledgments
Martin Jullum and Thordis Thorarinsdottir acknowledge the support of The Research Coun-cil of Norway through grant 240838 “Model selection and model verification for point processes”.Fabian Bachl was supported by EPSRC grant EP/K041053/1 “Modelling spatial distributionand change from wildlife survey data”. We thank Tor Arne Øig˚ard for helpful discussions andfor assisting with the photographic survey data which were provided by the Norwegian Instituteof Marine Research, Øivind Due Trier for retrieving the satellite image data, Dawid LawrenceMiller and Alex Lenkoski for helpful discussions regarding the various inference approaches.
Appendix A. The integrated nested Laplace approximation
The integrated nested Laplace approximation (INLA) methodology proposed by Rue et al.(2009), and implemented in the R -package INLA ( ), allows for computationallyfeasible approximate Bayesian inference for latent Gaussian models. In latent Gaussian models, n univariate observations y = ( y , . . . , y n ) (cid:62) are assumed to be conditionally independent given m latent Gaussian variables z = ( z , . . . , z m ) (cid:62) and a set of hyperparameters θ . More precisely,the INLA implementation covers models of the form p ( y | z , θ ) = n (cid:89) i =1 p ( y i | η i , θ ) , with η i = m (cid:88) j =1 c ij z j for fixed c ij ,p ( z | θ ) ∼ N( µ ( θ ) , Q ( θ ) − ) , (15) θ ∼ p ( θ ) , where the latent variables z may depend on additional (fixed) covariates, for a large class ofmodels for y .For computationally fast inference it is essential that the precision matrix Q ( θ ) is sparse andthat the parameter vector θ is of a fairly low dimension. This covers models where the latent STIMATING SEAL PUP PRODUCTION 21 field is a
Gaussian Markov random field (GMRF). For the inference, INLA utilizes severalnested Laplace approximations. That is, the posterior distribution of θ is approximated by p ( θ | y ) ≈ ˜ p ( θ | y ) ∝ p ( y , z , θ ) p G ( z | y , θ ) (cid:12)(cid:12)(cid:12)(cid:12) z = z ∗ ( θ ) , (16)where p G ( z | y , θ ) is a Gaussian approximation to the full conditional distribution of z , and z ∗ ( θ )is the mode of p ( z | y , θ ) for a given θ . The marginals of this low-dimensional posterior distribu-tion are typically computed by direct numerical integration. The marginals for the latent field, p ( z j | y ), are typically computed by first obtaining a Laplace approximation ˜ p ( z j | θ , y ) similarto (16), or a Taylor approximation of that distribution, and then solve (cid:82) ˜ p ( z j | θ , y )˜ p ( θ | y ) d θ bynumerical integration. See Rue et al. (2009) and Martins et al. (2013) for further details. ReferencesBachl, F. E. , Lindgren, F. , Borchers, D. L. & Illian, J. B. (2019). inlabru: an rpackage for bayesian spatial modelling from ecological survey data.
Methods in Ecology andEvolution , 760–766. Bellier, E. , Monestiez, P. & Guinet, C. (2010). Geostatistical modelling of wildlifepopulations: a non-stationary hierarchical model for count data. In
Geoenv Vii–geostatisticsfor environmental applications . Springer, pp. 1–12.
Blangiardo, M. & Cameletti, M. (2015).
Spatial and Spatio-temporal Bayesian Modelswith R - INLA . Wiley.
Christensen, O. F. & Waagepetersen, R. (2002). Bayesian prediction of spatial countdata using generalized linear mixed models.
Biometrics , 280–286. Cox, D. R. (1955). Some statistical methods connected with series of events.
Journal of theRoyal Statistical Society. Series B (Methodological) , 129–164.
De Oliveira, V. (2014). Poisson kriging: A closer investigation.
Spatial Statistics , 1–20. Diggle, P. J. (2013).
Statistical analysis of spatial and spatio-temporal point patterns . Chap-man and Hall/CRC.
Diggle, P. J. & Milne, R. K. (1983). Negative binomial quadrat counts and point processes.
Scandinavian journal of statistics , 257–267. Diggle, P. J. , Tawn, J. A. & Moyeed, R. (1998). Model-based geostatistics.
Journal ofthe Royal Statistical Society: Series C (Applied Statistics) , 299–350. Gneiting, T. & Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and esti-mation.
Journal of the American Statistical Association , 359–378.
Good, I. J. (1952). Rational decisions.
Journal of the Royal Statistical Society. Series B(Methodological) , 107–114.
Hastie, T. & Tibshirani, R. (1990).
Generalized additive models . Wiley Online Library.
ICES (2014). Report of the ICES/NAFO working group on harp and hooded seals (WGHARP).Tech. rep., ICES CM 2014/ACOM:20, Quebec City, Quebec, Canada.
Illian, J. , Penttinen, A. , Stoyan, H. & Stoyan, D. (2008).
Statistical analysis andmodelling of spatial point patterns , vol. 70. John Wiley & Sons.
Kingsley, M. , Stirling, I. & Calvert, W. (1985). The distribution and abundance of sealsin the Canadian High Arctic, 1980–82.
Canadian Journal of Fisheries and Aquatic Sciences , 1189–1210. Kovacs, K. M. & Lavigne, D. (1986). Cystophora cristata.
Mammalian Species , 1–9.
Krainski, E. T. , G´omez-Rubio, V. , Bakka, H. , Lenzi, A. , Castro-Camilo, D. , Simp-son, D. , Lindgren, F. & Rue, H. (2018).
Advanced Spatial Modeling with StochasticPartial Differential Equations Using R and INLA . Chapman and Hall/CRC.
Lindgren, F. , Rue, H. & Lindstr¨om, J. (2011). An explicit link between Gaussian fieldsand Gaussian Markov random fields: The stochastic partial differential equation approach.
Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 423–498. Martins, T. G. , Simpson, D. , Lindgren, F. & Rue, H. (2013). Bayesian computing withINLA: New features.
Computational Statistics & Data Analysis , 68–83. STIMATING SEAL PUP PRODUCTION 22
Mat´ern, B. (1971). Doubly stochastic poisson processes in the plane.
Statistical ecology ,195–213. Matheson, J. E. & Winkler, R. L. (1976). Scoring rules for continuous probability distri-butions.
Management science , 1087–1096. Møller, J. , Syversveen, A. R. & Waagepetersen, R. P. (1998). Log Gaussian Coxprocesses.
Scandinavian journal of statistics , 451–482. Møller, J. & Waagepetersen, R. P. (2003).
Statistical inference and simulation for spatialpoint processes . CRC Press.
Øig˚ard, T. A. , Haug, T. & Nilssen, K. T. (2014a). Current status of hooded seals inthe Greenland Sea. victims of climate change and predation?
Biological conservation ,29–36.
Øig˚ard, T. A. , Haug, T. & Nilssen, K. T. (2014b). From pup production to quotas:current status of harp seals in the Greenland Sea.
ICES Journal of Marine Science: Journaldu Conseil , 537–545. Øig˚ard, T. A. , Haug, T. , Nilssen, K. T. & Salberg, A.-B. (2010). Estimation of pupproduction of hooded and harp seals in the Greenland Sea in 2007: Reducing uncertaintyusing generalized additive models.
Journal of Northwest Atlantic Fishery Science , 103–123. Rue, H. , Martino, S. & Chopin, N. (2009). Approximate Bayesian inference for latentGaussian models using integrated nested Laplace approximations (with discussion).
Journalof the Royal Statistical Society, Series B , 319–392.
Rue, H. , Riebler, A. , Sørbye, S. H. , Illian, J. B. , Simpson, D. P. & Lindgren, F. K. (2016). Bayesian computing with INLA: A review. arXiv preprint arXiv:1604.00860 . Salberg, A.-B. , Haug, T. & Nilssen, K. T. (2008). Estimation of hooded seal (cystophoracristata) pup production in the Greenland Sea pack ice during the 2005 whelping season.
Polar Biology , 867. Salberg, A.-B. , Øig˚ard, T. A. , Stenson, G. B. , Haug, T. & Nilssen, K. T. (2009).Estimation of seal pup production from aerial surveys using generalized additive models.
Canadian Journal of Fisheries and Aquatic Sciences , 847–858. Sergeant, D. E. (1974). A rediscovered whelping population of hooded seals cystophoracristata erxleben and its possible relationship to other populations.
Polarforschung , 1–7. Sergeant, D. E. (1991). Harp seals, man and ice.
Canadian Journal of Fisheries and AquaticSciences , 1–153.
Simpson, D. , Illian, J. B. , Lindgren, F. , Sørbye, S. H. & Rue, H. (2016). Going off grid:Computationally efficient inference for log-gaussian cox processes.
Biometrika , 49–70.
Simpson, D. , Lindgren, F. & Rue, H. (2012). In order to make spatial statistics com-putationally feasible, we need to forget about the covariance function.
Environmetrics ,65–74. Simpson, D. , Rue, H. , Riebler, A. , Martins, T. G. , Sørbye, S. H. et al. (2017). Penal-ising model component complexity: A principled, practical approach to constructing priors.
Statistical Science , 1–28. Thorarinsdottir, T. L. & Schuhen, N. (2018). Verification: assessment of calibration andaccuracy. In
Statistical Postprocessing of Ensemble Forecasts . Elsevier, pp. 155–186.
Waagepetersen, R. , Guan, Y. , Jalilian, A. & Mateu, J. (2016). Analysis of multi-species point patterns by using multivariate log-gaussian cox processes.
Journal of the RoyalStatistical Society: Series C (Applied Statistics) , 77–96. Watson, D. F. (1981). Computing the n-dimensional delaunay tessellation with applicationto voronoi polytopes.
The computer journal , 167–172. Wolpert, R. L. & Ickstadt, K. (1998). Poisson/Gamma random field models for spatialstatistics.
Biometrika , 251–267. Wood, S. (2017).
Generalized Additive Models: An Introduction with R . Chapman andHall/CRC, 2nd ed.
Wood, S. N. (2003). Thin-plate regression splines.
Journal of the Royal Statistical Society(B) , 95–114. STIMATING SEAL PUP PRODUCTION 23
Wood, S. N. (2016). Just another Gibbs additive modeler: Interfacing JAGS and mgcv.
Journal of Statistical Software , 1–15. Yuan, Y. , Bachl, F. E. , Lindgren, F. , Borchers, D. L. , Illian, J. B. , Buckland,S. T. , Rue, H. , Gerrodette, T. et al. (2017). Point process models for spatio-temporaldistance sampling data from a large-scale survey of blue whales.
The Annals of AppliedStatistics11