Estimating Abundance from Counts in Large Data Sets of Irregularly Spaced Plots using Spatial Basis Functions

Monitoring plant and animal populations is an important goal for both academic research and management of natural resources. Successful management of populations often depends on obtaining estimates of their mean or total over a region. The basic problem considered in this paper is the estimation of a total from a sample of plots containing count data, but the plot placements are spatially irregular and non-randomized. Our application had counts from thousands of irregularly spaced aerial photo images. We used change-of-support methods to model counts in images as a realization of an inhomogeneous Poisson process that used spatial basis functions to model the spatial intensity surface. The method was very fast and took only a few seconds for thousands of images. The fitted intensity surface was integrated to provide an estimate from all unsampled areas, which is added to the observed counts. The proposed method also provides a finite area correction factor to variance estimation. The intensity surface from an inhomogeneous Poisson process tends to be too smooth for locally clustered points, typical of animal distributions, so we introduce several new overdispersion estimators due to poor performance of the classic one. We used simulated data to examine estimation bias and to investigate several variance estimators with overdispersion. A real example is given of harbor seal counts from aerial surveys in an Alaskan glacial fjord.


INTRODUCTION
Monitoring plant and animal populations is an important goal for both academic research and management of natural resources. Successful management of populations often depends on estimates of their mean or total over a region. Historically, this has been the purview of sampling theory using simple random sampling, stratified random sampling, etc., which are design-based methods. For design-based methods, sample units are chosen at random, measurements are made or observed from the sample units, and inference is derived from the inclusion probability for sample units (i.e., Horwitz- Thompson estimation). For overviews, see Cochran (1977) or Thompson (1992). An alternative approach developed in the early 1960's called geostatistics includes methods such as block kriging (Gandin 1959(Gandin , 1960Matheron 1963), which also estimates a regional total. These methods rely on an assumption about a stochastic process that generated the realized observations, and are hence "modelbased." Model-based inference relies on estimating parameters for the assumed model, and then forming probability statements (confidence intervals, prediction intervals, etc.) from the fitted model. In this paper we pursue the model-based approach because samples cannot always be drawn randomly. In particular, we consider counts from aerial photographs, which are difficult and inefficient to randomize, with the basic problem being the estimation of the total count for a region.
The goals and context of this paper are shown in Fig. 1. The situation for block kriging is shown in Fig. 1a. Let the spatial region of interest be R. Any particular location in R is given by x-and y-coordinates contained in the vector s = [s x , s y ] , and the random variable at the ith location is denoted Y (s i ). We assume a spatial random field {Y (s) : s ∈ R} (Cressie 1993, p. 30). The set {Y (s)} is continuous in space and hence infinite. We use the notation that for vector x = [x 1 , x 2 ] , x ≡ x 2 1 + x 2 2 , so s i − s j is the Euclidean distance between s i and s j . If the correlation between Y (s i ) and Y (s j ) goes to one as s i − s j → 0, then {Y (s i )} will form a smooth (differentiable) but random surface. Suppose that n observed values from the random surface are contained in the vector y = [y(s 1 ), . . . , y(s n )] (the solid circles in Fig. 1a). Block kriging uses a linear combination λ y to predict the average or total in block A; Y (A) ≡ A Y (u)du, where this integral is assumed to exist (Yaglom 1962, p. 23;Cressie 1993, p. 106). The salient feature of block kriging is that a model of the autocorrelation of the spatial random field can be estimated from the point-level data of the observations, and it is relatively easy to aggregate (through the integral) for an estimator of block A. However, extensions to count data have been difficult because data are modeled on a transformed space but the integral is desired on the original space (e.g., see Cressie 1993, p. 286). For example, Christensen and Waagepetersen (2002), Wikle (2002), Monestiez et al. (2006) develop maps from count data but do not attempt abundance estimates.
Counts are often obtained from plots, B i in Fig. 1b, that have substantial area, are in a regular grid, and exhaustively fill both R and A. Here, classical random sampling using design-based inferences are often employed, where a random sample of the observed count on the ith block, y(B i ), are used to estimate a total. To correctly estimate variance, a finite population correction factor is employed,v ar(Ŷ (A)) = (σ 2 /n)(1 − f ) whereσ 2 is the sample variance, f = n/N is the fraction of sampled units (n) to total sample units (N ) within A, and 1 − f is called the finite population correction factor. A model-based, finite population version of block kriging for the situation shown in Fig. 1b was developed by Ver Hoef (2000,2008). In the strict sense, distance between samples is not well-defined because of the non-point nature of each sample B i . However, models of autocorrelation are often built in this case by using the centroid of each plot. The main problem considered by Ver Hoef (2000,2008) was that the traditional formulation of block kriging as shown in  Hence, if one were to estimate an autocorrelation model for Fig. 1b, and then apply standard block kriging formulas (as has been done in the literature), there is no finite correction factor. For example, if we sampled all of the plots, then the prediction variance should be zero, but the traditional formulation of block kriging would have nonzero prediction variance. Hence, a finite population version was developed by Ver Hoef (2000,2008), and it performed well and had proper confidence intervals in a variety of situations (Ver Hoef 2002). Now consider the situation in Fig. 1c. Here, we would like to use counts from samples with substantial area, y(B i ), to predict at Y (A), but the sample units are not arranged in a regular grid that fills either R or A. Classical random sampling would usually be employed in this situation, with an estimator of the total being |A|Ȳ (B), whereȲ (B) is the total count divided by the sampled area, and the variance estimated by var(Ŷ (A)) = (σ 2 /n)(1 − f ), where here f = |a|/|A|, with |a| being the total area sampled. However, what if the samples Y (B i ) are not randomly placed, and may in fact be in some regular pattern that does not form a regular grid that does not exhaustively fill both D and A? The basic problem considered in this paper is the prediction of a total in region A from a sample of {Y (B i )} as in Fig. 1c, where sampling at random is not possible. This is the case for counts from photographs taken from aircraft. The National Marine Mammal Laboratory of the NOAA-NMFS Alaska Fisheries Science Center developed aerial survey methods to estimate and monitor harbor seal populations in glacial fjords in Alaska. We will give more details next as a motivating example.

