Modeling spatiotemporal abundance of mobile wildlife in highly variable environments using boosted GAMLSS hurdle models

Abstract Modeling organism distributions from survey data involves numerous statistical challenges, including accounting for zero‐inflation, overdispersion, and selection and incorporation of environmental covariates. In environments with high spatial and temporal variability, addressing these challenges often requires numerous assumptions regarding organism distributions and their relationships to biophysical features. These assumptions may limit the resolution or accuracy of predictions resulting from survey‐based distribution models. We propose an iterative modeling approach that incorporates a negative binomial hurdle, followed by modeling of the relationship of organism distribution and abundance to environmental covariates using generalized additive models (GAM) and generalized additive models for location, scale, and shape (GAMLSS). Our approach accounts for key features of survey data by separating binary (presence‐absence) from count (abundance) data, separately modeling the mean and dispersion of count data, and incorporating selection of appropriate covariates and response functions from a suite of potential covariates while avoiding overfitting. We apply our modeling approach to surveys of sea duck abundance and distribution in Nantucket Sound (Massachusetts, USA), which has been proposed as a location for offshore wind energy development. Our model results highlight the importance of spatiotemporal variation in this system, as well as identifying key habitat features including distance to shore, sediment grain size, and seafloor topographic variation. Our work provides a powerful, flexible, and highly repeatable modeling framework with minimal assumptions that can be broadly applied to the modeling of survey data with high spatiotemporal variability. Applying GAMLSS models to the count portion of survey data allows us to incorporate potential overdispersion, which can dramatically affect model results in highly dynamic systems. Our approach is particularly relevant to systems in which little a priori knowledge is available regarding relationships between organism distributions and biophysical features, since it incorporates simultaneous selection of covariates and their functional relationships with organism responses.