MOTIVATING EXAMPLE
To make the problem concrete, we consider an example that prompted the model development. Aerial surveys were flown over the ice haul-out area of harbor seals in Icy Bay, Alaska. Ice emanating from tidewater glaciers provides a dynamic expanse of floating ice on which the seals whelp and nurse pups and rest during the molting season. The aerial platform, a twin-engine Aero Commander Shrike, was flown at 1,000 ft and ca. 100 knots on transects with variable spacing that were oriented in two main directions to sample the two main arms of the bay (Fig. 2), covering about 79 sq km. A vertically mounted camera (Nikon D1X with a 60 mm lens) captured an image approximately every 2 s through a portal, each covering about 80 × 120 m at the surface of the water. This firing rate, and the spacing of the transects, allowed for a gap between images of about 30 m end-to-end, and transects varied in spacing from side-to-side, but largely ensured that images were separated from each other; i.e., seals were sampled only once. The camera was usually turned off when flying over large areas of open water where hauled-out seals would necessarily be absent. This survey was conducted on 11 August, 2008, in the afternoon (1,300-1,430 h) when seals typically haul-out in peak numbers. This is just one data set collected among dozens annually as part of a monitoring program for harbor seals. Images were georeferenced and embedded as a raster layer in an ArcGIS (ESRI 2009) project allowing individual seals to be spatially marked in a point layer by visually inspecting each image (n = 2080 images). Footprints showing the extent of each image were generated as polygons (Fig. 2) in a separate layer, and seal points were summed within and assigned to each centroid and exported for statistical analysis. The spatial extent of each image was assumed to be constant despite small random variation in altitude (max: ±30 m) during the survey. The total spatial extent of each day's survey effort, over which the intensity surface would be calculated, was delineated by creating a polygon that corresponded to: (1) the coastline of the bay (shorelines from Alaska Department of Natural Resources line shapefiles), (2) an estimate of the location of the face of the glaciers (by connecting points that marked the glacial terminus in each of the northernmost images from every transect, and (3) the extent of the ice field defined as the edge of the images where ice cover (by area) dropped to <5 % by visual estimation. Areas of open water (<5 % cover) were delineated by donutholes in the overall polygon where the spatial boundaries were defined by the outermost images in which ice cover increased to ≥5 %. In other words, to minimize problems with selection bias, any area that could not contain a seal was eliminated by creating the proper boundary.
The 2,080 images covered 25.3 % of the study area. Of the 2,080 images, 180 of them had nonzero counts, so about 91 % were zero. A total of 1,002 seals were observed in the 180 plots. A maximum of 44 seals were counted in a single photograph. The data are summarized spatially in Fig. 2.

PREVIOUS WORK
Most of the previous work in this area has used Bayesian models. The literature has concentrated on producing smoothed maps of relative abundance, although going from those smooth maps to an abundance estimate would not seem difficult. In particular, Wikle (2002) developed Poisson-lognormal models for a continuous surface, but the counts were at a scale that could be considered points, whereas we have counts in plots with substantial areas. Thogmartin et al. (2004) used a tessellation to create spatial conditional autoregressive (CAR) models (Besag 1974) for neighbors as a spatial random effects model, with the CAR random effects constant within the tessellation, but the model retained point-level data for covariates. Royle et al. (2007) model on a subsample of a systematic grid, and include detection models, but ultimately use a continuous Poisson-lognormal model for counts. Barber and Gelfand (2007) also use Poisson-lognormal models with known covariates to model the intensity surface. Note that none of these methods include a finite population correction factor, and while they are all very attractive, we do not adopt any of them for reasons that we describe next.

GOALS AND ORGANIZATION
Based on this introduction, we desire a total abundance estimator from data like the motivating example that will satisfy several practical conditions. (1) It must be fast to compute, robust, and require few modeling decisions, similar to classical survey methods. Annually, we compute dozens of estimates for data like the example and, depending on the size of the fjord, each may have thousands of photographs. (2) The estimator must use only counts within plots; actual spatial locations of animals are unknown. (3) We are interested in the actual number of seals, not the mean of some assumed process that generated the data. Thus, the estimator must make use of the actual number of seals, and predict for those areas that are unsurveyed. (4) The variance estimator should have a population correction factor that shrinks to zero as the proportion of the study area that gets sampled goes to one. In our real example, we surveyed approximately 26 % of the study area; in classical sampling, it directly reduces the variance by 26 %. Some fjords that we have surveyed are up to 50 % sampled. (5) The estimator should be approximately unbiased (demonstrated through simulations), and we want valid confidence intervals that cover the true number of seals the correct proportion of times; that is, we use this method dozens of times per year and desire confidence intervals in the frequentist sense. (6) It appears that there is nonstationary variance throughout the area, with large areas of zero counts (no seals). A variance estimator that accommodates this will be required. The goals of this manuscript are to develop an estimator to satisfy these criteria.
The rest of this paper is organized as follows. The estimator is developed from models for spatial point processes and generalized linear mixed models, so we begin with a brief review and then develop the estimator in Sect. 2. We provide some simulations to validate the estimator and compare variance estimators in Sect. 3. In Sect. 4, we use the estimator on the motivating example of aerial photographs taken of harbor seals in a glacial fjord in Alaska. We provide some concluding remarks in Sect. 5.