| INTRODUC TI ON
Understanding how the spatial distribution and abundance of an organism responds to biophysical features is fundamental to many aspects of ecology and conservation (Schröder & Seppelt, 2006).
Since continuous sampling of the entire range or population of a species is usually impossible, distribution mapping typically involves a series of steps that include surveying a representative subset of the area or population of interest at various time periods, fitting models to represent the relationships of observed data to environmental covariates, using these models to predict utilization of un-sampled areas or time periods based on biophysical habitat features, and finally validating predictions with onthe-ground observations (Borchers, Buckland, & Zucchini, 2002;Certain & Bretagnolle, 2008;Kinlan, Menza, & Huettmann, 2012;Nur et al., 2011). Although such model-based approaches are widely used, their implementation requires addressing complex statistical challenges (Guisan & Thuiller, 2005), particularly for mobile organisms whose distributions and habitat requirements may vary in space and time (Runge, Martin, Possingham, Willis, & Fuller, 2014).
Spatiotemporal variability and uncertainty surrounding both the distribution of a species and key environmental covariates represent frequent challenges to the development of predictive models. Landscape-or population-scale occupancy models may lack sufficient resolution to accurately address small-scale spatial variation in habitat use; conversely, small-scale models may be too precise to apply across landscapes (Johnson, 1980;Johnson, Nielsen, Merrill, McDonald, & Boyce, 2006;Johnson, Seip, & Boyce, 2004). Error can be introduced by spatial or temporal mismatches between occurrence estimates, environmental variables, and questions of interest (Austin & Van Niel, 2011;Guisan & Thuiller, 2005;Mainali et al., 2015). In addition to variability, uncertainty surrounding the biotic and abiotic factors determining the distribution of a species can also limit the development and implementation of model-based distribution estimates (Thuiller, 2004). Because a priori knowledge of how occupancy and abundance relate to biophysical features is often lacking, survey data themselves can be used to select key environmental covariates (Guisan & Thuiller, 2005). This selection process requires choosing appropriate habitat variables from among a suite of intercorrelated covariates while avoiding overfitting (Hoeting, 2009;Merow et al., 2014). Most predictive models involve assumptions about the form of the response function between the occurrence or abundance of an organism and individual biophysical features.
However, the information needed to inform these assumptions is typically unknown prior to analysis, which may lead to poor model performance (Austin, 2007;Mainali et al., 2015). Temporal variation in both the distribution and habitat preferences of a species can introduce further uncertainty, because organisms' responses to changes in habitat conditions may not be instantaneous and may vary across the annual cycle (Selonen, Varjonen, & Korpimäki, 2015;Yamanaka, Akasaka, Yamaura, Kaneko, & Nakamura, 2015).
Furthermore, both occupancy and abundance may respond not only to biophysical habitat features, but also to the distribution of other organisms such as conspecifics, competitors, predators, or prey (Blackburn, Cassey, & Gaston, 2006;Guisan & Thuiller, 2005).
Aside from their ecological complexity, survey data can present several statistical challenges to modeling and interpretation. Surveys can be modeled using occupancy (presence/ absence) or abundance (numerical) approaches, which measure different aspects of habitat use and have different distributions and response functions. Given the additional statistical complexity involved in interpreting count data, abundance estimates are often overlooked when mapping the distribution of organisms (He & Gaston, 2000); however, occupancy estimates alone may provide incomplete or misleading information about habitat quality (Pulliam, 2000). Variation in abundance is often a key component of a species' response to habitat quality and is crucial for accurately predicting species distributions (Howard, Stephens, Pearce-Higgins, Gregory, & Willis, 2014;Johnston et al., 2015). The statistically challenging features of count dataparticularly overdispersion, in which the variance of the data exceeds the mean-may in fact represent important biological responses to environmental features and conditions (McMahon, Purvis, Sheridan, Siriwardena, & Parnell, 2017;Richards, 2008).
Modeling count data also often requires accounting for zero-inflation (Martin et al., 2005), non-linear responses to covariates (Cunningham & Lindenmayer, 2005), and spatial and temporal autocorrelation (Hoeting, 2009), which require a highly flexible modeling approach with few assumptions about either underlying distribution or response functions.
Generalized additive models (GAMs: Hastie & Tibshirani, 1990) and their extension, GAMs for location, scale, and shape (GAMLSS: Rigby & Stasinopoulos, 2005) offer several features that make them well-suited for modelling complex, uncertain, or variable relationships between survey data and environmental covariates. GAMs do not assume linear effects on the response but flexibly adapt to the observed data, which makes them especially applicable to systems in which the form of the relationship between species occupancy, abundance, and underlying environmental conditions K E Y W O R D S distribution, GAM, gradient descent boosting, hurdle models, Nantucket Sound, sea ducks, surveys is often non-linear and unknown (Guisan & Zimmermann, 2000).
Moreover, unlike other generalized modeling approaches, GAMLSS allow both the mean and dispersion of the response to be modeled as a function of environmental covariates (Rigby & Stasinopoulos, 2005), which incorporates additional information about count data not reflected by mean values alone (McMahon et al., 2017). Despite these promising features, although GAM has recently gained popularity as a predictive distribution modeling approach (Miller, Burt, Rexstad, & Thomas, 2013), GAMLSS have yet to be widely adopted for modeling the spatial distribution of species based on biophysical features.
Our approach independently evaluates environmental covariates for both occupancy and abundance, while allowing response functions to vary. We generate a single distribution estimate based on both occupancy and abundance that can be applied to environments with high levels of spatiotemporal variation and uncertainty. As a case study, we applied our proposed modeling framework to sea ducks (tribe Mergini) in Nantucket Bay, MA, USA. Understanding the winter habitat use and distribution of sea ducks in southern New England is crucial for the siting of proposed offshore wind farms (Bradbury et al., 2014;Langston, 2013); however, the implicitly high spatial and temporal variability of sea duck distributions, as well as poor understanding of habitat factors driving temporal and spatial variation in their distribution, has previously limited fine-scale prediction of habitat use in the proposed construction area (Bowman, Silverman, Gilliland, & Leirness, 2015). We demonstrate an application of our modeling framework to informing conservation planning in the face of both high variability and ecological uncertainty by developing models from systematic aerial survey data of sea ducks in this system.

| MATERIAL S AND ME THODS
Our predictive approach incorporates five distinct methodological steps: (a) survey data collection, (b) separation of presence from abundance, (c) integration of environmental covariates, (d) synthesis of presence and abundance models, and (e) validation, which correspond to numbered sections in the model schematic ( Figure 1). We describe these five steps sequentially below, along with details of how we applied the modeling process to our case study of sea ducks in Nantucket Sound. All analyses were conducted in R (R Core Team, 2018) with the add-on packages gam-boostLSS , mboost (Hothorn, Buehlmann, Kneib, Schmid, & Hofner, 2017), and stabs ).

| Study system
We conducted fieldwork throughout Nantucket Sound in Massachusetts, USA ( Figure 2). Our study area encompassed ca.
1,500 km 2 of Nantucket Sound, was relatively shallow (generally <20 m deep), and included some of the most important sea duck wintering habitat in the western Atlantic (Silverman, Saalfeld, Leirness, & Koneff, 2013;White, Veit, & Perry, 2009). The primary species of sea ducks found in Nantucket Sound were Common Eider (Somateria  We flew at an average altitude of 152 m and speed of 167 km/hr (90 kts), the slowest speed at which the aircraft could safely fly. This altitude allowed us to identify most birds at the sea surface and reduced disturbance (i.e., flushing birds to another part of the study area and potential double counting). We conducted surveys only on days with wind speeds ≤15 kts and good visibility (>15 km). Surveys had an average duration of ~2.5 hr and occurred between 0900 to 1600 hours to ensure that birds had completed any post-dawn movements (Davis, 1997) but had yet to initiate pre-sunset movements from feeding to roosting areas; this time window also reduced glare for observers due to low sun angles.

| Survey design
On each survey flight, two observers used their unaided eyes to continuously detect individuals or flocks, identified sea ducks to species with the aid of binoculars as needed, and communicated the number, species, and behavior (on the water or flying) of observed ducks to a recorder who entered georeferenced locations using dLOG (Ford, 1999). Observers monitored the sea surface on their side of the plane in a ca. 91 m-wide transect between ca. 56 and 147 m from the plane. The narrow strip width ensured birds were detectable and identifiable with the naked eye and limited situations in which ducks were too abundant or spread over too wide an area to count accurately. We attempted to limit perception bias (i.e., to miss few individuals present to be counted; Marsh & Sinclair, 1989) by using low flight altitudes and narrow transect widths (Buckland et al., 2012;Certain & Bretagnolle, 2008); however, our survey methods did not address potential undercounting of individuals that were diving during flyovers, and therefore may not have been present to be detected. Transect dimensions resulted in the sampling of ~6% (median; 68.4 km 2 ) of the study area during a survey.
Due to the difficulties associated with identifying to species the three species of scoters, we pooled all scoter observations (hereafter, scoters), while we modeled Common Eider and Long-tailed Duck as separate species. While using pooled data from multiple scoter species reduces our ability to make inferences about species-specific ecology and habitat use, scoters overlap broadly in shared wintering habitat across the study region and are generally subject to common conservation and management regimes. We subsequently consolidated counts for each species (eider and Long-tailed Duck) or species group (scoters) into 2.25 km 2 segments ( Figure 2); this resolution (1.5 km × 1.5 km) corresponded approximately to the coarsest level of resolution of biophysical covariates (see below).

| Presence and abundance
We related spatiotemporal variation in sea duck occupancy (i.e., probability of presence) and abundance to potentially relevant biophysical and spatiotemporal covariates. Because we observed a high incidence of zero counts (e.g., no eiders were detected in 75% of segment observations), and we assumed our survey design led to a low incidence of false zeros in study segments (Certain & Bretagnolle, 2008), we applied a negative binomial hurdle model (Manté, Kidé, Yao-Lafourcade, & Mérigot, 2016). This approach allowed us to model presence/absence in all grid cells, and abundance only in cells where at least one individual was detected.
We first used a logistic regression model to represent the probability of occurrence of at least one individual (hereafter, the occupancy model) in a given segment (Figure 1:2). We then used a truncated negative binomial model to represent the abundance of sea ducks in that segment conditional on their presence (hereafter, F I G U R E 1 Schematic diagram of our modeling approach: (1) conducting initial transect surveys; (2) extrapolating occupancy probability and conditional abundance for each grid cell; (3) modeling relationships between occupancy, abundance, and habitat variables, (4) estimating unconditional abundance based on habitat characteristics, and (5) generating predictive estimates of occupancy and abundance over the full region. Detailed methodology for each step is provided in the corresponding numbered subsections in Section 2 F I G U R E 2 Aerial strip transect tracks (gray lines) conducted during winter (October-April, 2003 sea duck surveys (n = 30) in Nantucket Sound, Massachusetts, USA. The grid indicates the extent of the 1,100 km 2 study area and its division into 504 2.25 km 2 segments. The polygon (thick black line) in northwest Nantucket Sound indicates a 62 km 2 permitted wind energy development on Horseshoe Shoal TA B L E 1 Biophysical and survey covariates used to evaluate the distribution and abundance of Common Eider, Black, Surf, and Whitewinged Scoter, and Long-tailed

| Covariates
Distribution of large marine vertebrates is primarily a function of the distribution of preferred prey items. Since we did not have direct measurements of the availability of sea duck benthic prey (e.g., mollusks and crustaceans), we evaluated biophysical covariates expected to influence the distribution and abundance of these organisms, including water depth, sediment grain size, and primary productivity (Table 1). Additionally, we included interactions with time that allowed the effects of two ecological covariates (water depth and relative sea surface temperature) and all spatial covariates to vary over time within a given winter. We standardized (i.e., mean centered and scaled) all continuous covariates.

| Modeling approach
To relate occupancy and abundance data to environmental covariates, we used additive models (Figure 1:3). We implemented a GAM for occupancy that flexibly accommodated varying effects of covariates on presence/absence data (Hastie & Tibshirani, 1990;Wood, 2006). For abundance data, we used a GAMLSS approach (Rigby & Stasinopoulos, 2005). Using GAMLSS allowed us to independently model dispersion of count data in relation to biophysical features, as well as account for potential non-linear or heterogeneous responses of abundance to underlying environmental covariates (Stasinopoulos & Rigby, 2007).
We fitted GAM and GAMLSS using an iterative machine-learning approach, component-wise functional gradient descent boosting (Bühlmann & Hothorn, 2007;Hothorn et al., 2010;Mayr et al., 2012;Hofner, Boccuto, & Göker, 2015; in a cyclical framework (Thomas et al., 2018). The first step of this process was to compute the negative gradient of a pre-selected loss function, which acts as a working residual by giving more weight to observations not properly predicted in previous iterations. We used the binomial log-likelihood as the loss function for occupancy models and the truncated negative binomial log-likelihood as the loss function for count models. For GAM, we computed the negative gradient of the mean only. For GAMLSS, we computed the negative gradient separately for mean and overdispersion in each iteration while holding the other parameter as a fixed constant (Mayr et al., 2012).
In the next step, we fitted various functional forms of each covariate relative to each response, called base-learners (Hofner, Mayr, Robinzonov, & Schmid, 2014) to the negative gradients of the models. For each continuous covariate, we specified two possible base-learners: a linear base-learner and a base-learner for the smooth deviation from the linear effect via penalized splines (i.e., P-splines; Eilers & Marx, 1996;Schmid & Hothorn, 2008). This allowed the model to select the best alternative for each covariate between no effect, linear effect, and smooth effect.
For categorical covariates, we coded the categories and used a separate linear base-learner for each category, excepting a reference category (i.e., dummy-coding). To address potential spatial autocorrelation, we included a smooth surface function of the  spatial coordinates of segment centers (Kneib, Müller, & Hothorn, 2008), which accounted for underlying variance in sampling units similarly to a random effect term in a linear model. This surface comprised four base-learners-linear base-learners for the easting and northing, their linear interaction, and a penalized nonlinear tensor product (Kneib et al., 2008;Kneib, Hothorn, & Tutz, 2009;Maloney, Schmid, & Weller, 2012). We also allowed this surface to vary over time within a winter via an interaction. The decomposition of continuous covariates into linear and penalized nonlinear base-learners reduced bias and overfitting by preventing the preferential selection of smooth base-learners (Hofner, Hothorn, Kneib, & Schmid, 2011;Kneib et al., 2009). Thus, we restricted each base-learner to a single degree of freedom and omitted the intercept term from each base-learner (Hofner et al., 2011;Kneib et al., 2009) and added a linear base-learner to the overall model to represent the model intercept. Once all potential base-learners had been tested, the single best fitting base-learner was added to the current model fit. As only the single best-fitting base-learner was selected in each iteration, the algorithm integrated intrinsic selection of the most relevant covariates and their functional form.
In order to maximize predictive accuracy while avoiding model overfitting, we employed an early stopping mechanism (Maloney et al., 2012;Mayr et al., 2012) during variable selection by stopping the algorithm prior to convergence to maximum likelihood estimates. In other words, after adding the best-fitting base learner to the model, the negative gradient was then reevaluated at the current model fit and the procedure of testing, comparing, and adding base-learners ( Figure 1:3) was repeated until a pre-specified number of iterations was reached (Bühlmann & Hothorn, 2007). We used 25-fold subsampling to determine the optimal stopping iteration for each model.
Specifically, we randomly drew (without replacement) 25 samples of size n/2 from the original data set. We used the selected sample to estimate the model and the balance of the data in each sample to determine the out-of-bag prediction accuracy (empirical risk) measured by the negative log-likelihood of each model; the optimal stopping iteration (m stop ) is the iteration with the lowest average empirical risk.
In boosted GAMLSS models we used multi-dimensional subsampling to determine the stopping iteration for each of the GAMLSS parameters while allowing for potentially different model complexities in the parameters; a detailed explanation of this cross-validation (subsampling) scheme is given in Hofner, Mayr, et al. (2016).
Since boosting methods typically produce "rich" models relying to some extent on many base-learners (Hofner et al., 2015), we additionally used stability selection to compare the relative importance of covariates. Briefly, this process involved modeling subsamples of the data and measuring the frequency with which each covariate was included in the final models (Hofner et al., 2015). The methodology and results of this analysis are included in the Supporting information (Appendix S1).

| Synthesis
Both GAM and GAMLSS models took the following general form: For occupancy models (GAM), we modeled the occupancy probability of a given duck species in a segment g(π sea ducks ) as a function of all environmental covariates (Table 1), with g representing the logit link. Count models (GAMLSS) took two forms: the (conditional) mean count of sea ducks, g(μ sea ducks ), and the (conditional) overdispersion in sea duck counts, g(σ sea ducks ); g is the log link in both cases.
The same environmental covariates were included in both models (Table 1). f(·) indicates the penalized nonlinear deviations from the corresponding linear base-learners (e.g., f(time)) and were included for all non-categorical variables.
We included interaction terms between easting (xkm), northing (ykm), and day of year (time) to incorporate spatiotemporal effects.
The explicit intercept (int) was a necessary byproduct of our decomposition of base-learners (Hofner et al., 2011;Kneib et al., 2009). We included obs_window, our measure of survey effort (Table 1), as a covariate rather than an offset because small values in some segments impaired estimability.
Subsequent to their independent fitting, we consolidated occupancy and conditional count models into a single model (see Equation 6 in Zeileis, Kleiber, & Jackman, 2008) to generate unconditional, spatially-explicit estimates of sea duck abundance.

| Validation
Since additional test data were not available for our study area, we used a pseudo R 2 measure of the explained variation (Maloney et al., 2012;Nagelkerke, 1991) to evaluate the approximate explanatory power of our final models. We obtained the pseudo R 2 value by comparing the log likelihood values for our model-generated estimates of unconditional abundance for each species to log likelihood values obtained from null (intercept-only) models, giving us an estimate of the increase in explanatory power provided by our models.
g( ⋅ ) = int + covariate 1 + f(covariate 1 ) ⋯ + covariate n + f(covariate n ) F I G U R E 3 Marginal functional plots for stably selected covariates in the occupancy (probability of presence) and conditional abundance (mean and overdispersion) models for Common Eider (COEI), scoters (SCOT), and Long-tailed Duck (LTDU) in Nantucket Sound during three winters (2003-2004 to 2005-2006). Each plot illustrates the partial contribution of a covariate to the additive predictor (Y-axis), holding all other covariates at their mean. Within a model, plots share a Y-axis scale, enabling direct comparisons of effect sizes among covariates and species. Bivariate plots reflect the first (Y-axis) and second (X-axis) variables listed in the interaction; colors indicate the direction and magnitude of the partial contribution (blacks = negative, reds = positive; darker colors = larger effect). Northing by easting effects are given only at 31 December. For factor variables, only the general association (positive or negative) with the additive predictor is given. Covariate abbreviations correspond to Table 1 | 2355 SMITH eT al.

| RE SULTS
For each species or species group, we fitted independent models for occupancy and conditional count data (mean and overdispersion). Bootstrapped empirical risk suggested that occupancy models for all species converged to the maximum likelihood estimates (i.e., occupancy models failed to stop early; see Supporting information Appendix S2). Conversely, bootstrapped empirical risks prescribed early stopping for both the conditional mean and overdispersion parameter in all count models (see Supporting information Appendix S2). Final occupancy models and models for the conditional mean of count data included only a subset (12% to 38%) of the 48 base-learners initially specified for selection. Occupancy models generally contained more covariates and their interactions (8-10 of 23) than did count models (3-6 of 23), particularly among stably selected covariates and their interactions (Figure 3, see also Supporting information Appendix S3).

| Sea duck occupancy
Three covariates-grain size, sea surface temperature, and distance to land-were associated with probability of occurrence for all three sea duck species or species groups ( Figure 3). Standardized effects of each variable on the response can be compared among species and covariates within a model based on their range on the Y-axis.
For example, monthly sea surface temperature (SST m ) spanned a larger range of the Y-axis, and thus associated more strongly with eider occupancy, than did distance to land (d2land) (Figure 3). In contrast, monthly sea surface temperature (SST m ) associated much more strongly with occupancy of Long-tailed Duck than with eider or scoters ( Figure 3). Detailed comparisons of univariate, bivariate, and categorical effects for each species are included in the Supporting information (Appendix S3).
Spatiotemporal effects (i.e., occupancy associated with the xkmykm location of segments and the change over time within winter

| Sea duck conditional abundance and overdispersion
Spatial effects (ykm-xkm) were the dominant explanatory feature of conditional abundance estimates for scoters and Long-tailed Duck, but they were not selected in the eider model (Figure 3). In contrast with the corresponding occupancy model, scoter conditional abundance decreased with increasing sediment grain size (meanphi).
Additionally, the relationships between eider conditional abundance and dissolved organic material (cdom) and sea floor topography (SAR; Figure 3) were more complex than with eider occupancy. The conditional abundance of eider and scoter was also associated with relatively warm or cool sea surface temperatures (SST rel ; Figure 3).
Biophysical covariates associated with Long-tailed Duck conditional abundance exhibited general agreement with their counterpart in the occupancy models.
Spatially-explicit patterns of occupancy (cf. Figure 4, top row) did not necessarily reflect patterns of median conditional abundance ( Figure 5, top row). Some areas of Nantucket Sound exhibited mutually high conditional abundance and occupancy for a given species (e.g., eider in the southwest, scoter in the interior, and Long-tailed Duck in parts of the northeast). However, conditional abundance was low despite relatively high occupancy in some instances (e.g., eider in the northeast and Horseshoe Shoal, scoters in the northeast and southeast, and Long-tailed Duck along the northern margin). Conversely, other areas of Nantucket Sound exhibited lower occupancy but sea ducks, when present, were more abundant (e.g., eider along the eastern margin, and scoters and Long-tailed Duck in the southwest). As in occupancy models, sea ducks were relatively absent from the middle-western margin of Nantucket Sound (i.e., northeast of Martha's Vineyard; see Figure 3). In contrast to occupancy, which was less variable in areas of high abundance, areas of high conditional sea duck abundance typically also exhibited high relative variability over time ( Figure 5, bottom row).

| Expected sea duck abundance
Consolidated occupancy and conditional count models provided estimates of unconditional sea duck abundance in the study area over the survey period. Final models of expected sea duck abundance explained moderate amounts of variation in observed counts of eider, scoters, and Long-tailed Duck (pseudo R 2 = 0.31, 0.48, and 0.32, respectively). Conditional abundance ( Figure 5) strongly influenced the spatially-explicit patterns of expected abundance. Sea duck species exhibited relatively distinct patterns of abundance in Nantucket Sound. Eider were consistently most abundant in southwestern Nantucket Sound. They also were relatively abundant in the northeastern part of the study area but less consistently based on the relatively high MAD/median abundance over time ( Figure 6). Scoters were also most abundant, with occasional extremely large flocks, in southwestern Nantucket Sound. This was also the area of highest relative variation in scoter abundance; relatively high abundances of scoters also occurred in interior Nantucket Sound (Figure 6).

Long-tailed Ducks were consistently most abundant in northeastern
Nantucket Sound, as well as along its southern margin ( Figure 6). No species' highest abundances occurred in the permitted Nantucket Shoal area, although expected eider and scoters abundances were consistently elevated in some parts of the Shoal (west and southeast, respectively; Figure 6).
Summing the spatially-explicit estimates of unconditional sea duck abundance (i.e., Figure 6) provides an estimate of total abundance in a 1.5 km × 180 m transect through all segments in the study area. We calculated overall abundance by species throughout the study area by extrapolating these estimates across the full study area (Figure 7). Although absolute estimates differed between study years, patterns of abundance were similar across years, with scoter and long-tailed duck abundances highest early in the season, and eider abundance peaking in mid-winter. We also

| Temporal dynamics in wintering sea ducks
The MAD/median estimates (Figures 4-6, bottom rows; Supporting information Appendix S4) show that our spatially-explicit estimates of occupancy, abundance, and overdispersion invariably change over time, either explicitly via the selection of a within-or among-winter temporal effect (time and y2004/y2005, respectively) or implicitly via the selection of biophysical covariates that change within or among winters. The temporal dynamics of the wintering sea duck system in Nantucket Sound was one of its most striking attributes, and we illustrate these dynamics with an animation for scoter occupancy and abundance in the Supporting Information (Appendix S5).

| D ISCUSS I ON
We can be used to generate a variety of information on species abundance and distribution, including overall abundance estimates (this study), as well as estimates of prediction error, which can be generated via bootstrapping by repeated runs of the model (Hofner, Kneib, & Hothorn, 2016

| Environmental covariates that best explain sea duck distribution and abundance
The biophysical associations with sea duck occupancy derived from our models were relatively consistent among species, whereas their associations with sea duck conditional abundance were more species-specific. Distance to land, which was associated with both occupancy and abundance, tends to be positively associated with bathymetry and often has a strong influence on sea duck occupancy estimates (Flanders et al., 2015;Guillemette et al., 1993;Lewis, Esler, & Boyd, 2008;Winiarski et al., 2014). Sediment grain size can also affect prey availability for foraging sea ducks (Goudie & Ankney, 1988 (Logerwell & Hargreaves, 1996;Mannocci et al., 2017); thus, the smaller spatial scale of our analysis compared to previous studies may explain the apparent discrepancy between studies in the effect of chlorophyll a and NAO.
The unexplained variation in our models and the predominance of marginal spatiotemporal effects suggest that we likely missed important variable(s) relevant to the distribution of sea ducks in Nantucket Sound. This lack of explanatory power suggests a need for better biophysical proxies for the distribution of prey eaten by sea ducks or concurrent prey distribution information (Cervencl & Fernandez, 2012;Cervencl et al., 2014;Kaiser et al., 2006;Vaitkus & Bubinas, 2001;Žydelis, Esler, Kirk, & Boyd, 2009)

| The importance of spatial scale
The process whereby migratory animals such as sea ducks select a given area to inhabit during winter involves decisions at multiple spatial scales and the environmental attributes that determine this habitat selection often vary with spatial scale (Johnson, 1980;Johnson et al., 2006Johnson et al., , 2004. The majority of North American sea ducks migrate from high-latitude arctic and sub-arctic breeding areas to mid-latitude temperate wintering areas where they reside for most of the year (Bowman et al., 2015;Flanders et al., 2015;Silverman et al., 2013). At these large spatial scales, the distribution and abundance of sea ducks during winter may be affected by largescale ocean characteristics (Flint, 2013), climatic conditions (Zipkin et al., 2010), and static or persistent habitat features (e.g., bathymetry, substrate, current and frontal systems; Hyrenbach, Forney, & Dayton, 2000;Nur et al., 2011;Flanders et al., 2015). At regional and local scales, however, most species of sea ducks congregate in large flocks (e.g., up to tens of thousands of birds) at sites where prey are abundant and accessible (Flint, 2013;Loring et al., 2013) although the abundance and distribution of these prey, and thus predators, can be extremely ephemeral and dynamic (Cisneros, Smit, Laudien, & Schoeman, 2011;Hyrenbach et al., 2000).
Given that sea duck distribution and abundance is spatially and temporally dynamic, yet expected to be driven by biophysical covariates (Flanders et al., 2015;Oppel, Powell, & Dickson, 2009;Zipkin et al., 2010) that may differ in importance depending on spatial scale (Johnson, 1980;Johnson et al., 2006Johnson et al., , 2004, marine spatial planners must carefully consider the most appropriate information to use when deciding, for example, where to place marine protected areas or offshore wind energy developments to meet conservation and management goals. A larger-scale occupancy model developed by Flanders et al. (2015) suggested that eiders were relatively uniformly distributed across Nantucket Sound, whereas our higher resolution abundance models found that eiders were concentrated in southwestern and central, eastern areas within Nantucket Sound. While large-scale models (Flanders et al., 2015;Silverman et al., 2013) are useful to identify general geographic areas of importance to sea ducks, our models provide more detailed estimates of sea duck distribution and abundance within a specific area of interest. In terms of marine spatial planning, large-scale models can be used to inform the siting of lease areas or protected areas, followed by detailed modeling approach such as ours to select sites within these larger blocks that can be zoned for specific levels of development or protection.

ACK N OWLED G EM ENTS
We