MODEL DEVELOPMENT
From the Introduction and motivating example, note that data are observed at an aggregated support level, but we need a model at the point support level. The reason should be clear; because of the possibility of unbalanced spatial sampling (Fig. 1c), and a real example of it (Fig. 2), we need to predict and then integrate an abundance density surface continuously throughout the unsampled area. To achieve this, we develop a model-based estimator motivated by an inhomogeneous point process (IPP) that has been integrated to yield a Poisson regression model. Part of the attraction of this framework is that it allows inference on point-level support from data on areal support, and then we use the point-level support model to make our abundance estimate. We begin with a brief review of the IPP, describe how abundance is related to the intensity surface, and then draw the connection to Poisson regression.

INHOMOGENEOUS POINT PROCESS MODEL
Assume that locations s = [s x , s y ] of all individuals in A, say S + = (s 1 , . . . , s N ), are the result of an inhomogeneous Poisson process (IPP) with intensity function λ(s|θ) that varies with s, where θ is a vector of parameters controlling the intensity function. The intensity function is defined as where E(·) is expectation, T (R) is the total number of points in planar region R, and |R| is the area of R. In general, when analyzing IPP data, all of the individuals would be located within A, and then inference about θ could be made by maximizing the point process loglikelihood (e.g., Cressie 1993, p. 655). However, this is difficult in our case because A cannot be surveyed in its entirety and individual locations are unknown. For example, a simulated point pattern is shown in Fig. 3, and the plots form a disjointed window on the point pattern that is masked in the areas between the plots, and although we show the point locations within plots, we assume we only have one count for each plot. Ultimately, we will need to estimate the intensity function, but for now we proceed assuming that we have an estimate of the intensity function in the areas between plots, as shown in Fig. 3.

ESTIMATING ABUNDANCE
The primary quantity of interest is the abundance in a particular block A ⊆ R. Because the distribution of individuals is random under the model-based paradigm, there are two types of abundance to consider. First is the expected abundance in A, μ(A) = A λ(u|θ)du, and then there is the realized abundance T (A) for a given realization of S + from λ(s|θ). Assuming an inhomogeneous point process, T (A) ∼ Poi(μ(A)), which is Poisson distribution with mean μ(A). An estimate of the expected abundance isμ(A) = Aλ (u|θ )du, which is often based on plug-in methods from estimatesθ of θ ; i.e.,μ(A) = A λ(u|θ )du.
For an estimate of the realized abundance, consider that the total abundance can be partitioned into observed and unobserved. Assume there are n sample units B i ∈ A, and let B = ∪ n i=1 B i . Note that some B i could be outside A but within R. It is also possible that some B i straddle the boundary of A, though we will not consider that problem here. The region within A that was not sampled is

. Then T (B) is the number of observed points and T (U) the number of unobserved points and T (A) = T (B) + T (U). The total T (A) involves predicting
By substitutingθ for θ , one can use the estimatorμ(U) = U λ(u|θ)du, and, without any observations from U, we use the mean as a predictor T (U) = μ(U). Hence an estimator of the total is for making inference to the realized abundance T (A). Note that, as , so an estimator of this form satisfies condition 3 in Sect. 1.3. Making such inferences involves first estimating the intensity function λ(s|θ), and also incorporating the uncertainty of estimating λ(s|θ). So our immediate goal is to infer λ(s|θ) from data on areal support, which we describe in the next sections.

FROM IPP TO POISSON REGRESSION
Suppose that we have a smooth spatial surface λ(s|θ) that varies with spatial location s and is controlled by parameters θ ; this is the intensity surface. This surface may be integrated over some compact region, such as the plot B i , and this forms the mean of an IPP.
( 2) Now, let s i be the centroid of plot B i . If the area of B i is small compared to the survey area A, and if λ(u|θ) is smooth (i.e., changing slowly within B i ), then Berman and Turner (1992) show that a reasonable approximation for (2) is, where |B i | is the area of B i . This is an important assumption and is part of the general problem of change-of-support; see Gotway and Young (2002), Banerjee et al. (2004, Chapter 6), and Wikle and Berliner (2005). For example, Brillinger (1990Brillinger ( , 1994 shows an early attempt at creating a continuous surface from count data in census tracts. The mean of Y (B i ) can then be modeled with a log link function, forming a GLM with offset log(|B i |), Now we use spatial radial basis functions to model λ(s i |θ). Let s i be the centroid of plot B i . Then where z(s i ) is a vector of covariates at location s i and γ is a parameter vector of fixed effects. The spatial basis functions will form the values of z(s i ).
There has been increasing interest lately in spatial models that use radial basis functions. Suppose there is a set of fixed points in the study area, {κ j ; j = 1, . . . , K }, which are called "knots." Let z(s i ) be a row vector where the jth item contains a radial basis function value C( s i − κ j ) ; ρ). For example, we will use C(h; ρ) = exp(−h 2 /ρ); ρ > 0, which is a Gaussian basis function. A flexible surface is created by taking a linear combination of the radial basis functions. The surface value at location s i depends on parameters γ and ρ as z ρ (s i ) γ , and we attach the subscript to show that values in z depend on ρ. Using radial basis functions can be viewed as a semiparametric approach to spatial modeling (Ruppert et al. 2003), and they have been used for models with non-Euclidean distance measurements (see, e.g., Wang and Ranalli 2007) and for computational efficiency for large data sets (Cressie and Johannesson 2008).
To make the model more flexible, following Cressie and Johannesson (2008), we considered radial basis functions at two scales. Let the "coarse" scale knots be {κ C, Note that Cressie and Johannesson (2008) use 3 scales with approximately 3 times as many knots at the next finer scale. Here, because we only have two scales, we use 4 times as many knots at the finer scale. The knots are generally spread out more or less regularly throughout the study area; more details on an algorithm for knot locations are given in Sect. 2.4.
Consider the log-linear model , and the jth column of Z C has C( s i − κ C, j ) ; ρ C ) as the ith element, and the jth column of Z F has C( s i − κ F, j ) ; ρ F ) as the ith element. We will not consider any covariates in our model, allowing all spatial variation to be modeled through the spatial basis functions, although future development could easily accommodate covariates here. From (4), we only consider an overall constant, x(s i ) β = β 0 . Also, we assume all plots are the same size, |B i | = |B|. Then we can write, where z(s i ) is the ith row of Z. Model (6) will form the basis for estimation and prediction throughout the rest of this paper.

KNOT SELECTION
To place coarse-scale knots, a systematic grid of points was generated within A, and K-means clustering (MacQueen 1967) on the coordinates was used to create K C groups. Because K-means clustering minimizes within-group variance while maximizing amonggroup variance, the centroid of each group tends to be regularly spaced; i.e., it is a spacefilling design that can work well when the region A has an irregular boundary, as in our example data set (Sect. 1.1). We also used K-means clustering placing K F fine-scale knots, but the systematic grid was generated within a minimum convex polygon that contained all nonzero counts intersected with A; this polygon was defined on the centroids of plots with nonzero counts, and an example can be seen in Fig. 3. We found that this helped ensure convergence of the algorithm. If there are too many basis functions with a small range centered in a large area that is all zeros, the fitting algorithm that we describe next would fail to converge. The effect of knot numbers, both K C and K F , are examined in the simulation experiments in Sect. 3. Other methods for spatial knot placement could be used; for example see Nychka and Saltzman (1998). The software PROC GLIMMIX (SAS Institute Inc 2008) generates spatial knots using vertexes of a k-d tree (Friedman et al. 1977). Regarding the number of spatial knots, Ruppert et al. (2003, p. 255) recommend K C + K F = n/4, with no less than 20 and no more than 150.

PARAMETER ESTIMATION
Recall that the ith plot B i is very small in relation to A, and we let s i be the centroid of the ith plot. The count in the ith plot is random, denoted as Y (B i ), and starting from an inhomogeneous Poisson process, from (6) . This is Poisson regression, more generally formed as a generalized linear model (GLM) (McCullagh and Nelder 1989), where , and X and θ are defined after (5). Conditional on fixed ρ values contained in the Z part of X, and iteratively weighted least squares (IWLS) (Nelder and Wedderburn 1972) provides maximum likelihood estimation for θ . Recall that the negative log-likelihood for Poisson regression is (5) with the specific case being (6). Here, we show the dependence of that row on ρ values. An iterative algorithm using block-wise coordinate descent for minimizing the negative likelihood is, • condition on ρ = [ρ C , ρ F ] and use IWLS to estimate θ , • embed the IWLS estimation in a numerical optimization of (8) for ρ.
This optimization routine over just two parameters, ρ, converges quickly and can use existing Poisson regression software for the IWLS update, so it satisfies the speed requirement of condition 1 in Sect. 1.3. To help ensure convergence, we constrained ρ F to be between 0.5 and 3 times the minimum distance between any two knots in {κ F }, and constrained ρ C to be greater than ρ F but less than 3 times the minimum distance between any two knots in {κ C }.
Optimization used the glm() and optim() functions in R (R Core Team 2014), where optim() used the Nelder-Mead optimization algorithm (Nelder and Mead 1965). To ensure boundary conditions, say a as a lower bound and b as an upper bound for one of the elements in ρ, we used a transformation ρ = a + (b − a) exp(ρ * )/(1 + exp(ρ * )), and then optimized for unconstrained ρ * (note that a was a sliding lower boundary for ρ C , but it would stabilize as ρ F found its optimum). Also, note the connection to the Janossy density for IPP (see, e.g., Cressie 1993, p. 655). For some area B with Y ∈ 1, 2, . . . points at locations {s k ; k = 1, 2, . . . , Y } within B, the Janossy likelihood is, From Sect. 2.3, we are assuming that the plots are small enough so that the intensity function is approximately constant within plot, with the intensity value taken from the intensity surface at the centroid of the plot. Using this approximation, then from (9) the negative log-likelihood for all plots is and, when using model (6) for λ(s i |θ , ρ), this makes it apparent that minimizing (8) for θ and ρ is an approximation to maximizing (9). This connection is important because, in Sect. 2.7, we use results from maximum likelihood estimation of the Janossy density in IPP literature to obtain variance estimates. Note that other approaches may be taken, including penalized splines (Ruppert et al. 2003) or Bayesian approaches (see Sect. 1.2).

PLUG-IN ABUNDANCE ESTIMATOR
Denoteθ andρ as the maximum likelihood estimates from Sect. 2.5. Going back to our estimator, recall that T (A) = T (B) + T (U), and we will use our parameter estimates from Sect. 2.5 to obtain the predictor T (U) = μ(U) = U λ(u|ρ,θ )du, where λ(u|ρ,θ ) = exp(xρ(u) θ ). The integral can be approximated with a dense grid of n p points within u j ∈ U, where |U i | is a small area around each u j . We generally assume all |U i | are equal to |U|/n p , yielding a 2-dimensional Riemann integral approximation, which is sufficient if n p is large. Better approaches using numerical integration by quadrature could also be used.

VARIANCE ESTIMATION
The mean-squared prediction error of (1) is Note that as U ∩ A → ∅, then we count all animals in A, and M( T (A)) → 0, so that this estimator satisfies condition 4 in Sect. 1.3. Thus, a finite population correction factor is automatically embedded in the variance estimator. Also,T (U) depends on random counts in B, while T (U) depends on random counts in U. Under the IPP assumption, these will be independent from each other. Further, assume thatT (U) is an unbiased predictor, so For the IPP, var[T (U)] = μ(U), and this is estimated witĥ over the same fine grid of points used in (10).
Define a vector c where the ith element of c is We approximate this integral with where the sum is over a dense grid of n p prediction points in the unsampled area. Using the delta method (Dorfman 1938 A similar result is given by Johnson et al. (2010) in a distance sampling context. Then, as shown by Rathbun and Cressie (1994), ifθ is a maximum likelihood estimator from the Janossy density for the IPP, then an estimator of iŝ Assuming that |B i | = |B| ∀ i is small, this can be approximated aŝ Note that, in (14), variances may become large if the dimension of θ is too high (due to overfitting from too many knots). Through simulations, we will investigate the following variance estimator,M whereμ(U) is given by (12), elements of c are given by (13), andˆ is given by (14). Equation (15) has a nice interpretation by decomposing the variance into the prediction of the total due to fixed intensity surfaceμ(U) (given the regression parameters θ ), plus the variance in estimating the regression parameters θ . Note that we have not taken into account the estimation of ρ. While this would be desirable, we useρ as plug-in estimators for now. This is similar to geostatistical models where covariance parameters are first estimated from the data, and then used for subsequent predictions (see, e.g., Schabenberger and Gotway 2005, p. 263). While this is not ideal, and can be the subject of further research, our simulations show that it has little consequence for the type of data that we analyze.

OVERDISPERSION
Animals (as well as other spatially patterned points) are often clustered at very fine spatial scales. For animals, this might occur as mother-offspring pairs, clustering around locally desirable habitats, etc. The inhomogeneous intensity surface estimated in the foregoing discussion will be unlikely to capture this fine-scale clustering, which will contribute to the overall variance, and without considering it, the confidence intervals on abundance estimates will be too short. Various estimators of overdispersion for count models have been proposed, and the negative binomial and quasi-Poisson are commonly used (e.g., Ver Hoef and Boveng 2007); see Hinde and Demétrio (1998) for an overview. Here, we consider quasi-type models, where, if the mean is φ, then the overdispersion is constant multiplier, ω, so the variance is ωφ. As we demonstrate next, some form of robust estimation or further modeling is required because overdispersion changes through space. In the negative binomial context, robust but nonspatial estimation can be found in Moore and Tsiatis (1991), and nonparametric estimation is found in Gijbels et al. (2010). Our situation is different than general robustness because we want to either trim residuals based on data with low expected values, or downweight them. We describe several estimators next, and compare them in simulations.
Let φ i = E(Y i |β) = g −1 (x iβ ) = exp(x iβ ) be the fitted intensity surface value for the ith plot, where x i is the ith row of X. Denoting y i as the observed value for the ith plot, we considered four different ways to estimate overdispersion: • The traditional estimator: where q is the rank of X.
• A linear regression estimator. Under the Poisson model, the variance is equal to the mean. By regressing the squared residuals against the fitted value, any slope greater than one would be evidence of overdispersion. The linear regression is set up with a zero intercept, so the model is (y i − φ i ) 2 = ωφ i . We used weighted least squares to obtain the estimator, where √ φ i were the weights. Notice that generally, this may not be a desirable estimator. Values with small expectations have virtually no effect on the slope, whereas values with larger expectations will have a great deal of leverage. In our case, this is a desirable feature, as discussed earlier. In fact, we create additional weight for values with large expectation by using √ φ i .
• Estimator based on a trimmed mean of squared Pearson residuals from the upper quantile of fitted values. Let F = {φ 1 , φ 2 , . . . , φ n } be an unordered set of expected values for the n observed counts, and {φ (1) , φ (2) , . . . , φ (n) } be the set of ordered values, from smallest to largest, where φ (1) = min(F) and φ (n) = max(F). Also, if φ (i) = φ j , then y (i) = y j ; that is, the observed values are ordered by their fitted values as well. Let 0 ≤ p < 1 be some proportion, then where x rounds x down to the nearest integer. That is, the proportion p of the squared Pearson residuals with the lowest fitted values are trimmed from the overdispersion computation.
Examples of the overdispersion estimators are shown in Fig. 4, which were taken from the data seen in Fig. 3. The traditional estimator can be viewed as a constant fit (the average value) through all of the squared Pearson residuals for all fitted values, so this is shown as a horizontal solid line, the one that is below the short-dashed line (whose value is constant at one) in Fig. 4a. Note especially the wide divergence in squared Pearson residuals for low expected values. This is not surprising because we are dividing by very small numbers, so any count greater than zero will have a very large residual. This instability, along with the fact that these values do not really contribute much to overall abundance, leads to estimator ω TG ( p). Here, we trim off the lowest expected values. Trimming off the lowest 75 %, and averaging the rest, can be viewed as a constant fit (horizontal line) through the squared Pearson residuals for the upper 25 % of fitted values, and is shown as the long-dashed horizontal line that is above the short-dashed line in Fig. 4a. The other idea is to treat raw squared residuals, (y i − φ i ) 2 as a response variable in a zero-intercept regression, where the predictor variable is the fitted value φ i . This is shown as the solid line in Fig. 4b, which is above the one-to-one line. The estimated slope of this line is taken as the overdispersion estimate. Similar to trimming in ω TG ( p), the regression estimator ω WR ( p) downweighs residuals with small expected values by forcing the line through zero, and it eliminates division by very small numbers. In fact, we considered weighted regression to add even more weight to higher fitted values. After some trail and error, we used weights √ φ i , but this is clearly an area for further research.
With these three overdisperson estimators, we have several variance estimators of the abundance estimator (1) at our disposal, where k = OD, WR, or TG, and ω TG has the additional trimming parameter p. We include one more estimator using the same logic applied to the IPP variance estimator as the trimmed overdispersion estimator ω TG ( p). If φ np represents the smallest fitted value summed in ω TG , then we computed (14) using only those i sites whose values satisfied exp(xρ(s i ) θ ) > φ np . Let us call this , which when substituted into (15) and combined with ω TG yields var(T (A)) TL ≡ ω TG (μ(U) + c c).
For confidence intervals, note that from (10), the estimate is a sum of a large number of lognormal variates. That is, we can assume thatθ are normal because they are maximum likelihood estimates. If each summand in (10) was independent, then (10) would converge to normality because of the central limit theorem, but due to correlation, the distribution is unknown and may be asymmetric. We investigated this by simulating (10) usingˆ from (14) as estimated from various data sets. In all cases, (10) was skewed, and a log transformation made the distribution approximately normal. Thus, we recommend computing confidence intervals on the log scale, and then back-transforming. Using the delta method (Dorfman 1938;Ver Hoef 2012), an approximate 100(1 − α) % level confidence interval is for k = OD, WR, TG, or TL, where z α/2 is the upper α/2 percentage point of a standard normal distribution. Note that this also yields the desirable property that the lower bound of the confidence interval is always greater than 0.

SIMULATION EXPERIMENTS
We simulated data under four different conditions to examine the performance of the abundance estimator and the variance estimators of abundance. In all experiments, data were simulated in the region A = {(x, y) : x ∈ [0, 10] × y ∈ [0, 10]}. For each experiment, data sets were simulated 1,000 times, with the number of points and their spatial locations changing, but the sample units were held fixed as shown in Fig. 5.

EVALUATING THE EXPERIMENTS
For each experiment described below, the expected number of simulated points was near 1,000. Let T t be the actual number of points simulated in the tth simulation, letT t be the estimator of the total from (10), and letv k,t be a variance estimator given in (16) or (17) for k = O D, W R, T G, or T L, for the tth simulation. The performance of the abundance estimator was evaluated in three ways: • Bias for an experiment was computed as, 1 1000 1000 t=1T t − T t .
• Root-mean-squared prediction error (RMSPE) was computed as, • Coverage of 90 % confidence interval (CI90) was computed as, where I (·) is the indicator function, equal to one if the argument is true, otherwise it is zero. Note that (CI90) can be computed for each k in (16) and (17), which we denote as CI90 k in the tables that summarize the experiments.
For all experiments we used ω TG (0.75), but investigated the effect of p in Experiment 4. We also investigated knot density by changing K C and K F , but always using the algorithm described in Sect. 2.4. For each experiment, we included bias, RMSPE, and CI90 for simple random sampling (SRS), as described in the Introduction. We realize that SRS is inappropriate for these data, but it provides a convenient benchmark for comparison.

EXPERIMENT 1: INHOMOGENEOUS SPATIAL POINT PROCESS WITH REGULAR SAMPLING
The study area A was a square starting at (0,0) and 10 units on each side. Data were simulated using rejection sampling. Consider the intensity surface λ(x, y) = x/20 + y/20, Table 1. Results for bias, RMSPE, confidence interval coverage, and failure rate for simulation experiment 1. Fail rate f 0.000 0.000 0.000 0.000 0.000 a 90 % confidence interval coverage based on standard errors without overdispersion. b 90 % confidence interval coverage using classical overdispersion. c 90 % confidence interval coverage using a weighted regression overdisperion estimator. d 90 % confidence interval coverage using a global trimmed mean overdisperion estimator. e 90 % confidence interval coverage using a local trimmed mean overdisperion estimator. f Failure of the estimator due to lack of convergence or excessively large estimates or standard errors.
The number of coarse-scale knots used is given by K C , and K F is the number of fine-scale knots. An example of a single simulated data set is given in Fig. 5.
which increases linearly from 0 at (0,0) to 1 at (10,10) within A. A location was simulated with x * , y * , and z * , where each was drawn from Unif(0, 1). If z * < λ(x * , y * ), then the simulated location was retained. The set (x * , y * , z * ) was drawn 2000 times, so the expected number of locations retained was 1,000 per simulated data set, but note that the actual number varied randomly among simulations. For each simulated data set, sample units as square plots that measured 0.3 on a side were systematically placed in a 16 × 16 grid as shown in Fig. 5. The 256 plots covered 23.04 % of the study area A.
The results of experiment 1 are given in Table 1. Note that for SRS and various knot proportions, there was little bias, which is <1 % of the average total. All methods had very similar RMSPE as well, and 90 % confidence intervals generally had the appropriate coverage, no matter which overdispersion method was used. Of course, the data were not simulated with overdispersion. The only exception occurred when using many knots for the trimmed overdispersion estimators CI90 TL and CI90 TG , where variance was overestimated.

EXPERIMENT 2: INHOMOGENEOUS SPATIAL POINT PROCESS WITH IRREGULAR SAMPLING
For simulation experiment 2, locations were simulated exactly as they were in Experiment 1 (Sect. 3.2). However, we created unbalanced spatial sampling by removing one column and two rows of sample units, as shown in Fig. 5. In this case, 210 plots covered 18.9 % of the study area A. Table 2. Results for bias, RMSPE, confidence interval coverage, and failure rate for simulation experiment 2. Details on column and row labels are given in Table 1 and the text. An example of a single simulated data set is given in Fig. 5. Row names are described in Table 1.
The results of experiment 2 are given in Table 2, which should not be surprising for SRS. Indeed, the goal of this research was to find good estimators when high (or low) abundance areas were oversampled, and SRS makes no weighting adjustments for this. Consequently, it had a large bias, whereas T (A) in (10) remained relatively unbiased with RMSPE only slightly larger than experiment 1 (note, too, that sample sizes were smaller here). The same basic patterns appeared for CI90, with generally valid confidence intervals (except for SRS), except when many knots are used for the trimmed overdispersion estimators CI90 TL and CI90 TG .

EXPERIMENT 3: DOUBLE SPATIAL CLUSTER PROCESS ON INNER RECTANGLE
One difficult situation for estimating abundance from counts occurs when there is a large number of zeros and there is fine-scale clustering, so we tested the abundance estimator under both of those conditions. For this experiment, seed points were simulated within a rectangle within the study area. Again, the study area A was a square starting at (0,0) and 10 units on each side. The inner rectangle itself had random boundaries, where the lower x-axis and y-axis boundaries were each randomly drawn from Unif(3.5, 4.5), and the upper x-axis and y-axis boundaries were each drawn from Unif(7.5, 8.5). Next 100 parent seed points were uniformly simulated over the inner rectangle (such a rectangle is shown with solid lines in Fig. 3), and then each parent had a random number, Poi(15), of children that were uniformly distributed on a square with sides of length 2 centered on each parent. A second finer scale cluster process was added by creating 25 more parent points uniformly distributed over the inner rectangle, where each parent had a random number, Poi(9), of children that were uniformly distributed on a square with sides 0.4 centered on each parent. After simulating all points, they were thinned using the same function as experiments 1 and 2, by simulating z * i ∼ Unif(0, 1) at each simulated location with coordinates x i and Table 3. Results for bias, RMSPE, confidence interval coverage, and failure rate for simulation experiment 3.
Knots Details on column and row labels are given in Table 1 and the text. An example of a single simulated data set is given in Fig. 5. Row names are described in Table 1.
y i , and keeping that location if z * i < x i /20 + y i /20. From 1,000 simulated experiments, this yielded an average of 1034 points per simulation. An example of one simulation is given in Fig. 5. The sample units were placed in the same positions as in Experiment 2 (Sect. 3.3).
The results of experiment 3 are given in Table 3. Once again, because of the unbalanced sampling, SRS had a large RMSPE and was highly biased, whereas T (A) in (10) remained unbiased at generally less than 0.5 % of the total. RMSPE was larger than experiments 1 and 2, but this is not surprising given the smaller area with positive count values. This experiment showed some poorer performance for CI90, especially for CI90 OD , which underestimated variance with coverage near 80 % rather than 90 %. Both CI90 WR and C I 90 TG performed more poorly with increasing numbers of knots. CI90 TL coverage is about 3 % low for the fewest number of knots, but generally improves with the number of knots. Notice also that there is almost a 2 % chance that the parameter estimation algorithm will fail when there are many knots.

EXPERIMENT 4: DOUBLE SPATIAL CLUSTER PROCESS ON DOUBLE INNER RECTANGLES
For this experiment, the seed points were simulated in a manner similar to experiment 3. However, the sample units were made smaller, with length 0.14 on a side, on a 26 × 26 grid in a study area, A, which was again a square starting at (0,0) and 10 units on each side. This time there were two inner rectangles. One inner rectangle had random boundaries where the lower x-axis and y-axis boundaries were each randomly drawn from Unif(5.8, 6.2), and the upper x-axis and y-axis boundaries were each drawn from Unif(7.8, 8.2). Next 75 parent seeds were uniformly simulated over this rectangle, and each parent seed had a random number, Poi (14), of children uniformly distributed in a box with sides of length 2 centered on each parent. A finer scale cluster process was added by creating 25 more parent points Table 4. Results for bias, RMSPE, confidence interval coverage, and failure rate for simulation experiment 4. Details on column and row labels are given in Table 1 and the text. An example of a single simulated data set is given in Fig. 5. Row names are described in Table 1. uniformly distributed over this inner rectangle, where each parent had a random number, Poi(8), children that were uniformly distributed on a square with sides 0.4 centered on each parent. A second inner rectangle had random boundaries where the lower x-axis boundary was randomly drawn from Unif(0.8, 1.2) and the upper x-axis boundary was drawn from Unif(3.8, 4.2), while the lower y-axis boundary was randomly drawn from Unif(4.8, 5.2), and the upper y-axis boundary was drawn from Unif(7.8, 8.2). Here, 25 parent seeds were uniformly simulated over this rectangle, and each parent seed had a random number, Poi (14), of children uniformly distributed in a box with sides of length 1 centered on each parent. A finer scale cluster process was achieved by creating 10 more parent points uniformly distributed over this inner rectangle, where each parent had a random number, Poi(8), of children that were uniformly distributed on a square with sides 0.4 centered on each parent. After simulating all points, they were thinned using the same function as experiments 1-3, by simulating z * i ∼ Unif(0, 1) at each simulated location with coordinates x i and y i , and keeping that location if z * i < x i /20+ y i /20. From 1,000 simulated experiments, this yielded an average of 1012 points per simulation. An example of one simulation is given in Fig. 5. We created unbalanced spatial sampling by removing one column and two rows of sample units, as shown in Fig. 5. In this case, 600 plots covered 11.8 % of the study area A.
The results of experiment 4 are given in Table 4. Once again, because of the unbalanced sampling, SRS had a large RMSPE and was highly biased, whereas T (A) in (10) remained unbiased for smaller number of knots, but was >1 % of the total for the largest number of knots. RMSPE was less than experiments 3, likely due to a larger area with nonzero counts and smaller plots, leading to a larger sample size. Once again, CI90 OD underestimated variance with coverage near 83 % rather than 90 %. CI90 WR was about 4 % high for small number of knots, but improved with numbers of knots. CI90 TG remains about 2 % high for all combinations of knot numbers, while CI90 TL is 1-2 % low for all combinations of knot numbers. Here, notice also that there is over 24 % chance that the parameter estimation algorithm will fail when there are many knots. The effects of varying p in CI90 TG and CI90 TL are shown in Fig. 6, using 1,000 simulated data sets as described for experiment 4. Notice that confidence coverage for CI90 TG was in the "valid" zone between p = 0.4 and p = 0.8, and CI90 TL was in the "valid" zone when p ≥ 0.5. In a real setting, such a p value must be chosen by data examination (at least without further research). Looking at the simulated data in Fig. 5, one would try to estimate the area dominated by zero counts, and then trim them. A p value anywhere from 0.6 to 0.8 seems reasonable, and would lead to nearly correct confidence intervals using either CI90 TG or CI90 TL . Clearly, CI90 TG is a more conservative strategy, but trimming aggressively with CI90 TL appears viable. Another consideration is that we expect that smaller p values would be more efficient by using more data.

EXAMPLE: AERIAL SURVEYS OF HARBOR SEALS
The study area boundary, the locations of all 2,080 plots, and the observed counts of seals within plots are shown in Fig. 2, and the data were summarized in Sect. 1.1. We used K C = 4 knots in the coarse grid and we used K F = 15 knots in the fine grid shown in Fig. 7. Notice that the fine grid of knots is contained in a bounding rectangle around only those plots with nonzero counts; with the reason explained in Sect. 3.1. The model fit took 17.56 s on a Intel Xeon 2.66 GHz processor running under the linux operating system. The estimated range parameter ρ F was 1.81 km, and ρ C was 4.04 km. The fitted intensity surface is shown in Fig. 7. The estimate obtained from integrating this surface, along with the observed count, using T (A) in (10)  We tried different combinations of knots and p in ω TG and ω TL . No systematic attempt is made to present those results, and knot selection and overdispersion is a topic for future research. However, we did note that the abundance estimate was little changed for knots up to K C = 8 and K F = 32. However, there were changes in overdispersion estimates. Also, we noticed that for the knots that we selected, increasing p in ω TG and ω TL decreased overdispersion, rather than the increasing values seen in Fig. 6. Conceptually, it is possible that an area with high counts could have less variability than a larger area that includes both high counts and low counts. Also, note that ω OD = 17.26, which is much larger than one, and larger than all other overdispersion estimators, in contrast to the example provided by Fig. 4. However, the explanation is also provided by Fig. 4 because for this data set there were a few nonzero counts with very low expected values that dominated ω OD . We chose this example in part because it showed some exceptions for the overdispersion parameters. Based on dozens of similar real examples, most of them have ω OD = 1, and this appears to be an unstable estimator for our purposes.

DISCUSSION AND CONCLUSIONS
Our objectives were to develop a model-based estimator of a total by using counts from irregularly spaced plots. We wanted this estimator to have goals, properties, and performance similar to classical sampling: to estimate the realized total, not the mean of an assumed process, to have a finite population correction factor, to be unbiased with valid confidence intervals that must be robust to nonstationary mean and variance, and to be fast to compute. The problem was made more difficult by features of the data that were seen in the real example. First, the nonzero counts were highly clustered in space, and there were large areas of zeros. Secondly, the counts, where they occurred, showed overdispersion. Finally, sample sizes were quite large, so we needed to use computationally efficient methods. The general approach that we took was to assume that the data came from an inhomogeneous point process with overdispersion. We modeled the intensity surface of the inhomogeneous point process using spatial basis functions in a generalized linear mixed model framework.
To test the method, we started with fairly benign conditions in experiment 1. We then added complexity to the simulations that matched the complexity seen in the real data, by simulating spatially unbalanced sampling, data with overall trend in point density, several areas of clustered points, overdispersion within the cluster areas, and yet large areas with zeros, culminating in experiment 4. Overall, our method worked well, especially when using some of the newly introduced overdisperion estimators. One of the main contributions of this manuscript was the introduction of overdispersion estimators when the overdispersion appears to be varying spatially (a nonstationary overdispersion).
One of the interesting findings from our research was the effect of knot placement and proportion (Tables 1, 2, 3, and 4). Initially, we found a lack of convergence when fitting these models if there were too many knots, with short ranges, over areas containing all zeros. In retrospect, that may not seem surprising, but we have not found it reported in the literature. For that reason, we created spatially restricted knots over areas with nonzero values (dashed line in Fig. 3; also see Fig. 7). For placing knots, we used a K-means clustering algorithm on the spatial coordinates to create a space-filling design. Other approaches that could be used are mentioned in Sect. 2.4. Our results show that the bias does not depend very much on the number of knots, but the standard errors are quite sensitive. An encouraging result is that standard errors yield better confidence intervals when the number of knots is quite small, making the algorithm very fast. The whole issue of knot selection needs further research.
Besides more research on overdispersion estimators, and knot placement and number, there are numerous modifications that could be applied to our basic approach. We chose spatial basis functions that were Gaussian kernels at two scales, and there are many obvious modifications. Because we annually analyze dozens of data sets like our example for harbor seals, we took a maximum likelihood approach to estimate parameters and total abundance; however, Bayesian approaches could also be used (Wikle 2002;Christensen and Waagepetersen 2002). The ability to estimate the intensity function, essentially pointlevel information, from data at an aggregated level, i.e., counts from plots, depends on the plots being small in relation to changes in the intensity function. If plots are very large, then methods in this paper could still be used, but spatial points would need to be mapped within plots, and more traditional methods from the spatial point process literature could be used for estimating the intensity surface. For example, the R (R Core Team 2014) package spatstat (Baddeley and Turner 2005) can be used when the realized point patterns have a complicated "mask" comprised of many disjoint plots (Baddeley and Turner 2006). Alternatively, area-to-point geostatistical methods (Kyriakidis 2004) could be used. The main point is that once the intensity surface is estimated, (10) can be used to estimate total abundance. The variance of the total abundance can be estimated with (15) if maximum likelihood methods are used, along with one of the overdispersion factors (Sect. 2.8) that make sense for the problem under consideration.