OCS Study

BOEM 2018-010

Modeling At-Sea Density of Marine Birds 

to Support Atlantic Marine Renewable 

Energy Planning

Final Report 

US Department of the Interior

Bureau of Ocean Energy Management

Office of Renewable Energy Programs

OCS Study 

BOEM 2018-010

Modeling At-Sea Density of Marine Birds 

to Support Atlantic Marine Renewable 

Energy Planning

Final Report

Authors

Arliss J. Winship, Brian P. Kinlan, Timothy P. White, Jeffery B. Leirness, John Christensen

Prepared under NCCOS IAA MOA-2013-046-8696, BOEM OCS Study 2018-010, and NCCOS 

BOEM IAA M13PG00005

by

U.S. Department of Commerce

National Oceanic and Atmospheric Administration 

National Ocean Service

National Centers for Coastal Ocean Science

Marine Spatial Ecology Division

Biogeography Branch

1305 East-West Hwy, SSMC-4, N/SCI-1 

Silver Spring, MD 20910

Published by

U.S. Department of the Interior

Bureau of Ocean Energy Management

Office of Renewable Energy Programs

June 2018
ii

DISCLAIMER

This study was funded, in part, by the US Department of the Interior, Bureau of Ocean Energy 

Management (BOEM), Environmental Studies Program, Washington, DC, through Inter-Agency 

Agreement Number M13PG00005 with the US Department of Commerce, National Oceanic and 

Atmospheric Administration, National Ocean Service, National Centers for Coastal Ocean 

Science, Silver Spring, MD. This report has been technically reviewed by BOEM and it has been 

approved for publication. The views and conclusions contained in this document are those of the 

authors and should not be interpreted as representing the opinions or policies of the US 

Government, nor does mention of trade names or commercial products constitute endorsement 

or recommendation for use.

REPORT AVAILABILITY

To download a PDF file of this Environmental Studies Program report, go to the US Department 

of the Interior, Bureau of Ocean Energy Management, Environmental Studies Program 

Information System website and search on OCS Study BOEM 2018-010. You may request the 

report from the BOEM Office of Renewable Energy Programs. The contact information is:

U.S. Department of the Interior

Bureau of Ocean Energy Management

Office of Renewable Energy Programs

45600 Woodland Road, VAM-OREP

Sterling, VA 20166

CITATION

A.J. Winship, B.P. Kinlan, T.P. White, J.B. Leirness, and J. Christensen. 2018. Modeling At-Sea 

Density of Marine Birds to Support Atlantic Marine Renewable Energy Planning: Final Report. 

U.S. Department of the Interior, Bureau of Ocean Energy Management, Office of Renewable 

Energy Programs, Sterling, VA. OCS Study BOEM 2018-010. x+67 pp. 

ABOUT THE COVER

Cover photo (Herring Gull) courtesy of David Pereksta (BOEM). Used with permission.
iii

ACKNOWLEDGEMENTS

We thank the many scientists who collected and contributed the survey data analyzed in this 

study. We are grateful to Allan O’Connell, Mark Wimer, Allison Sussman, Tim Jones, Kaycee 

Coleman, Kyle Dettloff, and Robert Fowler for processing and providing the data in the 

Northwest Atlantic Seabird Catalog (formerly Avian Compendium) database and to Carina 

Gjerdrum (Canadian Wildlife Service, Environment and Climate Change Canada) for providing 

the Eastern Canada Seabirds at Sea survey data. We thank Peter Cornillon, Michael Coyne, 

Peter Miller, and Matthew Poti for providing and helping with the processing of environmental 

predictor data. We thank Robert Rankin for initial development of the modeling framework and 

computer code, and we thank Zhifa Liu for execution of the models on high performance 

computing systems. We are grateful to Todd Callaghan, Andrew Gilbert, Carina Gjerdrum, Wing 

Goodale, Holly Goyert, Christopher Haney, Kevin Powers, Caleb Spiegel, and David Wiley for 

helpful feedback on the maps that improved this report. We thank David Bigger and Mary 

Boatman for comments and suggestions that improved this report. This project was funded by 

the Bureau of Ocean Energy Management (BOEM) through Intra-agency Agreement 

M11PG00059 with the United States Geological Survey (USGS) and Inter-agency Agreement 

M13PG00005 with the U.S. Department of Commerce, National Oceanic and Atmospheric 

Administration, National Ocean Service, National Centers for Coastal Ocean Science (NCCOS), 

and by the United States Geological Survey through Inter-agency Agreement G13PG00008 with 

NCCOS. Arliss Winship, Timothy White, and Jeffery Leirness were supported by NOAA 

Contracts No. DG133C07NC0616 and EA-133C-14-NC-1384 with CSS, Inc. (formerly CSSDynamac). 
iv

Executive Summary 

Marine birds have the potential to be affected by human activities in the ocean environment such as 

offshore wind energy development. This report describes a project that developed maps of the spatial 

distributions of marine bird species in U.S. Atlantic Outer Continental Shelf (OCS) waters that can be 

used to aid ocean planning in the region and guide future data collection efforts.

Sighting survey data from over three decades contained in the ‘Northwest Atlantic Seabird Catalog’ 

database, along with Eastern Canada Seabirds at Sea data from Canadian Wildlife Service, Environment 

and Climate Change Canada, were analyzed to derive seasonal maps of the spatial distributions of 47 

marine bird species in U.S. Atlantic OCS and adjacent waters from Florida to Maine.

Spatial predictive modeling was applied to the survey data to account for spatial and temporal 

heterogeneity in survey effort, platform, and protocol. An ensemble machine-learning technique, 

component-wise boosting of hierarchical zero-inflated count models, was used to relate the relative 

density of each species to multiple spatial and temporal predictor variables while accounting for survey 

heterogeneity and the aggregated nature of sightings. Dynamic spatial environmental predictor variables 

were formulated as long-term climatologies. The modeling technique allowed for complex non-linear 

relationships between response and predictor variables and interacting effects among predictors. 

Bootstrapping was used to derive estimates of the uncertainty in model predictions. 

Model predictions are presented as seasonal maps of the relative density of each study species throughout 

the study area. The maps were reviewed by experts with experience and knowledge of marine birds in the 

study area and their comments were incorporated in this report. The maps indicate where species are 

likely to be more or less abundant. The analysis was not designed to estimate the actual number of 

individuals/density of a given species that would be expected in any location, so the maps should not be 

interpreted that way. Also, the maps represent the spatial distributions of birds averaged over time (e.g., 

across days within a season and across years for a given season). The analysis was not designed to 

provide predictions of the density of birds that would be expected in a specific location at a specific date 

or time, so the maps should also not be interpreted that way.

Two indications of the uncertainty associated with the model predictions are provided. First, a hatched 

overlay is included on the maps of predicted relative density to indicate areas with no survey effort. 

Model predictions in areas with no survey effort should be interpreted with extreme caution. Predictions 

in these areas were often questionable or unrealistic, so we recommend additional field surveys in these 

areas to validate the model predictions. Second, estimates of the precision of model predictions are 

presented as maps of the coefficient of variation (CV) of predicted relative density. Less precise 

predictions (i.e., higher CV) should be interpreted with more caution. The maps of predicted relative 

density should always be considered in conjunction with these two indications of uncertainty. 

The relative importance of different predictor variables is also presented, indicating which variables most 

influenced the predicted distributions for each species in each season. While the primary objective of this 

study was not to determine the ecological drivers and mechanisms behind the spatial distributions of 

marine bird species in the study area, our model results may provide useful hypotheses for future studies 

aimed more at ecological inference.
v 

Modeling At-Sea Density of Marine Birds to Support Atlantic Marine 

Renewable Energy Planning: Final Report

Contents

Executive Summary ........................................................................................... iv

List of Figures .................................................................................................... vii

List of Tables .................................................................................................... viii

List of Appendices ............................................................................................. ix

Abbreviations and Acronyms ............................................................................ x

1. Introduction ..................................................................................................... 1

2. Methods ............................................................................................................ 2

2.1 Overview................................................................................................................ 2

2.2 Survey data ........................................................................................................... 2

2.3 Species modeled ................................................................................................... 4

2.4 Predictor variables................................................................................................. 4

2.5 Statistical modeling framework .............................................................................. 6

2.5.1 Likelihoods and model components ........................................................................... 6

2.5.2 Effort offset ................................................................................................................ 6

2.5.3 Base-learners............................................................................................................. 7

2.5.4 Stochastic gradient boosting ...................................................................................... 7

2.5.5 Boosting offsets ......................................................................................................... 7

2.5.6 Tuning of learning rate and number of boosting iterations.......................................... 7

2.5.7 Model selection and performance .............................................................................. 8

2.5.8 Spatial prediction ....................................................................................................... 9

2.5.9 Variable importance ................................................................................................... 9

2.5.10 Uncertainty............................................................................................................... 9

2.5.11 Implementation .......................................................................................................10

2.6 Map display ......................................................................................................... 10

2.7 Expert review of maps ......................................................................................... 11
vi

2.8 Predicted proportional relative abundance by BOEM wind energy planning/lease 

area ........................................................................................................................... 11

3. Results ........................................................................................................... 11

3.1 Model selection and performance........................................................................ 11

3.2 Predicted spatial distributions .............................................................................. 12

3.2.1 Spatial and seasonal patterns ...................................................................................12

3.2.2 Uncertainty................................................................................................................14

3.2.3 Problematic predictions.............................................................................................14

3.3 Predictor variable relative importance.................................................................. 16

4. Discussion ..................................................................................................... 17

4.1 Interpretation of maps.......................................................................................... 17

4.2 Data limitations and information gaps.................................................................. 18

4.3 Combining seasonal maps .................................................................................. 19

4.4 Species identification........................................................................................... 20

4.5 Comparison to other marine bird distribution models........................................... 20

4.6 Comparison to tracking studies ........................................................................... 22

4.7 Conclusion........................................................................................................... 23

5. Literature Cited .............................................................................................. 24

Glossary ............................................................................................................. 65
vii

List of Figures

Figure 1. Study area with BOEM Wind Energy Lease and Planning Areas overlaid (approximate 

boundaries current as of 2017-04-10)................................................................................53

Figure 2. Area surveyed by year. ..............................................................................................54

Figure 3. Total area surveyed each season within the study area, binned into 10 x 10 km cells. 

White areas represent no survey effort. See Appendix A for maps of individual survey 

effort..................................................................................................................................55

Figure 4. Schematic overview of statistical modeling process. See Section 2 Methods for details.

..........................................................................................................................................56

Figure 5. Example maps of predicted relative density and its coefficient of variation (CV) for one 

species (Atlantic Puffin) and season (spring).....................................................................57

Figure 6a. Relative importance of predictor variables for p (red) and μ (blue) components of best 

spring (March-May) models. Areas of circles are proportional to importance.....................58

Figure 6b. Relative importance of predictor variables for p (red) and μ (blue) components of best 

summer (June-August) models. Areas of circles are proportional to importance................59

Figure 6c. Relative importance of predictor variables for p (red) and μ (blue) components of best 

fall (September-November) models. Areas of circles are proportional to importance.........60

Figure 6d. Relative importance of predictor variables for p (red) and μ (blue) components of best 

winter (December-February) models. Areas of circles are proportional to importance. ......61

Figure 7. Example annual maps of predicted relative density and its coefficient of variation (CV) 

for one species (Atlantic Puffin). ........................................................................................62

Figure 8. Example map of predicted relative density for Coastal Waterfowl species group during 

winter. Coastal Waterfowl included Black Scoter, Common Eider, Common Loon, Longtailed Duck, Red-throated Loon, Surf Scoter, and White-winged Scoter. ...........................63

Figure 9. Maps of predicted relative density for phalarope multi-species models for spring and 

fall. These models were fit to pooled sighting data for Red Phalarope, Red-necked 

Phalarope, and unidentified phalaropes. ...........................................................................64
viii

List of Tables

Table 1. Survey datasets, sample sizes, and surveyed area analyzed. All datasets (and their ID 

codes) were extracted from the USFWS Northwest Atlantic Seabird Catalog except ECSAS 

which was provided by CWS-ECCC. Continuous data were segmented as described in 

Section 2.2. A description of each individual survey is presented in Appendix A. Area 

surveyed does not necessarily represent unique area surveyed; some surveys covered 

similar areas over multiple days, months, or years. Seasons were defined as spring 

(March-May), summer (June-August), fall (September-November), and winter (DecemberFebruary). Sample sizes and surveyed area analyzed for corrected MassCEC datasets 

(see Section 2.2) are shown in parentheses......................................................................32

Table 2. List of species-season combinations modeled with sample sizes and counts. Only 

combinations with ≥50 sightings were modeled. Seasons were defined as spring (MarchMay), summer (June-August), fall (September-November), and winter (DecemberFebruary)...........................................................................................................................36

Table 3. Predictor variables used in models with native resolutions, year ranges, and sources.38

Table 4a. Pairwise Spearman rank correlation coefficients for spatial predictor variables (spring 

– March-May). High correlations are highlighted in yellow (>0.7), orange (>0.8), and red 

(>0.9).................................................................................................................................42

Table 4b. Pairwise Spearman rank correlation coefficients for spatial predictor variables 

(summer – June-August). High correlations are highlighted in yellow (>0.7), orange (>0.8), 

and red (>0.9)....................................................................................................................43

Table 4c. Pairwise Spearman rank correlation coefficients for spatial predictor variables (fall – 

September-November). High correlations are highlighted in yellow (>0.7), orange (>0.8), 

and red (>0.9)....................................................................................................................44

Table 4d. Pairwise Spearman rank correlation coefficients for spatial predictor variables (winter 

– December-February). High correlations are highlighted in yellow (>0.7), orange (>0.8), 

and red (>0.9)....................................................................................................................45

Table 5. Base-learners employed in the boosted generalized additive modeling framework. 

Base-learner names are from the ‘mboost’ package for R (Hothorn et al. 2015; R Core 

Team. 2017)......................................................................................................................46

Table 6. Best models with model performance metrics. ............................................................47
ix

List of Appendices

Appendix A. Description of individual survey datasets analyzed

Appendix B. Predictor variable plots and grids

Appendix C. Maps of predicted relative density and their coefficient of variation (CV) for each 

species and season modeled

Appendix D. Tables of predicted relative abundance as a proportion of total relative abundance 

in the study area by BOEM wind energy planning/lease area and season for each species 

modeled
x 

Abbreviations and Acronyms

AMO Atlantic Multidecadal Oscillation

AUC area under the ROC curve

BOEM Bureau of Ocean Energy Management

CV coefficient of variation

CWS-ECCC Canadian Wildlife Service, Environment and Climate 

Change Canada

ECSAS Eastern Canada Seabirds at Sea

ID identity

MEI Multivariate El Niño-Southern Oscillation Index

NAO North Atlantic Oscillation

NOAA National Oceanic and Atmospheric Administration

OCS Outer Continental Shelf

PDE percent deviance explained

RISAMP Rhode Island Ocean Special Area Management Plan

RMSE root mean square error

ROC receiver operating characteristic

SD standard deviation

SE standard error

SST sea surface temperature

TNI Trans-Niño Index

USGS United States Geological Survey

USFWS United States Fish and Wildlife Service

ZINB zero-inflated negative binomial

ZIP zero-inflated Poisson
1 

1. Introduction

Marine birds spend much of their time in coastal waters and on the open ocean. As a result, these species 

have the potential to be affected by human activities in the marine environment such as offshore wind 

energy development. A prerequisite for quantifying that potential is knowledge of the spatial distributions 

of marine birds at sea. This report describes a project aimed at producing maps of the spatial distributions 

of marine bird species in U.S. Atlantic Outer Continental Shelf (OCS) waters (Fig. 1) that can be used to 

inform marine spatial planning in the region.

Some of the best available information about the at-sea distributions of marine birds comes from visual 

sighting and photographic surveys conducted aboard boats and aircraft. For U.S. Atlantic OCS waters 

many data from past sighting surveys have been compiled in the ‘Northwest Atlantic Seabird Catalog’, 

hereinafter referred to as the Catalog. The Catalog was originally developed by the U.S. Geological 

Survey (USGS) Patuxent Wildlife Research Center as the ‘Compendium of Avian Occurrence 

Information for the Continental Shelf waters along the Atlantic Coast of the U.S.’ (O’Connell et al. 2009).

The Catalog is currently maintained by the U.S. Fish and Wildlife Service (USFWS). The project 

described here analyzed sighting data from the Catalog along with sighting data from Eastern Canada 

Seabird at Sea surveys (ECSAS) by Canadian Wildlife Service, Environment and Climate Change Canada 

(CWS-ECCC) to derive maps of the spatial distributions of 47 marine bird species in U.S. Atlantic OCS 

and adjacent waters from Florida to Maine.

The survey data analyzed represent numerous surveys spanning nearly four decades. Survey coverage and 

intensity was highly variable geographically and temporally. Furthermore, a range of survey platforms, 

observers, and protocols were used. This heterogeneity complicates the quantification of the at-sea 

distribution of marine birds from these data, and biases simple data summaries. To deal with this 

heterogeneity, the project described here employed spatial predictive modeling. An ensemble machinelearning technique was used to model counts of each species as a function of multiple predictor variables 

while accounting for heterogeneous survey effort. The fitted models were then used to predict the spatial 

distribution of relative density of each species throughout the study area. 

The distributions of marine birds at sea are a result of interactions between their behavior (e.g., foraging) 

and the environment. Atmospheric and oceanographic features and processes across a range of spatial and 

temporal scales influence the environmental conditions and prey availability experienced by marine birds, 

and thus ultimately determine their at-sea distributions. The spatial predictive modeling framework 

employed here relied on a wide suite of spatial and temporal environmental predictor variables to explain 

and predict the distributions of marine birds. In particular, static environmental variables (e.g., 

bathymetry) and long-term climatologies of dynamic environmental variables (e.g., sea surface 

temperature) were considered to explain spatial patterns of relative density. 

The project described here was designed to provide broad-scale spatial information that can be used to 

guide future data collection efforts and aid marine spatial planning in the region. It is important to note 

that the results presented in this report represent the spatial distributions of birds averaged over time. The 

project was not designed to provide precise predictions of the actual number of individuals of a given 

species that would be expected in a specific location at a specific time. The project was also not designed 

to determine the ecological drivers of marine bird distributions, although the results provide related 

hypotheses for future research.
2 

2. Methods

2.1 Overview

A statistical modeling framework was used to relate bird sighting data from historical surveys to a range 

of temporal and spatial environmental predictor variables. The estimated relationships between the 

counts of birds and the predictor variables were then used to predict the spatial distributions of birds 

across the entire study area. Separate models were developed for each combination of species and season 

for which there were sufficient data. Seasons reflected major transitions in environmental conditions in 

the study region: spring (1 March – 31 May), summer (1 June – 31 August), fall (1 September – 30 

November), and winter (1 December – 28/29 February). Note that these seasons do not necessarily align 

with the timing of migration for individual species. As a result, summer (June-August) and winter 

(December-February) models may at least partially represent distributions during migration (e.g., 

Northern Gannet in summer), and spring (March-May) and fall (September-November) models may not 

be entirely representative of migration per se.

2.2 Survey data 

Ninety-two survey datasets were analyzed (Table 1, Appendix A). The majority of the data were obtained 

from the ‘Northwest Atlantic Seabird Catalog’ maintained by USFWS through Intra-agency Agreement 

M14PG00014 with the US Bureau of Ocean Energy Management (BOEM) (Table 1, Appendix A). The 

Catalog was formerly known as the ‘Compendium of Avian Occurrence Information for the Continental 

Shelf waters along the Atlantic Coast of the U.S.’ or ‘Avian Compendium’, and was originally created 

and maintained by USGS Patuxent Wildlife Research Center through Intra-agency Agreement 

M11PG00059 with BOEM. The Catalog is a relational database containing a large number of datasets on 

the occurrence and distribution of marine birds along the Atlantic coast of the U.S. (O’Connell et al. 

2009). These datasets were collected by a range of entities including government agencies, nongovernmental organizations, academic researchers, and other individuals.

The quality of the datasets in the Catalog varies from systematic scientific surveys with rigorous sampling 

protocols to opportunistic sightings. We used only science-quality, geographically-referenced datasets. 

Opportunistic sighting data (e.g., eBird) often do not have associated quantitative measures of survey 

effort and often lack explicit information about observed absences. Without data on survey effort it is not 

possible to determine whether large numbers of sightings reflect high density or simply reflect more 

survey effort. Without data on observed absences it is not possible to calibrate the estimated probability of 

occurrence or density. For these reasons opportunistic sighting data were excluded from our analysis. We 

also excluded data from survey effort where chumming was used to attract birds. We aimed to estimate 

the typical spatial distributions of birds, therefore we did not use data from survey effort where the local 

density of birds was intentionally increased temporarily by the survey itself. 

Additional survey data were obtained from CWS-ECCC. CWS-ECCC provided ECSAS data from 

surveys in eastern Canadian and US Atlantic waters. An accidental duplication of data occurred because 

one of the Catalog datasets (NewEnglandSeamount06) was also included in the ECSAS data; however, 

the duplicated effort represented only 0.02% of the total surveyed area analyzed so the duplication likely 

had a negligible impact on our analysis.

The datasets analyzed spanned 1978-2016 with most survey effort occurring from 1978-1988 and from 

2002 onward, especially between 2010 and 2014 (Fig. 2). A total of 181,140 km2 was surveyed with more 
3 

survey effort during summer and fall than during winter and spring (Table 1). Most survey effort was in 

coastal waters and over the continental shelf with much less effort beyond the shelf break (Fig. 3). Survey 

effort was greatest offshore during the summer and least during the winter. Year-round concentrations of 

survey effort occurred near the mouths of Delaware Bay and Chesapeake Bay, in Nantucket Sound, and to 

the southwest of Nantucket Sound. The datasets with the largest combined sample size and widest 

geographic coverage were collected by Manomet Bird Observatory in coordination with the US National 

Oceanic and Atmospheric Administration (NOAA) and other cruises between 1978 and 1988 (datasets 

CSAP and NOAAMBO7880). More recent aerial surveys by USFWS and boat surveys by NOAA as part 

of the Atlantic Marine Assessment Program for Protected Species (AMAPPS) also covered U.S. Atlantic 

coastal (aerial) and offshore (boat) waters. Other multi-year survey efforts covered large sections of the 

coast including NOAA ecosystem monitoring cruises from North Carolina to the Gulf of Maine (datasets 

EcoMon), pelagic surveys off Georgia, South Carolina, and Florida (dataset GeorgiaPelagic), and the 

Canadian ECSAS surveys. The remaining datasets are more localized, often from New England and the 

Gulf of Maine or the mid-Atlantic, but sometimes they had large sample sizes over multiple years (e.g., 

CapeWind, DOEBRI, HerringAcoustic, MassAudNanAerial, and MassCEC). 

The survey data analyzed were counts of birds mainly from aerial and boat-based visual surveys at sea. 

Three datasets (DOEBRIAerial) represented high-resolution digital video aerial surveys. The original data 

took the form of species-specific counts along strip transects. The width of the strip transect was reported 

for most surveys, but for 5 boat-based surveys where it was not, a 300-m wide strip transect was assumed, 

which is the standard for boat-based seabird surveys (Tasker et al. 1984). 

Counts were sometimes recorded continuously, and other times had been binned into discrete transect 

segments (e.g., 15 minute duration). Binned data were only from boat surveys. If the distance travelled 

during each transect segment was not available then it was estimated from the duration of the recording 

period and the vessel speed. The geographic midpoint for each segment was also estimated from the 

duration of the recording period and the vessel speed and direction. For consistency with the binned data 

the continuously recorded data were discretized into transect segments of approximately 4 km in length, 

and counts for each species were summed within each transect segment. Because the length of any given 

transect was not necessarily divisible by 4 km, any remainder distance was either treated as its own 

segment (if it was >2 km) or was added to another segment (if it was <2 km). The placement of the 

resulting short or long segment along a transect was randomized to avoid it always occurring at the end of 

a transect. Any transects that were originally <6 km in length were treated as single transect segments. 

The distance travelled during each transect segment was calculated assuming straight-line travel between 

recorded vessel locations, and the geographic midpoint of each transect segment was calculated as the 

point along the trackline where half the distance had been travelled. Transect segments <1 km in length 

were excluded from analysis.

An error in the segmentation of the MassCEC2011-2012, MassCEC2013, and MassCEC2014 datasets 

was discovered after analysis. The error resulted in the geographic misplacement of sightings. The error 

was corrected, these datasets were re-processed, and models were re-run for species for which these 

datasets represented a large proportion of the total count or whose sighting data and model predictions 

were noticeably affected by the error (Table 2).

The survey datasets analyzed varied in terms of the type of survey platform used, observer identity and 

expertise, species focus, and environmental conditions. The type of survey platform used (i.e., boat or 

plane), characteristics of the survey platform (e.g., observation height), detection method (visual or 
4 

photographic), and observer expertise all influence the probability that individual birds will be detected 

and correctly identified to the species level. For some surveys the primary target was marine birds, while 

for others birds were surveyed along with other components of the ecosystem. Even among bird-focused 

surveys there was variability in the species of primary focus (e.g., sea ducks). The level of focus on 

marine birds and specific species also influences the probability that individual birds will be detected and 

correctly identified. Sighting conditions (e.g., sea state and visibility) can also greatly affect sighting 

rates. Ideally these factors would be explicitly accounted for in any model of density. Unfortunately, data 

on many of these factors were not readily available for all of the survey datasets analyzed here. 

Nevertheless, the analytical framework used attempted to account for this heterogeneity across survey 

datasets through effects of survey platform type, transect, and survey (see Section 2.4).

It is important to recognize that because the influence of these factors on sighting rates was not modeled 

explicitly, the spatial distributions of density presented here should be interpreted as relative predictions, 

not predictions of actual densities.

2.3 Species modeled 

Species with ≥100 transect segments with sightings of ≥1 individual in at least one season were selected 

for analysis. For each of those species, all seasons with ≥50 transect segments with sightings of ≥1 

individual were modeled. By these criteria 47 species and 140 species-season combinations were 

analyzed representing sightings of 2,687,745 individual birds (Table 2). These species represented 13 

families. The five species-season combinations with the greatest numbers of transect segments with

sightings were Wilson’s Storm-Petrel in summer, Herring Gull in fall, Northern Gannet in winter, and 

Great Shearwater in fall and summer.

2.4 Predictor variables

A wide range of predictor variables were used to model variation in the number of birds counted per 

transect segment and to predict the spatial distributions of birds throughout the study area (Table 3, 

Appendix B). Predictor variables fell into one of six categories: survey, temporal, geographic, 

bathymetric, oceanographic, and atmospheric. 

Survey predictor variables were designed to account for variation in counts arising from heterogeneity in 

the type of survey platform, characteristics of the survey platform (e.g., observation height), observer 

identity and expertise, species focus, and sighting conditions. These factors influence the probability that 

individual birds will be detected (e.g., Heinänen et al. 2017) and correctly identified to the species level. 

Of these factors, only the type of survey platform (boat-visual, aerial-visual, or aerial-photographic) was 

consistently available for all datasets, and thus was directly usable as a predictor variable. We attempted 

to account for the effects of the remaining factors through two random-effect predictor variables 

representing survey identity (ID) and transect ID, respectively. The exact definition of transect ID differed 

somewhat between datasets, but unique transect IDs generally represented pre-defined survey transects or 

individual days of effort. 

Temporal predictor variables were designed to account for variation in counts over time. Day of the year 

was used to account for changes in the numbers of birds in the study area over time within a season, for 

example arising from migratory movements in and out of the study area. Year was used to account for 

changes in the number of birds in the study area across years, for example arising from changes in 

population abundance or distributional shifts. Effects of day of the year and year were modeled as smooth 
5 

continuous changes over time. Four climate indices (Table 3) were also included as temporal predictor 

variables to account for variation in counts across years and across months within seasons arising from 

linkages between the environment and population abundance and distribution. For each climate index two 

values were included as predictor variables: the value for the month and year of a given transect segment 

and the value for the same month one year previous.

Geographic predictor variables were designed to account for variation in counts arising from spatial 

location per se. Projected longitude and latitude were included as predictor variables and their effects 

were modeled two ways. The first longitude-latitude predictor term allowed for smooth changes in 

numbers across the study area arising from spatial factors not captured by the other predictor variables 

(e.g., colonization history). The second longitude-latitude predictor term was formulated using radial 

basis functions with the intent of capturing some of the spatial autocorrelation in the data after accounting 

for the effects of other predictor variables. Distance to land was also included as a geographic predictor 

variable. 

Bathymetric predictor variables were designed to account for variation in counts arising from the direct 

and indirect effects of bathymetry on bird distributions. A depth predictor variable was developed by 

combining information from six different bathymetric datasets (Table 3). Other bathymetric variables 

were derived from depth including slope, slope of slope, and planform and profile curvature. 

Oceanographic and atmospheric predictor variables were designed to account for variation in counts 

arising from the direct and indirect effects of the physical state and dynamics of the ocean and air above 

the ocean as well as biological productivity. Seventeen oceanographic and atmospheric predictor 

variables were developed from a range of data sources (Table 3). Remote sensing data were used to 

characterize sea surface height, temperature, chlorophyll-a, turbidity, and wind stress. Other variables 

were derived from the remotely sensed variables including sea surface height and temperature variability, 

probabilities of cyclonic and anticyclonic eddy rings, probability of sea surface temperature fronts, wind 

divergence, and an index of upwelling. Estimates from a data-assimilating ocean dynamics model were 

used to characterize water currents, and divergence and vorticity were derived from current velocities. 

All of the oceanographic and atmospheric variables that we considered are dynamic. We formulated these 

predictor variables to characterize long-term spatial patterns in average values and variability. Long data 

time series ranging in length from 10-25 years were used (Table 3). To characterize average values, 

monthly mean climatologies across years were developed and then integrated to create seasonal 

climatologies. To characterize variability, standard deviations or probabilities (frequencies) were 

calculated from the native temporal resolution of the corresponding predictor variables. 

Geographic, bathymetric, oceanographic, and atmospheric predictor variables were spatially explicit. 

Each variable was calculated on a standard study grid with a spatial resolution of 2 km and an oblique 

Mercator projected coordinate system (origin = 35°N 75°W; azimuth = 40°; scale = 0.9996; NAD83 earth 

datum). When the native spatial resolution of a predictor variable was finer than that of the study grid, 

predictor values were averaged within study grid cells. When the native spatial resolution of a predictor 

variable was similar to or coarser than that of the study grid, bilinear interpolation was used to derive 

predictor values at the center of study grid cells. Each survey transect segment was matched to the 

predictor variable values from the study grid cell that contained the midpoint of that segment. 

Some of the spatially explicit predictor variables were highly correlated with each other (Table 4). 

Predictor variables were chosen to avoid correlations >0.9, although summer and fall each had three 
6

pairwise correlations >0.9. Four of the correlations >0.9 involved spatial coordinate variables that were a 

key structural component of our model. The other two correlations >0.9 were between chlorophyll-a and 

turbidity and between the mean and standard deviation of sea surface temperature during the summer. 

Because of the high correlations between some predictor variables, inferences regarding relative variable 

importance should be made with caution. The accuracy of predictions should be less affected by 

collinearity among predictor variables.

2.5 Statistical modeling framework

A boosted generalized additive modeling framework (Bühlmann and Hothorn 2007; Hofner et al. 2012) 

was used to estimate relationships between the numbers of birds counted per transect segment and the 

predictor variables (Fig. 4). Those relationships were then used to predict the relative density of each 

species throughout the study area in each season. Our main objective was to provide accurate predictions 

so we chose a modeling framework that allowed for flexible relationships and multiple interactions 

between predictor variables while accounting for sampling heterogeneity between and within datasets.

2.5.1 Likelihoods and model components

The number of individuals of a given species counted per transect segment was modeled using zeroinflated Poisson (ZIP; Eq. 1) and zero-inflated negative binomial (ZINB; Eq. 2) likelihoods to account for

the overdispersed nature of wildlife count data. Each component/parameter of the likelihood was modeled 

as a separate function of the predictor variables (Schmid et al. 2008; Mayr et al. 2012). For the ZIP

likelihood, the two model components were the probability of an ‘extra’ zero (p) and the mean of the 

Poisson distribution (μ):

[1] L(p, μ; y) = ∏ [p + (1 − p)e

−μ

]

Iyi=0

[(1 − p)

μ

yie

−μ

yi!

]

Iyi>0

n

i=1

The same components were modeled for the ZINB likelihood in addition to the dispersion parameter of 

the negative binomial distribution (θ):

[2] L(p, μ, θ; y) = ∏ [p + (1 − p) (

θ

θ+μ

)

θ

]

Iyi=0

[(1 − p)

Γ(yi+θ)

yi!Γ(θ)

(

θ

θ+μ

)

θ

(

μ

θ+μ

)

yi

]

Iyi>0

n

i=1

The probability of an extra zero (p) was modeled on the logit scale while the mean of the 

Poisson/negative binomial distribution (μ) and the dispersion parameter of the negative binomial 

distribution (θ) were modeled on the log scale.

In Eqs 1 and 2, y represents the vector of yi which represent the total count for transect segment i, n

represents the total number of segments, and Iyi=0 and Iyi>0 are indicators of whether yi is equal to or 

greater than zero, respectively (I=1 when the condition is true and I=0 when the condition is false).

2.5.2 Effort offset

To account for variation in the lengths and widths of strip transect segments an effort offset was included 

in the model. Specifically, the area surveyed per transect segment was used as an effort offset in the μ

model component, enforcing a proportional relationship between the area surveyed and the expected 

number of individuals counted on a transect segment conditional on that segment not being an ‘extra’ zero 
count. Model predictions correspond to counts per 1 km2

 of survey effort and therefore are referred to as 

density. 

2.5.3 Base-learners

Within the boosting framework, each model component was essentially modeled as a function of an 

ensemble of ‘base-learners’. Each base-learner represented a specific functional relationship between a 

model component and one or more predictor variables. We utilized a suite of base-learners each 

representing different predictor variables, and different sets of base-learners were employed for different 

model components (Table 5). 

All spatially explicit predictor variables except longitude and latitude were included together in a single 

tree base-learner. The trees for that learner had a maximum depth of 5, which allowed for interacting 

effects among the spatially explicit predictor variables. Longitude and latitude appeared in two base 

learners, and those variables always entered the model as a pair. The remaining survey and temporal 

predictor variables entered the model individually, either through their own base-learners or in the case of 

climate indices one at a time through a tree base-learner with a maximum depth of 1. Our model structure 

did not allow for interactions between temporal and spatial predictor variables.

2.5.4 Stochastic gradient boosting

Stochastic gradient boosting was used to fit models whereby a sub-sample of the data (4/5 during tuning, 

see Section 2.5.6, or 2/3 otherwise) was fitted in each boosting iteration (Friedman 2002). Rather than resampling the data for each boosting iteration, a set of 25 or 50 random sub-samples was created before 

boosting, and one sub-sample was randomly drawn from this set for each boosting iteration. Mean square 

error was used to select the base-learner that gave the best fit to the gradient in each boosting iteration. 

2.5.5 Boosting offsets

Model component estimates were initialized (‘offset’ in boosting terminology; Hofner et al. 2012) using 

one of two methods. The first method was to fit an intercepts-only ZIP or ZINB model with an effort 

offset to the data (‘pscl’ package for R; Zeileis et al. 2008, Jackman 2015), then use the estimated 

intercepts for each component of that model as the boosting offsets. The second method was to fit a 

simple zero-inflated generalized additive model (Poisson or negative binomial) with an effort offset to the 

data (‘mgcv’ and ‘pscl’ packages for R; Wood 2006, Zeileis et al. 2008, Jackman 2015) and then use the 

estimated intercepts for each component of that model as the boosting offsets. The model for the second 

method incorporated effects of projected longitude and latitude through a bivariate thin-plate regression 

spline with 3 or 4 degrees of freedom. 

2.5.6 Tuning of learning rate and number of boosting iterations

A stratified (by transect ID) 5-fold cross-validation approach was used to determine the number of 

boosting iterations and component-wise learning rates that resulted in the best predictive performance. 

Three different learning rates were considered for each model component: 0.001, 0.01, and 0.1. For each 

unique combination of component-wise learning rates (e.g., 32 = 9 for 2 model components) and crossvalidation fold, the model was fit to the in-bag data (i.e., 4/5 of the data) and the number of boosting 

iterations at which the negative log-likelihood of the out-of-bag data (i.e., 1/5 of the data) was minimized 

was determined. This cross-validation was repeated four times. For each combination of component-wise 

learning rates the optimal numbers of boosting iterations and their corresponding performance were 

averaged across cross-validation folds and replicates. Finally, the combination of component-wise 

learning rates with the best average performance and its corresponding average optimal number of 

7 
8 

boosting iterations were identified. The maximum number of boosting iterations allowed was 20,000, so 

models with an optimal number of boosting iterations near 20,000 should be interpreted with caution as 

their performance may have improved with additional boosting iterations.

2.5.7 Model selection and performance

Four different models were fitted for each species-season combination. The four models represented all 

combinations of the two likelihoods (ZIP and ZINB) and the two methods of calculating the boosting 

offsets (intercept-only model and generalized additive model). 

The performance of the four models was evaluated from a suite of nine performance metrics. The first 

performance metric was percent deviance explained (PDE), which is essentially the percentage of 

variation in the data explained by the fitted model beyond the amount of variation explained by a simpler 

model without predictor variables. PDE indicates overall model fit. PDE is somewhat analogous to the 

more familiar R2 metric for a linear regression. The second performance metric was the area under the 

receiver operating characteristic curve (AUC), which indicates how well a model predicts binary data. We 

calculated AUC by converting the count data to presence/absence data. The AUC metric indicates how 

well the models predicted the observed presence of a species, but not necessarily how well the models 

predicted the relative density of a species. The third and fourth performance metrics were the Spearman 

and Gaussian (Boudt et al. 2012; Bodenhofer et al. 2013) rank correlation coefficients between the 

observed and predicted data. These metrics indicate how well the model predictions of relative density for 

each transect segment correlate with the number of individuals counted. The fifth and sixth performance 

metrics were the median and mean absolute residual error as a percentage of the mean count. The seventh

and eighth performance metrics were the median and mean residual error (predicted minus observed) as a 

percentage of the mean count, referred to as ‘bias’. Negative bias values indicate that predicted density 

was lower than observed, while positive bias values indicate that predicted density was higher than 

observed. The final performance metric was root mean square error (RMSE), the square root of the 

average squared residual error.

It is important to recognize that the model performance metrics mainly reflect the statistical fit of the 

models to the data. They reflect only the data that were analyzed, and they do not reflect the quality of 

model predictions away from the data. For example, the survey data did not cover the full extent of the 

study area, so some model predictions are essentially interpolations/extrapolations from data in other parts 

of the study area. The accuracy of those predictions is not necessarily reflected by the model performance 

metrics. Nevertheless, the performance metrics give an indication of how accurately a model was able to 

predict the observed data, and good performance provides a measure of confidence in the modeled 

distributions, especially within the temporal and spatial coverage of the observed survey data.

The ‘best’ of the four models for each species-season combination was chosen on the basis of a subset of 

the performance metrics: PDE, AUC, Gaussian rank correlation, median and mean absolute residual error, 

and RMSE. The values of each of these performance metrics were ranked across the four models, and the 

model with the highest sum of ranks across these performance metrics was chosen as the best model. In 

the case of a tie the model with the highest PDE was chosen.

A small percentage of the ‘best’ models (4%; 5 models) had highly problematic predictions: Common 

Eider winter, Dovekie winter, Red Phalarope summer, Roseate Tern fall, and Sooty Shearwater fall. 

Problems included highly skewed distributions of predicted density across bootstrap replicates and space, 

flat distributions of predicted density across space, spatial patterns that did not make sense ecologically, 

and erroneous predicted density values in some bootstrap replicates (e.g., unrealistically high values). The 
9 

problematic models for Common Eider, Dovekie, Red Phalarope, and Sooty Shearwater were ZIP

models. These species often form large aggregations, which could have contributed to the observed 

problems, especially with Poisson models that do not allow for overdispersed count data. For three of the 

five problematic species-season combinations we chose the second best model instead and for two we 

chose the fourth best model. 

2.5.8 Spatial prediction

The final fitted model for each species and season was used to predict relative density throughout the 

study area. It is important to recognize that the model predictions do not represent actual density. During 

visual surveys, individual birds may be missed either because they are below the surface of the water 

(availability bias) or simply because observers failed to notice them (perception bias) (Barlow 2015). The 

failure to count some individuals biases estimates of density downward relative to actual density. Animal 

movement can also bias estimates of density. Birds may be attracted or repelled by ships, small boats, and 

planes biasing estimates upward or downward, respectively. Flying birds can also bias estimates, with the 

direction of the bias depending on the speed and direction of the animals’ movement relative to those of 

the survey platform (Spear et al. 1999). Furthermore, non-randomized survey data can result in biased 

estimates of abundance and density in species distribution models (Conn et al. 2017). Our model 

predictions should only be interpreted as indices of density. The model predictions can be interpreted as 

being proportional to the expected number of individuals per km2

, but they do not represent the actual 

expected number of individuals per km2

. 

The predicted relative density in a given grid cell corresponds to predictions for a transect segment whose 

mid-point falls within that grid cell. Spatially explicit predicted values were calculated for each cell of the 

study grid from the values of the spatially explicit predictor variables for that cell. All other predictor 

variables were set to their mean values. The predictions integrated the zero-inflated and Poisson/negative 

binomial components of the likelihood. 

2.5.9 Variable importance

While determining the ecological drivers and mechanisms behind the spatial distributions of marine bird 

species in the study area was not our primary objective, our model results do provide some indication of 

which variables were most useful for predicting those distributions. Those variables may provide useful 

starting points for future studies aimed more at ecological inference.

We calculated the relative importance of a given predictor variable in the final fitted models by summing 

the decrease in the negative log-likelihood in each boosting iteration attributable to that predictor variable. 

Thus, variable importance reflects the frequency with which a given predictor variable occurred in the 

selected base-learners across boosting iterations and that variable’s ability to explain variation in the data 

when it was selected. When multiple predictor variables occurred in the selected base-learner for a given 

boosting iteration, the decrease in the negative log-likelihood was divided evenly among those predictor 

variables. Relative variable importance was re-scaled so that it summed to 1 across predictor variables.

2.5.10 Uncertainty

Uncertainty in model predictions was estimated using a non-parametric bootstrapping framework. For 

each bootstrap iteration, the set of unique transect IDs was resampled with replacement, and the data for

each transect ID were assigned weights proportional to the frequency of that ID in the sample. These data 

weights were then applied when fitting the model during that bootstrap iteration. Predictor variables that 

were not included in the final model were excluded from the bootstrap analysis. Two hundred bootstrap 
10 

iterations were conducted producing a sample of predictions from which we calculated means, standard 

errors (SE), coefficients of variation (CV), 5%, 25%, 50%, 75%, and 95% quantiles, and 50% and 90% 

confidence interval widths to characterize uncertainty in the predictions. The 50% and 90% confidence 

intervals were defined as the intervals between the 25% and 75% quantiles and between the 5% and 95% 

quantiles, respectively.

As with the model performance metrics, the estimated uncertainty in the model predictions is conditional 

on the model and the data. It does not capture all of the uncertainty associated with our model predictions. 

Nevertheless, the estimated uncertainty is an important indication of the precision of the model 

predictions, and it should be an integral consideration when using the model predictions. 

The bootstrapped median predictions (50% quantile) were generally better (e.g., smoother density 

gradients and fewer artifacts) than the bootstrapped mean predictions and the non-bootstrapped 

predictions, so we chose the bootstrapped median predictions as the best representation of the predicted 

spatial distributions of relative density.

A small number of bootstrap iterations in seven different models produced unreasonably large predictions 

in a few grid cells that caused numerical issues when calculating the bootstrapped mean, SE, and CV. 

These extreme outlying values were removed when calculating these bootstrapped quantities for these 

models. The models that exhibited this issue were: Brown Pelican winter, Common Eider fall, Great 

Black-backed Gull winter, Herring Gull winter, Leach’s Storm-Petrel spring, Long-tailed Duck spring, 

and Red-throated Loon winter. All but one of these models were ZIP models suggesting that the issue 

may have been related to bootstrapping models with that particular likelihood. 

2.5.11 Implementation

The analysis was conducted with R 64-bit (R Core Team 2017) version 3.3.3 (or 3.4.4 for model re-runs; 

see Section 2.2) running under CentOS Linux 7 using custom scripts that relied on multiple existing 

contributed packages, including ‘boot’ version 1.3‐17 (Canty and Ripley 2015), ‘DBI’ version 0.6-1 (RSIG-DB et al. 2017), ‘fields’ version 8.10 (Nychka et al. 2015), ‘gsubfn’ version 0.6-6 (Grothendieck 

2014a), ‘lattice’ version 0.20-31 (Sarkar 2008), ‘maps’ version 3.1.1 (Becker et al. 2016), ‘maptools’ 

versions 0.8‐36 (Bivand and Lewin‐Koh 2015), ‘MASS’ versions 7.3‐42 (Venables and Ripley 2002), 

‘Matrix’ version 1.2‐8 (Bates and Maechler 2017), ‘mboost’ version 2.4‐2 (Hothorn et al. 2015), ‘mgcv’ 

version 1.8-17 (Wood 2006), ‘modeltools’ version 0.2‐21 (Hothorn et al. 2013), ‘mvtnorm’ version 1.0-2 

(Genz et al. 2014), ‘nlme’ version 3.1-131 (Pinheiro et al. 2017), ‘party’ version 1.0‐21 (Hothorn et al. 

2006), ‘pROC’ version 1.8 (Robin et al. 2011), ‘proto’ version 0.3-10 (Grothendieck et al. 2012), ‘pscl’ 

version 1.4.9 (Jackman 2015), ‘raster’ version 2.4‐15 (Hijmans 2015), ‘reshape’ version 0.8.5 (Wickham 

2007), ‘rgdal’ version 1.0‐4 (Bivand et al. 2015), ‘rgeos’ version 0.3‐11 (Bivand and Rundel 2015), 

‘rococo’ version 1.1.2 (Bodenhofer et al. 2013), ‘RSQLite’ version 1.0.0 (Müller et al. 2014), ‘sandwich’ 

version 2.3-3 (Zeileis 2006), ‘sp’ version 1.1‐1 (Pebesma and Bivand 2005), ‘spam’ version 1.3-0 (Furrer 

and Sain 2010), ‘sqldf’ version 0.4‐10 (Grothendieck 2014b), ‘stabs’ version 0.5-1 (Hofner and Hothorn 

2015), ‘strucchange’ version 1.5-1 (Zeileis et al. 2002), ‘VGAM’ version 0.9‐8 (Yee 2015), and ‘zoo’ 

version 1.8-0 (Zeileis and Grothendieck 2005).

2.6 Map display 

Spatial predictions within the study area are displayed as a pair of maps for each species and season (Fig. 

5). Seasons were defined as spring (March-May), summer (June-August), fall (September-November), 

and winter (December-February). The first map for each species and season displays the predicted 
11 

relative density (bootstrapped median predictions), and the second map displays the CV of predicted 

relative density, both at a 2-km spatial resolution. The first map also contains an inset map of observed 

relative density: the total count divided by the total area surveyed at a 10-km spatial resolution. 

Color spectrums are employed to visualize spatial variation in values across the study area. For relative 

density the colors range from blue (lower density) to red (higher density), and for the CV the colors range 

from light brown (lower CV, more precise) to dark brown (higher CV, less precise). The number ranges 

corresponding to each color are indicated in the map legends. The break points between the number 

ranges were evenly distributed on the natural log scale. To emphasize the relative nature of the model 

predictions, the number ranges displayed in the legend of each predicted relative density map were rescaled by dividing them by the maximum predicted value for that map. Thus, predicted relative density 

values on each map range from 0-1, representing the lowest to the highest predicted relative density for 

that map. For the inset observed density maps, values greater than the maximum predicted value were 

included in the highest legend category. We encourage others making maps from the model predictions 

provided by this project to similarly re-scale the legends to emphasize the relative nature of the 

predictions and to prevent misinterpretation of the predictions as actual density. 

Both maps also display a hatched overlay indicating areas without survey effort at a 10-km spatial 

resolution. Model predictions in areas without survey effort are essentially extrapolations and should be 

interpreted with extreme caution. We recommend additional field surveys in these areas to validate model 

predictions. 

2.7 Expert review of maps 

We solicited feedback on the accuracy of the maps from 44 experts with experience and knowledge of 

marine birds in the study area. These experts were from a range of organizations including federal and 

state government agencies, academic institutions, non-profits, and consultants. We received comments 

from nine of these experts, and their feedback was incorporated into this report.

2.8 Predicted proportional relative abundance by BOEM wind energy 

planning/lease area 

To facilitate comparisons of predicted relative density among BOEM wind energy planning/lease areas

the predicted relative abundance (density × area) in each BOEM area as a proportion of the total relative 

abundance in the entire study area was tabulated by BOEM area and season for each species. These 

calculations accounted for any partial coverage of study grid cells. Computer files defining the BOEM 

areas for this analysis were downloaded from https://www.boem.gov/Renewable-Energy-GIS-Data/, 

https://www.boem.gov/OCS-A-502/, and https://www.boem.gov/OCS-A-503/ in September 2017. It is 

important to recognize that these proportional relative abundance results are relative to the study area.

They do not account for birds outside of the study area. If the entire range of a species were to be 

considered then the proportional relative abundance in BOEM areas could decrease.

3. Results

3.1 Model selection and performance 

The final selected models varied across species and seasons (Table 6). ZINB models were selected more 

frequently than ZIP models (85 versus 55 models), and models with boosting offsets derived from a 
12 

generalized additive model were selected more frequently than models with boosting offsets derived from 

a simple intercepts-only model (104 versus 36 models). A higher proportion of selected ZINB models had 

boosting offsets derived from a generalized additive model than did selected ZIP models. Models with 

different likelihoods and different boosting offset methods were often selected for different seasons for 

the same species.

The statistical performance of the final selected models was highly variable among species and seasons

(Table 6). Over half of the models had an optimal number of boosting iterations >19,000, close or equal 

to the maximum of 20,000. PDE ranged from 14-86%. AUC ranged from 0.81 to 1.00. Spearman and 

Gaussian rank correlation coefficients ranged from 0.05-0.59 and 0.06-0.60, respectively. The median and 

mean absolute residual error ranged from 0-56% and 39-262% of the mean count, respectively. The 

median and mean residual error (bias) ranged from 0-56% and -100 to 122% of the mean count, 

respectively. RMSE ranged from 0.05-426. 

Performance metrics were often correlated with each other, usually in the expected direction, although the 

correlations were often not very strong. The Spearman and Gaussian rank correlation performance metrics 

were very strongly correlated with each other (Spearman rank correlation coefficient rs = 0.99) as were 

the median residual error and median absolute residual error (rs = 0.98). PDE was strongly correlated (|rs| 

> 0.5) with AUC (positive), median/mean absolute residual error (negative), and median residual error 

(negative). AUC was strongly negatively correlated with median absolute residual error and median 

residual error. Spearman/Gaussian rank correlation coefficient were strongly negatively correlated with 

mean absolute residual error. Mean absolute residual error was strongly positively correlated with median 

residual error and median absolute residual error. Unexpectedly, Spearman/Gaussian rank correlation 

were positively correlated with RMSE and to a lesser extent negatively correlated with AUC. 

The best performing models overall were Black-capped Petrel fall and winter, which had 7 of 9 

performance metrics in the top 10 across all species and seasons. Audubon Shearwater winter had 5 

performance metrics in the top 10. These three models were all ZIP models with boosting offsets derived 

from a generalized additive model. The models with the worst performance overall were Parasitic Jaeger

spring and Pomarine Jaeger summer which had 7 of 9 performance metrics in the bottom 10 across all 

species and seasons. Parasitic Jaeger summer and fall had 6 performance metrics in the bottom 10, and 

Red-breasted Merganser winter and Roseate Tern fall had 5 metrics in the bottom 10. Some models had 

performance metrics in the bottom 10 and top 10. For example, Parasitic Jaeger spring and summer and 

Pomarine Jaeger summer had RMSE in the top 10. It is possible that low sample size contributed to the 

poor performance of some models. For example, Parasitic Jaeger spring had the lowest number of 

segments with sightings, the lowest number of individuals counted, and the lowest count per segment 

when sighted (Table 2). Parasitic Jaeger summer had the second lowest total count. 

3.2 Predicted spatial distributions 

3.2.1 Spatial and seasonal patterns

Given the large number of species modeled the predicted spatial distributions of relative density varied 

widely (Appendix C). Broad patterns in predicted distributions generally aligned with observed sighting 

locations and densities and usually matched what is known about the distributions of these marine bird 

species. Predicted distributions of Brown Pelican, Double-crested Cormorant, Horned Grebe, loons, Redbreasted Merganser, sea ducks, and several gull and tern species were relatively coastal with highest 

relative density near shore. Black-capped Petrel, Bridled Tern, Northern Fulmar, shearwaters, Sooty Tern, 
13 

and storm-petrels had highest densities further offshore. Alcids tended to have offshore distributions 

except during the summer (June-August) for species that breed in the study area (Atlantic Puffin, Black 

Guillemot, and Razorbill). Jaegers and skuas also had offshore distributions except Parasitic Jaeger during 

summer when the highest density was nearer to shore overlapping with the area of highest density for 

some breeding species (e.g., Arctic Tern, alcids). Areas of highest density for Black-legged Kittiwake, 

Great Black-backed Gull, Herring Gull, and Northern Gannet generally spanned nearshore to offshore. 

Distributions of phalaropes tended to be closer to shore during summer than during spring (March-May)

and fall (September-November). 

The predicted distributions of some species within the study area were relatively restricted in area, at least 

seasonally. For example, most of the highest density areas for Arctic Tern, breeding alcids, and 

phalaropes during summer were relatively nearshore off Maine and Nova Scotia around the mouth of the 

Bay of Fundy. Distributions of sea ducks often exhibited localized areas of highest density especially 

during winter (December-February), for example in and around Nantucket Sound, and near the mouths of 

Delaware Bay and Chesapeake Bay. The distribution of the endangered Roseate Tern exhibited relatively 

small areas of highest density in and around Cape Cod and Nantucket Sound during spring-fall but also 

further south along the east coast of Florida during summer. Most other species’ distributions were less 

restricted in size but their highest densities were usually still limited to parts of the study area. The bulk of 

many species’ distributions was in the northern part of the study area within the Gulf of Maine and over 

surrounding shelf waters out to the shelf break or beyond and south as far as Cape Hatteras (e.g., Atlantic 

Puffin, Black-legged Kittiwake, Dovekie, Great Black-backed Gull, Great Shearwater, Herring Gull, 

Northern Fulmar, and Sooty Shearwater). Other species had more southerly distributions, either coastal 

(e.g., Brown Pelican, Royal Tern) or offshore (e.g., Black-capped Petrel, Sooty Tern).

The predicted spatial distributions of many species changed seasonally reflecting seasonal movements 

and migrations. As mentioned above, the seasonal distributions of alcids that breed in the study area were

more concentrated nearshore around the colonies during summer transitioning to more widely dispersed 

offshore distributions during winter, although the area of highest density for Razorbill during winter and 

spring was around and to the south of Cape Cod and Nantucket Sound. Areas of relatively high density 

extended further north during the summer and/or further south during the winter for many species (e.g., 

Common Loon, Laughing Gull, Northern Fulmar, and Northern Gannet). The distributions of Great 

Black-backed Gull and Herring Gull were closer to shore, more concentrated, and more northerly during 

summer and were further from shore, more dispersed, and more southerly during winter. The distributions 

of both phalarope species exhibited similar seasonal changes with highest densities along the shelf break 

during spring, closer to shore off Maine and Nova Scotia during the summer, and from the mouth of the 

Bay of Fundy through the Gulf of Maine and the Northeast Channel between Georges Bank and Browns 

Bank extending along the shelf break during winter. Migrations of species that mostly leave the study 

area altogether were reflected in the absence of models for certain seasons (e.g., sea ducks during summer 

and jaegers, phalaropes, skuas, storm-petrels, and terns during winter).

Predicted relative abundance in individual BOEM wind energy planning/lease areas expressed as a 

proportion of total relative abundance in the study area reflected the broader predicted spatial 

distributions, but was additionally influenced by the size of individual BOEM areas with larger areas 

having higher relative abundance for the same density (Appendix D). Relative abundance in individual 

BOEM areas was usually a small proportion of the total relative abundance within the study area, 

although for 9 species this proportion was ≥1% for one or more BOEM areas in one or more seasons. 

Often these BOEM areas were New Jersey areas (Common Loon spring, Common Tern summer, Surf 
14 

Scoter spring) or the South Carolina Grand Strand area (Bonaparte’s Gull spring, Common Loon winter, 

Red-throated Loon spring, Royal Tern spring-summer). Twenty-two species had total relative abundance 

across all BOEM areas that was ≥1% of total relative abundance in the study area for one or more 

seasons. The highest of these proportions was for Razorbill and Common Loon during spring where total 

relative abundance across all BOEM areas was >5% of the total relative abundance in the study area, 

mostly in Rhode Island and Massachusetts areas for Razorbill and in New Jersey areas and the South 

Carolina Grand Strand area for Common Loon. 

3.2.2 Uncertainty

Generally speaking, the estimated uncertainty in the model predictions was high. The minimum CV for 

predictions in individual spatial grid cells ranged from 0.02-1.34 across species and seasons, and the 

maximum CV ranged from 1.23-14.14. Thus, all models had at least some predictions with CVs 

exceeding 1 indicating that the SE was greater than the mean. These CVs reflect multiple aspects of the 

data and model including the amount of survey effort (i.e., sample size), sampling variability and the 

degree to which birds of a given species aggregate, extrapolation to predictor variable values outside the 

range of values associated with the observed survey effort, and un-modeled variation in the spatial 

distributions of birds (e.g., changes in spatial distributions over time). It is difficult to quantify the 

contribution of each of these sources of variability to the CV of any particular prediction. The spatial 

patterns in the estimated CV of predictions varied widely across species and seasons (Appendix C). In 

some cases the CV was relatively lower in areas with predicted higher density (e.g., Northern Fulmar 

summer), but in other cases the CV was higher in areas with predicted higher density (e.g., Black 

Guillemot summer). Higher CVs sometimes reflected an absence of survey effort in areas with predicted 

low (Ring-billed Gull winter, offshore) or high relative density (Atlantic Puffin summer, off the coast of 

Nova Scotia). Many of the problematic predictions (see Section 3.2.3) were associated with relatively 

high CVs (e.g., Brown Pelican offshore in winter, Herring Gull over northeast seamounts in winter). The 

spatial patterns in the estimated CV often reflected the influence of particular predictor variables, even 

when that influence was not as apparent in relative density itself (e.g., Red-throated Loon winter, 

probability of anticyclonic eddies). The estimated SE and confidence interval widths of predicted relative 

density were highly correlated with the relative density values, usually resulting in very similar spatial 

patterns. As indicated by the magnitude of the CVs, the values of the SE and confidence interval widths 

were generally large relative to the corresponding predictions of relative density.

3.2.3 Problematic predictions

Predictions of relative density for some species in some areas were questionable, especially in areas 

without survey effort. These predictions were often associated with high CVs. There was much less 

survey effort beyond the shelf break than coastally and over the shelf, especially during winter (Fig. 3). 

The predicted densities of many species were relatively high in this offshore area, some of which are 

plausible but others of which are doubtful. For example, several coastal species had moderate relative 

densities extending far offshore during winter that are unrealistic (e.g., Brown Pelican, Horned Grebe, 

Red-breasted Merganser, Red-throated Loon). Predictions of even moderate relative density of sea ducks 

far offshore (e.g., Long-tailed Duck, Surf Scoter, and White-winged Scoter in fall) are biologically 

implausible given their benthic foraging habitat (White et al. 2009; Veit et al. 2016; White and Veit

2018). Even for species that are generally found far from shore, areas of moderate to high relative density 

offshore may be larger than is realistic given the limited survey coverage there (e.g., Black-capped Petrel, 

Audubon Shearwater, Manx Shearwater). Predicted moderate to high relative densities of other species in 

some offshore areas are also unrealistic or questionable at best (e.g., Common Murre in winter, Dovekie 

along the edge of the Blake Plateau in spring and fall, Leach’s Storm-Petrel in fall, Least Tern in fall, 
15 

South Polar Skua in fall). The extreme southern part of the study area (Straits of Florida) had very little 

survey effort (Fig. 3). Several northern species had moderate to high relative densities in this area that are

unrealistic (e.g., Black Scoter, Dovekie, Northern Fulmar, Surf Scoter, Thick-billed Murre). Several 

species had high relative densities close to the southern shore of Nova Scotia in an area without survey 

effort (Arctic Tern, Atlantic Puffin, Black Scoter, Razorbill, Surf Scoter). There was predicted moderate 

density of Brown Pelican in Canadian Bay of Fundy waters during summer that is further north than is 

realistic for this species’ current distribution.

Predicted spatial distributions sometimes exhibited patterning or extremely localized areas of high relative 

density that are likely unrealistic. Unrealistic patterning often reflected the undue influence of individual 

predictor variables on model predictions, especially in areas with little or no survey effort. These 

predicted features were often associated with higher CVs. For example, it is apparent that the predicted 

distribution of Brown Pelican in winter was heavily influenced by the probability of anticyclonic eddies 

with large discrete patches of higher density offshore corresponding to areas with high eddy probabilities. 

The distribution of Black-legged Kittiwake in spring exhibits an obvious effect of distance to land 

resulting in a line of discrete change in relative density paralleling the coast. Neither of these patterns are 

realistic. The distribution of Bridled Tern in summer exhibits a circular area of high relative density in the 

South Atlantic Bight corresponding to an area of high probability of cyclonic eddies. Some predictions 

were clearly extrapolations of the relationship between bathymetric predictor variables and density to 

areas without survey effort. For example, the predicted distributions for Audubon Shearwater in fall,

Cory’s Shearwater in fall, Great Shearwater in spring, and Herring Gull in winter exhibit well-defined 

areas of high relative densities over seamounts in the northeast part of the study area. It is probable that at 

least the discreteness of these distribution features are unrealistic. Predictions for Parasitic Jaeger, 

Wilson’s Storm-Petrel, and other species exhibit ribbons of moderate to high density that match high 

slope transitions extending from north to south throughout the entire study area. These predictions are 

more likely realistic along the shelf break in the north where there was more survey effort than they are 

further south along the edge of the Blake Plateau further offshore in an area with very little survey effort.

The predicted distribution for Long-tailed Duck in the spring exhibited fairly high density over Georges 

Bank in an area with some survey effort but no sightings. That prediction was primarily a result of

similarities in spring chlorophyll-a and turbidity between Georges Bank and the area of observed high 

Long-tailed Duck density in Nantucket Sound and over Nantucket Shoals. 

The very highest densities were sometimes limited to one or a few very small areas suggesting localized 

long-term concentrations of abundance that are almost certainly unrealistic. These predictions were 

sometimes associated with higher CVs but not always. Some example species and seasons exhibiting this 

behavior were Arctic Tern (summer, spot over Georges Bank), Black-Legged Kittiwake (fall), Common 

Murre (spring), Herring Gull (winter), Manx Shearwater (summer), Red-necked Phalarope (spring), and 

Thick-billed Murre (winter). It is difficult to ascribe these highly localized distribution features to a 

specific cause, but in at least some cases they may partially reflect large temporal and spatial aggregations 

of birds that coincided with survey effort rather than average long-term spatial patterns per se. It is 

unlikely that these are persistent areas of higher density than adjacent areas. While our spatial predictive 

modeling framework theoretically accounts for effort and attempts to account for the aggregated nature of 

bird distributions and sightings, limited sample size combined with extreme aggregations can unduly 

influence model predictions.

In addition to over-predictions of relative density, it is possible that the model predictions 

underrepresented the relative density of some species in some areas. For several species that frequent 
16 

Stellwagen Bank and Cape Cod in the southwestern Gulf of Maine the models did not always indicate 

substantial relative density there (e.g., Common Murre in winter and spring, Cory’s Shearwater in 

summer and fall, Common Tern in fall, Great Shearwater in summer, Pomarine Jaeger in fall, and Roseate 

Tern fall). If the models indeed underrepresented relative density in this area, it could have been a result 

of survey timing and limited sightings in the data analyzed here or an issue with model fit. It is also 

possible that in some cases the use of 3-month seasons diluted the signal from more short-term 

concentrations of abundance (e.g., migrating birds that are only present for the first month or so of fall). 

As another example, Razorbill are known to occur in high densities in the Bay of Fundy during winter 

(Huettmann et al. 2005), yet predicted relative density in that area during winter was low potentially 

suggesting an under-prediction. There was no survey effort in that area during winter, which possibly 

explains the predicted low relative density.

Some of the predicted distributions for sea duck species were flatter than expected with moderate relative 

density extending further offshore and further north or south than the sighting data would suggest. The 

distributions for Common Eider in spring and winter most clearly manifested this issue. The areas of 

highest relative density corresponded with the highest observed densities in and around Nantucket Sound, 

but moderate relative densities extended offshore throughout the Gulf of Maine (winter) or the entire 

study area (spring) where there were few if any sightings. The predicted relative density values in these 

problematic areas were quite high (e.g., as high as 10 birds per km2

 or more in winter), emphasizing that 

the model predictions should not be interpreted in absolute terms. Other sea duck distributions exhibited 

the same issue to a lesser extent (e.g., Long-tailed Duck winter). This issue was likely a result of the 

highly skewed distributions of observed counts for these species. For example, the maximum observed 

density on a single transect segment was >40,000 birds/km2

 for Common Eider in winter and spring and 

Long-tailed Duck in spring while the mean observed densities when sighted ranged from only 73-313 

birds/km2

. Although the model likelihoods allowed for overdispersion in the count data they might not 

have been sufficient for species and seasons with extremely large aggregations.

These problematic predictions sometimes affected the predicted proportional relative abundance by 

BOEM wind energy planning/lease area. For example, the relative abundance of Common Eider in 

individual BOEM areas was a small proportion of the total relative abundance within the study area 

(Appendix D), but the proportional relative abundance values for Common Eider during spring were more 

evenly distributed across BOEM areas from north to south than is realistic. There were no sightings south 

of New Jersey.

3.3 Predictor variable relative importance 

The modeling framework used in this study was designed to provide accurate predictions. It was not 

designed to identify which environmental predictors were most ecologically relevant in determining the 

distribution of birds. Ecological inferences from the variable importance results should be cautious. 

Nevertheless, these results may suggest interesting hypotheses for future research. Also, as discussed in 

Section 3.2.3 extrapolations of the predicted relationships between density and predictor variables to areas 

with little or no survey effort were often questionable or unrealistic.

The most important predictor variables varied across species and seasons (Fig. 6a-d). Day of year, mean 

sea surface temperature, mean chlorophyll-a, turbidity, and distance to land were in the top eight most 

important predictor variables averaged across all species in all four seasons for both the p and μ model 

components. Year was almost always in the top eight except for the p model component of fall models.

Day-of-year was the most important predictor on average in spring and fall models likely reflecting the 
17 

movement of migratory species in and out of the study area during the transitions between summer and 

winter. The importance of the year effect reflects inter-annual variability in species’ seasonal abundance 

within the study area. Mean sea surface temperature was a relatively important predictor for most models, 

reflecting its ability to capture broad latitudinal and Gulf Stream-related patterns in distribution. For 

example, all seasonal Sooty Shearwater models predicted an increase in density as sea surface 

temperature decreased corresponding to the occurrence of this species mainly in the northern part of the 

study area. The effects of mean chlorophyll-a and turbidity tended to be most important in models for 

coastal species (e.g., Brown Pelican, Double-crested Cormorant, Horned Grebe, some gull species, loons, 

Red-breasted Merganser, sea ducks, and some tern species), reflecting these predictors’ ability to capture 

coastal distributions, especially concentrations of abundance in certain areas (e.g., in and around 

Nantucket Sound). For example, Black Scoter models usually predicted a strong increase in density as 

chlorophyll-a and turbidity increased. Distance to land was often an important predictor for coastal 

species but was also often important for species with more offshore distributions (e.g., shearwaters).

Other predictor variables that were relatively important for multiple species in multiple seasons were 

transect ID, sea surface temperature standard deviation (SD), and wind stress (particularly in the east-west 

direction). Transect ID was an especially important predictor in the μ model component, and its 

importance indicates substantial variation in the count data that was not explained by the temporal and 

spatial predictor variables. Sea surface temperature SD was highly correlated with mean sea surface 

temperature (Table 4) so it is difficult to separate the effects of these two predictors. Wind stress was 

more eastward and/or stronger further north and further offshore (Appendix B), so some seasonal models 

with more southern spatial distributions predicted a negative relationship between density and wind stress 

in the east-west direction (e.g., Brown Pelican winter, Least Tern fall, Royal Tern fall) while some 

seasonal models with more northern spatial distributions predicted a positive relationship between density 

and wind stress in the east-west direction (e.g., Common Murre winter, Leach’s Storm-Petrel fall). Wind 

stress was more often important in fall and winter models when the westerlies are stronger and the northsouth gradient in east-west wind stress is greater.

Generally speaking the climate index predictor variables were relatively unimportant. They were more 

important in the spring and fall models than in the summer and winter models. The Atlantic Multidecadal 

Oscillation (AMO) index and the Trans-Niño Index (TNI) were the most important of the four indices on 

average including their 1-year-lagged versions, although the Multivariate El Niño-Southern Oscillation 

Index (MEI) and the North Atlantic Oscillation (NAO) index were important in some models.

4. Discussion

4.1 Interpretation of maps 

This report presents maps of the seasonal spatial distributions of 47 marine bird species in U.S. Atlantic 

OCS waters (Appendix C) that can be used to inform marine spatial planning in the region. The maps of 

predicted relative density are accompanied by corresponding maps of the estimated CV of model 

predictions. It is important that these uncertainty maps are considered alongside the relative density maps. 

In many cases the CVs are very large indicating substantial statistical uncertainty and variability 

associated with the corresponding predictions, and those predictions should be interpreted cautiously. 

Although model predictions are at a 2-km resolution, interpretation of the maps presented here to inform 

spatial planning is probably more reliable at scales of 10-100 km. The maps of predicted relative density 

are also accompanied by inset maps of observed relative density. These inset maps indicate where species 
18 

were seen during a given season by the surveys analyzed here, and this information should also be 

considered when interpreting model predictions.

The maps represent spatial predictions of long-term average relative density. They do not provide 

predictions of the actual number of individuals of a given species that would be expected in a given area; 

they only indicate where a given species may be more or less abundant. Also, the maps do not provide 

predictions of density at a specific time; they only indicate seasonal distributions averaged across the 

timeframe of the survey dataset (1978-2016). The spatial distributions of marine birds often change 

during the course of a season and from year to year in response to fluctuations in environmental 

conditions, prey distributions, and fisheries (Oedekoven et al. 2001; Ballance et al. 2006; Renner et al. 

2013; Nisbet and Veit 2015; Gjerdrum and Bolduc 2016; Goyert et al. 2016; Suryan et al. 2016; Veit et 

al. 2016; White and Veit 2018), so it is likely that the spatial distributions of many if not all of the 

modeled species have changed over time. Furthermore, long-term average spatial distributions do not 

necessarily highlight areas that are consistently used by a large number of birds but only for a short period 

during the year (e.g., hours or days); for example, movement corridors. By definition, areas of high longterm relative density should reflect relatively large numbers of birds using those areas for longer periods.

The maps represent predicted density relative to the study area only. The study area excluded very 

nearshore waters and coastal sounds and bays (e.g., Bay of Fundy, Long Island Sound, Delaware Bay, 

Chesapeake Bay) that are used by multiple marine bird species, sometimes with large numbers of birds 

occurring in these areas at certain times of year (e.g., wintering waterfowl). More generally, areas of high 

density occur adjacent to the study area for some species, and these areas are not represented on the maps.

Thus, the areas with the highest predicted relative densities are not necessarily the areas of highest density 

for a species were its entire range to be considered. Depending on data availability, habitat-based 

predictive modeling could theoretically be extended to very nearshore waters and coastal sounds and 

bays, but such modeling would entail additional considerations. For example, relationships between bird 

density and environmental predictors may be different in these areas than in more offshore areas. Another 

consideration would be an appropriate spatial modeling technique for complex coastlines (Scott-Hayward 

et al. 2014). 

4.2 Data limitations and information gaps 

The maps presented here represent model predictions that ultimately rely on the underlying survey data 

that the models were fit to. The distribution of survey effort was highly uneven across the study area and 

over time (Figs 2 and 3), so some regions and years were much better sampled than others. For example, 

there were many fewer data offshore than nearshore. Model predictions in areas with no survey effort are 

indicated by hatched overlays on the maps. Predictions in these areas should be interpreted with caution, 

and we recommend additional field surveys to validate them. More generally, additional survey effort in

under-sampled areas (e.g., the southern part of the study area; Fig. 3) would help to improve future model 

predictions. In order for future survey data to be amenable to an analysis such as that presented here, 

surveys should adhere to standardized survey protocols for marine birds at sea and document key 

information related to survey effort (e.g., complete spatial descriptions of transects, duration of 

observation, strip width, etc.). 

Future habitat-based predictive modeling could also benefit from data on additional environmental 

predictor variables. Most of the environmental predictors used in this study would be indirect drivers of 

marine bird distributions. Data on variables that birds are responding to more directly (e.g., prey) would 

likely improve model predictions, particularly in areas with little or no survey effort. Unfortunately, 
19 

comprehensive climatologies of prey densities across the entire study area were not available for this 

study. 

4.3 Combining seasonal maps 

The seasonal maps presented in Appendix C represent periods of three months of the year for individual 

species. For some applications it may be of interest to consider an annual average distribution for 

individual species or a seasonal/annual average distribution across species. Such combining of seasonal 

maps should be done cautiously while acknowledging important assumptions. 

To derive an annual map for a single species, we suggest averaging the seasonal predicted grids for that 

species assuming zero relative density for seasons that were not modeled. Ideally for bootstrapped 

predictions, a bootstrapped sample of annual predicted grids would be derived first and then desired 

quantities calculated from that sample (e.g., median, quantiles, etc.). The averaging of seasonal grids 

entails several important assumptions. First, species are often present in the study area at low density 

during seasons that were not modeled, so the assumption of zero density during those seasons is an 

approximation. Second, the predictions in each seasonal grid correspond to the average year of the survey 

data for that season which varied across seasons (spring – 2002, summer – 2003, fall – 2004, winter – 

2006). Thus, the seasonal maps represent predictions for different years, so any variation in relative 

density across years will influence predictions averaged across seasons. Third, as discussed previously the 

predictions are of relative density and do not explicitly account for several potential biases (e.g., 

detectability). Any variation in these biases across seasons will influence predictions averaged across 

seasons. Fig. 7 shows annual maps for an example species, Atlantic Puffin. By averaging the seasonal 

predicted grids the annual map most reflects the seasons with the highest predicted relative density values,

in this case summer and winter. 

To derive multi-species maps, we suggest first normalizing the predicted grids for individual species and 

then averaging grids across species. Predicted relative density values are not comparable between species,

so it would be inappropriate to directly average predicted grids across species. Normalization scales each 

grid the same, which makes averaging of grids across species more appropriate. Predicted grids can be 

normalized by dividing the predicted relative density value for each grid cell by the sum of predicted 

relative density values in all grid cells. If the grid dimensions are identical among species, as they are 

here, dividing by the mean predicted density value will achieve an equivalent normalization. Averaged 

normalized predicted values do not have a straightforward interpretation, so as with all predictions of 

relative density presented in this report it is the relative differences in predictions across space that are 

relevant, not the actual numbers themselves. It is important to note that because of normalization each 

species will essentially contribute equally to the average grid. Such equal weighting for all species may be 

desirable in some applications, but in others it may be desirable to weight species differently, for example 

by their population abundance. The results presented here could theoretically be used to derive weighted 

averages, but such an approach would require supplementary information to inform the weighting (e.g., 

independent estimates of population abundance). Fig. 8 shows an example multi-species winter map for a 

‘coastal waterfowl’ species group as defined by Curtice et al. (2016). The group includes Black Scoter, 

Common Eider, Common Loon, Long-tailed Duck, Red-throated Loon, Surf Scoter, and White-winged 

Scoter.
20 

4.4 Species identification 

A fundamental assumption of the analysis presented here is that all species, when present, were recorded.

This assumption might have been violated on one or more surveys if observers were focused on particular 

groups of species and failed to record occurrences of other species. A related issue is the identification of 

observed birds to the species level. From some survey platforms, or again if observers were more focused 

on some species than others, some birds might have been less likely to be identified to species. Birds that 

can be difficult to identify to species include but are not limited to alcids (e.g., Common Murre and 

Razorbill), gulls, jaegers, loons, phalaropes, scoters, storm-petrels, and terns. Birds that were not 

identified to species were not included in the main analysis presented here.

Our statistical modeling framework allowed for differences in the expected count of a given species 

among survey platforms and transects, so theoretically the models could account for some differences 

arising from failure to record or identify species. However, if a geographic area was covered by a limited 

number of surveys or platforms, then it would have been difficult or impossible for the model to 

determine whether differences in counts in that area were because of fewer birds in that area or because of 

differences in species recording and identification in that area. 

Phalaropes provide an illustrative example of the effect of species identification on predicted relative 

density. Phalaropes are relatively small birds that can be difficult to identify to species from typical aerial 

surveys. Unidentified phalaropes accounted for 18%, 33%, and 44% of all phalaropes counted in spring, 

summer, and fall, respectively. A ‘multi-species’ model fit to pooled sighting data for Red Phalarope, 

Red-necked Phalarope, and unidentified phalaropes predicted areas of high relative density similar to 

those predicted by the individual species’ models but also predicted areas of high relative density that 

were not evident in the individual species’ maps. The spring multi-species model predicted a band of 

moderate to high relative density off South Carolina and North Carolina and the fall multi-species model 

predicted a band of moderate relative density off the east coast of northern Florida and Georgia (Fig. 9), 

neither of which were apparent in the spring and fall maps for Red Phalarope and Red-necked Phalarope 

(Appendix C). These areas of predicted moderate to high relative density corresponded to multiple 

sightings of phalaropes that were not identified to species, almost entirely from three aerial surveys: 

AMAPPS_FWS_Aerial_Spring2012, AMAPPS_FWS_Aerial_Fall2012, and 

AMAPPS_FWS_Aerial_Fall2013 (Table 1). The areas of predicted moderate to high relative density in 

spring ran along and between BOEM areas for South Carolina and North Carolina, with particularly high 

predicted relative density along the eastern edge of the South Carolina Grand Strand area. The results of 

this example multi-species model for phalaropes highlights that models fit to individual species data may 

fail to identify an area of high relative density when the number of surveys in that area is limited and 

identification to species is difficult. Pooling data across species, as was done in this example, is one 

option for identifying these areas, but inference is then restricted to the group of species as a whole.

Alternative modeling techniques in combination with appropriate survey data may provide some ability to 

make individual species inferences from sightings of unidentified species (e.g., Johnston et al. 2015). 

4.5 Comparison to other marine bird distribution models 

A wealth of information is available about the distribution of marine birds in US Atlantic waters (e.g., see 

Nisbet et al. 2013 for a recent review). Here we compare our model predictions to results of other recent 

habitat-based modeling studies in this region. Habitat-based modeling of the spatial distributions of 

marine birds is an active field of research (Tremblay et al. 2009). Multiple studies have developed habitatbased distribution models for some of the species modeled here in parts of our study area often using 
21 

some of the same survey data (Zipkin et al. 2010; Winiarski et al. 2013, 2014a, b; Goyert 2014; Lieske et 

al. 2014; Flanders et al. 2015; Balderama et al. 2016; Goyert et al. 2016; Sollmann et al. 2016). By 

focusing on a smaller area, fewer predictor variables, fewer species, or a subset of the survey data with 

more homogenous sampling protocols these studies were able to employ alternative modeling techniques 

like distance sampling (e.g., Winiarski et al. 2013), occupancy modeling (Flanders et al. 2015), 

hierarchical multi-species community models (e.g., Goyert et al. 2016), and temporally varying spatial 

distributions (e.g., Balderama et al. 2016). Many of these studies were more focused on inference about 

ecological drivers of marine bird distributions so they employed modeling techniques more appropriate 

for that objective, although one of the three techniques considered by Lieske et al. (2014) was a machine 

learning technique with some similarity to the framework used here. Goyert (2014) was able to use data 

on prey abundance that were not available for our entire study area and timeframe improving their ability 

to make ecological inferences. While the combination of geographical scope, spatial resolution, number 

and heterogeneity of survey data, and number of species modeled here precluded many of these 

alternative modeling techniques, there are some elements that could theoretically be incorporated into the 

modeling framework used here. For example, there are statistical likelihood distributions that may better 

capture the extremely skewed count distributions for some species than the negative binomial likelihood 

used here (Zipkin et al. 2014, 2015; Balderama et al. 2016). More fundamental changes to the modeling 

framework would require sacrificing other aspects. For example, allowing spatial distributions to vary 

over time would limit either the number of predictor variables that could be used or the timespan of 

survey data that could be analyzed because most predictors do not have complete data from the entire 

timespan of the survey data (Mannocci et al. 2017). Fewer predictors or fewer survey data could reduce 

the predictive ability of the models.

Despite methodological differences, there were some similarities between the predicted spatial 

distributions of marine birds presented here and those found by other studies, but there were also 

differences especially at a fine spatial scale. Winiarski et al. (2013, 2014a) predicted higher Common 

Loon wintering densities in the western part of the Rhode Island Ocean Special Area Management Plan

(RISAMP) area than in the eastern part, similar to the predictions presented here, although there are fine 

scale differences. The predicted winter distribution of Northern Gannet density in this area presented by 

Winiarski et al. (2014b) was less similar to our predicted distribution. Predicted spring distributions of 

Common Tern and Roseate Tern in Massachusetts waters presented by Goyert (2014) were quite different 

from those presented here, although there was an area of predicted high density of Roseate Tern east of 

Martha’s Vineyard in both studies. Goyert (2014) predicted areas of high density west of Martha’s 

Vineyard for both species, at the north end of Cape Cod for Common Tern, and inside Cape Cod Bay for 

Roseate Tern that were not apparent in our model predictions. Predicted seasonal density distributions for 

several species in the lower Mid-Atlantic Bight presented by Goyert et al. (2016) were qualitatively 

similar to the predicted distributions presented here. For example, predicted distributions of Common 

Tern (summer), Royal Tern (summer), Laughing Gull (fall), Dovekie (winter), Red-throated Loon (winter 

and spring), and Wilson’s Storm-Petrel (summer) all exhibited similar inshore-offshore gradients in 

density. However some of the predicted distributions presented by Goyert et al. (2016) exhibited finer 

scale patterning that was not evident in the predictions presented here (e.g., Common Tern summer 2013, 

Razorbill winter 2013-2014). Broad similarities between the predicted distributions presented here and 

predicted distributions from other modeling studies in the same area are perhaps not surprising especially 

when there is overlap in the data analyzed. It is more difficult to ascribe the observed differences to 

specific causes. One of the likely common causes of differences is the shorter timespan of data analyzed 

in other studies combined with intra-seasonal and inter-annual variation in marine bird distributions. For 
22 

example, the spring survey data analyzed by Goyert (2014) came from a 10-day period in May over two 

years, whereas our predicted spring distributions reflect the average distribution from March, April, and 

May over many years.

4.6 Comparison to tracking studies

Electronic tracking of marine birds and habitat-based modeling of their distributions from tracking data is 

a burgeoning field of research (Tremblay et al. 2009; Wakefield et al. 2009, 2017). By following 

individual animals for extended periods of time tracking studies produce dense longitudinal data on 

animal locations that provide information generally much more difficult if not impossible to obtain from 

at-sea surveys. For example, tracking data more easily allow for the identification of areas that are used 

by animals consistently but for only short periods of time (e.g., movement corridors). One can infer 

behavioral states (e.g., foraging or migration) from tracking data (Jonsen et al. 2005) and thereby partition 

space usage by behavior, although direct behavioral data are often collected during at-sea surveys. Many 

electronic tracking technologies monitor animals wherever they go and so are not limited in geographical 

coverage like many survey platforms. Common limitations of tracking data are that results can be

unrepresentative at the population level due to small sample size and that some species can be difficult to 

track (e.g., too small to carry electronic devices), although as sample sizes continue to increase and 

technology continues to improve these limitations are becoming less of an issue. Electronic tracking 

studies provide valuable information on the spatial distribution of marine birds that is complementary to 

the information provided by the analysis of at-sea survey data. 

Multiple species modeled in this study have been tracked within the study area. Spiegel et al. (2017) 

reported on tracking of Red-throated Loon, Surf Scoter, and Northern Gannet in mid-Atlantic U.S. waters 

and compared their estimated distributions to seasonal predicted distributions from an earlier phase of this 

project (Kinlan et al. 2016). There was general agreement between the distributions of all three species in 

winter, but the distributions were less similar in spring and fall. In at least one case, Red-throated Loon 

during fall, the predicted distribution presented here is a better match with the distribution presented by 

Spiegel et al. (2017), but differences remain for other species and seasons. Spiegel et al. (2017) outline 

several possible reasons for these differences, but perhaps the most obvious reason is our exclusion of 

coastal sounds and bays and the heavy usage of those areas by tracked individuals (e.g., Surf Scoter in 

Chesapeake Bay).

Loring et al. (2014) and Beuth et al. (2017) reported on tracking of sea ducks, Black Scoter and Common 

Eider, respectively, in waters off southern New England. The winter distribution of tracked Black Scoter 

individuals (Fig. 4 of Loring et al. 2014) bore some broad similarities with the predicted distribution 

presented here; for example, areas of high usage/relative density between Cape Cod and Nantucket, 

between Nantucket and Martha’s Vineyard, and off the east end of Long Island, with lower usage/relative 

density in the eastern central part of the RISAMP area. The winter distribution of Common Eider tagged 

along the coast of Rhode Island (Fig. 2 of Beuth et al. 2017) was quite different from the predicted 

distribution presented here. Areas of high usage by tracked individuals along the coast and near Block 

Island were not evident in our model predictions which were fairly flat across the RISAMP area in 

comparison to localized areas of predicted high relative density further east in and around Nantucket 

Sound. 

Hatch et al. (2016) and Powers et al. (2017) reported on tracking of Great Shearwater in the Gulf of 

Maine from June-December. Broadly speaking, the estimated distribution (Hatch et al. 2016; their Fig. 

1b) and locations (Powers et al. 2017; their Fig. 3) of tracked individuals overlapped with areas of 
23 

predicted highest density presented here for summer and fall. However, Powers et al. (2017) identified a 

migration pathway out of the Gulf of Maine in August and September (their Fig. 4) that is not particularly 

evident in the predicted distributions presented here, highlighting the advantage of tracking data for 

identifying movement corridors. 

Jodice et al. (2015) reported on tracking of Black-capped Petrel tagged at a breeding site in the Caribbean. 

During the post-breeding period (summer-fall) locations of tracked individuals overlapped with sighting 

locations from the at-sea survey data and corresponding areas of predicted highest density presented here, 

for example along the western edge of the Gulf Stream and in the South Atlantic Bight (see Fig. 4 of 

Jodice et al. 2015 for visual comparison). Tracked individuals also visited waters further offshore, east of 

the Gulf Stream and around and beyond the Blake Spur. There was little at-sea survey effort in these 

areas, but the models predicted high relative density of Black-capped Petrel there, especially during 

summer. 

Several past and current studies have tracked Brown Pelican within and adjacent to our study area (Poli

2015; Lamb 2016; Smithsonian Migratory Connectivity Project 

http://www.migratoryconnectivityproject.org/livetracks/). The distribution of tracked individuals along 

the U.S. Atlantic coast (Poli 2015; their Fig. 2.2) exhibited localized areas of high usage that 

corresponded to areas of predicted high relative density in the seasonal distributions presented here (e.g., 

near Charleston, SC; near Cape Canaveral, FL; and along the coast of the Delmarva peninsula and the 

mouth of Chesapeake Bay, especially during the summer breeding season). The distribution of tracked 

individuals presented by Poli (2015) was highly coastal, further confirming the spuriousness of predicted 

areas of moderate density far offshore for this species during winter.

In addition to the studies discussed above, there are other past and current tracking projects of marine 

birds in and adjacent to our study area. For example, at the north end of our study area, alcids have been 

tracked from Machias Seal Island in the Gulf of Maine (e.g., Clarke et al. 2010). The movements of 

Common Tern and Roseate Tern along the northeast U.S. Atlantic coast have been studied using 

geolocators (Nisbet et al. 2011; Mostello et al. 2014) and using VHF transmitters and automated radio 

telemetry stations (Loring 2016). As tracking data continue to accumulate for marine bird species along 

the U.S. Atlantic coast they will provide an invaluable resource for describing and understanding the 

spatial distributions of these species, complementing the information provided by the analysis of at-sea 

survey data (Perrow et al. 2015). We echo the call by Spiegel et al. (2017) for better integration of these 

two types of information, especially the development of a rigorous analytical framework that can be used 

to inform management decisions.

4.7 Conclusion 

The maps presented in this report can be used to inform planning for ocean activities such as offshore 

wind energy development. The maps identify areas where the densities of marine bird species are likely to 

be higher or lower on a seasonal basis. This information can contribute to assessments of relative risk 

exposure aimed at minimizing the impacts of activities on marine birds (Winiarski et al. 2014b; Fifield et 

al. 2017). As discussed in Section 4.2, the survey data analyzed here had some limitations. Any full 

assessments of risk exposure in specific areas at specific times should augment the information provided 

here with more targeted analyses and potentially additional survey (or tracking) data to inform those 

analyses, especially when there are few existing data. Nevertheless, our maps provide a starting point for 

assessing relative risk exposure, which when combined with information about other aspects of risk 

exposure and vulnerability (e.g., flight height, avoidance/attraction, demographic impacts), can inform 
24 

assessments of overall risk for marine birds in U.S. Atlantic OCS waters (Garthe and Hüppop 2004; 

Robinson Willmott et al. 2013; Bailey et al. 2014; Bradbury et al. 2014). 

5. Literature Cited

Bailey, H., Brookes, K. L., and Thompson, P. M. 2014. Assessing environmental impacts of offshore 

wind farms: lessons learned and recommendations for the future. Aquatic Biosystems 10:8. 

Balderama, E., Gardner, B., and Reich, B. J. 2016. A spatial–temporal double-hurdle model for extremely 

over-dispersed avian count data. Spatial Statistics 18:263-275. 

Ballance, L. T., Pitman, R. L., and Fiedler, P. C. 2006. Oceanographic influences on seabirds and 

cetaceans of the eastern tropical Pacific: A review. Progress in Oceanography 69:360-390. 

Barlow, J. 2015. Inferring trackline detection probabilities, g(0), for cetaceans from apparent densities in 

different survey conditions. Marine Mammal Science 31:923‐943.

Bates, D. and Maechler, M. 2017. Matrix: Sparse and dense matrix classes and methods. R package 

version 1.2-8. URL: https://CRAN.R-project.org/package=Matrix (accessed March 1, 2018). 

Becker, R. A., Wilks, A. R., Brownrigg, R., Minka, T. P., and Deckmyn, A. 2016. maps: Draw 

geographical maps. R package version 3.1.1. URL: https://CRAN.R-project.org/package=maps (accessed 

March 1, 2018).

Beuth, J. M., Mcwilliams, S. R., Paton, P. W., and Osenkowski, J. E. 2017. Habitat use and movements of 

common eiders wintering in southern New England. Journal of Wildlife Management 81:1276-1286. 

Bivand, R., Keitt, T., and Rowlingson, B. 2015. rgdal: Bindings for the Geospatial Data Abstraction 

Library. R package version 1.0-4. URL: https://CRAN.R-project.org/package=rgdal (accessed March 1, 

2018).

Bivand, R. and Lewin-Koh, N. 2015. maptools: Tools for reading and handling spatial objects. R package 

version 0.8-36. URL: https://CRAN.R-project.org/package=maptools (accessed March 1, 2018).

Bivand, R. and Rundel, C. 2015. rgeos: Interface to Geometry Engine - Open Source (GEOS). R package 

version 0.3-11. URL: https://CRAN.R-project.org/package=rgeos (accessed March 1, 2018).

Bodenhofer, U., Krone, M., and Klawonn, F. 2013. Testing noisy numerical data for monotonic 

association. Information Sciences 245:21-37. 

Boudt, K., Cornelissen, J., and Croux, C. 2012. The Gaussian rank correlation estimator: robustness 

properties. Statistics and Computing 22:471-483. 

Bradbury, G., Trinder, M., Furness, B., Banks, A. N., Caldow, R. W. G., and Hume, D. 2014. Mapping 

seabird sensitivity to offshore wind farms. PLoSOne 9:1-17. 

Bühlmann, P. and Hothorn, T. 2007. Boosting algorithms: Regularization, prediction and model fitting. 

Statistical Science 22(4):477-505. 
25 

Canty, A. and Ripley, B. 2015. boot: Bootstrap R (S‐Plus) functions. R package version 1.3-17. URL: 

https://cran.r-project.org/web/packages/boot/index.html (accessed March 1, 2018).

Clarke, T. C. R., Diamond, A. W., and Chardine, J. W. 2010. Origin of Canadian Razorbills (Alca torda) 

wintering in the outer Bay of Fundy confirmed by radio-tracking. Waterbirds 33:541-545. 

Conn, P. B., Thorson, J. T., and Johnson, D. S. 2017. Confronting preferential sampling when analysing 

population distributions: diagnosis and model-based triage. Methods in Ecology and Evolution 8:1535-

1546. 

Curtice, C., Cleary J., Shumchenia E., and Halpin P.N. 2016. Marine-life Data and Analysis Team 

(MDAT) technical report on the methods and development of marine-life data to support regional ocean 

planning and management. Prepared on behalf of the Marine-life Data and Analysis Team (MDAT). 

URL: http://seamap.env.duke.edu/models/MDAT/MDAT-Technical-Report-v1_1.pdf (accessed March 1, 

2018). 

Fifield, D. A., Hedd, A., Avery-Gomm, S., Robertson, G. J., Gjerdrum, C., and McFarlane Tranquilla, L. 

2017. Employing predictive spatial models to inform conservation planning for seabirds in the Labrador 

Sea. Frontiers in Marine Science 4:149. 

Flanders, N. P., Gardner, B., Winiarski, K., Paton, P. W. C., Allison, T., and O'Connell, A. F. 2015. Key 

seabird areas in southern New England identified using a community occupancy model. Marine Ecology 

Progress Series 533:277-290. 

Friedman, J. H. 2002. Stochastic gradient boosting. Computational Statistics & Data Analysis 38:367-

378. 

Furrer, R. and Sain, S. R. 2010. spam: A sparse matrix R Package with emphasis on MCMC methods for 

Gaussian Markov random fields. Journal of Statistical Software 36(10).

Garthe, S., and Hüppop, O. 2004. Scaling possible adverse effects of marine wind farms on seabirds: 

developing and applying a vulnerability index. Journal of Applied Ecology 41:724-734. 

Genz, A., Bretz, F., Miwa, T., Mi, X., Leisch, F., Scheipl, F., and Hothorn, T. 2014. mvtnorm: 

Multivariate normal and t distributions. R package version 1.0-2. URL: http://CRAN.Rproject.org/package=mvtnorm (accessed March 1, 2018).

Gjerdrum, C. and Bolduc, F. 2016. Non-breeding distribution of Herring Gull (Larus argentatus) and 

Great Black-backed Gull (Larus marinus) in Eastern Canada from ship-based surveys. Waterbirds 

(Special Publication 1):202-219. 

Goyert, H. F. 2014. Relationship among prey availability, habitat, and the foraging behavior, distribution, 

and abundance of common terns Sterna hirundo and roseate terns S. dougallii. Marine Ecology Progress 

Series 506:291-302. 

Goyert, H. F., Gardner, B., Sollmann, R., Veit, R. R., Gilbert, A. T., Connelly, E. E., and Williams, K. A. 

2016. Predicting the offshore distribution and abundance of marine birds with a hierarchical community 

distance sampling model. Ecological Applications 26:1797-1815. 
26 

Grothendieck, G., Kates, L., and Petzoldt, T. 2012. proto: Prototype object-based programming. R 

package version 0.3-10. URL: https://CRAN.R-project.org/package=proto (accessed March 1, 2018). 

Grothendieck, G. 2014a. gsubfn: Utilities for strings and function arguments. R package version 0.6-6. 

URL: https://CRAN.R-project.org/package=gsubfn (accessed March 1, 2018). 

Grothendieck, G. 2014b. sqldf: Perform SQL selects on R data frames. R package version 0.4-10. URL: 

https://CRAN.R-project.org/package=sqldf (accessed March 1, 2018). 

Hatch, J. M., Wiley, D., Murray, K. T., and Welch, L. 2016. Integrating Satellite-Tagged Seabird and 

Fishery-Dependent Data: A Case Study of Great Shearwaters (Puffinus gravis) and the U.S. New England 

Sink Gillnet Fishery. Conservation Letters 9:43-50. 

Heinänen, S., Žydelis, R., Dorsch, M., Nehls, G., and Skov, H. 2017. High-resolution sea duck 

distribution modeling: Relating aerial and ship survey data to food resources, anthropogenic pressures, 

and topographic variables. Condor 119:175-190. 

Hijmans, R. J. 2015. raster: Geographic data analysis and modeling. R package version 2.4-15. URL: 

https://CRAN.R-project.org/package=raster (accessed March 1, 2018).

Hofner, B., Mayr, A., Robinzonov, N., and Schmid, M. 2012. Model-based boosting in R: A hands-on 

tutorial using the R package mboost. Technical Report 120, Department of Statistics, University of 

Munich, Germany. 

Hofner, B. and Hothorn, T. 2015. stabs: Stability selection with error control. R package version 0.5-1. 

URL: https://CRAN.R-project.org/package=stabs (accessed March 1, 2018).

Hothorn, T., Hornik, K., and Zeileis, A. 2006. Unbiased recursive partitioning: A conditional inference 

framework. Journal of Computational and Graphical Statistics 15:651-674. 

Hothorn, T., Leisch, F., and Zeileis, A. 2013. modeltools: Tools and classes for statistical models. R 

package version 0.2-21. URL: https://CRAN.R-project.org/package=modeltools (accessed March 1, 

2018). 

Hothorn, T., Buehlmann, P., Kneib, T., Schmid, M., and Hofner, B. 2015. mboost: Model-based boosting. 

R package version 2.4-2. URL: https://CRAN.R-project.org/package=mboost (accessed March 1, 2018).

Huettmann, F., Diamond, A. W., Dalzell, B., and Macintosh, K. 2005. Winter distribution, ecology and 

movements of Razorbills Alca torda and other auks in the outer Bay of Fundy, Atlantic Canada. Marine 

Ornithology 33:161-171. 

Jackman, S. 2015. pscl: Classes and methods for R developed in the Political Science Computational 

Laboratory, Department of Political Science, Stanford University. Stanford, CA. R package version 1.4.9. 

URL: https://CRAN.R-project.org/package=pscl (accessed March 1, 2018).

Jenness, J. 2013. DEM Surface Tools. Jenness Enterprises. URL: 

http://www.jennessent.com/arcgis/surface_area.htm (accessed March 1, 2018).

Jodice, P. G. R., Ronconi, R. A., Rupp, E., Wallace, G. E., and Satgé, Y. 2015. First satellite tracks of the 

Endangered black-capped petrel. Endangered Species Research 29:23-33. 
27 

Johnston, A., Thaxter, C. B., Austin, G. E., Cook, A. S., Humphreys, E. M., Still, D. A., Mackay, A., 

Irvine, R., Webb, A., and Burton, N. H. 2015. Modelling the abundance and distribution of marine birds 

accounting for uncertain species identification Journal of Applied Ecology 52:150-160. 

Jonsen, I. D., Flemming, J. M., and Myers, R. A. 2005. Robust state-space modeling of animal movement 

data. Ecology 86:2874-2880. 

Kinlan, B. P., Winship, A. J., White, T. P., and Christensen, J. 2016. Modeling at-sea occurrence and 

abundance of marine birds to support Atlantic marine renewable energy planning: Phase I report. U.S. 

Department of the Interior, Bureau of Ocean Energy Management, OCS Study BOEM 2016-039. 

Sterling, VA. 113 pp. 

Lamb, J. S. 2016. Ecological drivers of brown pelican movement patterns and reproductive success in the 

Gulf of Mexico. Ph.D. thesis. Clemson University, Clemson, SC. 216 pp. 

Lieske, D. J., Fifield, D. A., and Gjerdrum, C. 2014. Maps, models, and marine vulnerability: Assessing 

the community distribution of seabirds at-sea. Biological Conservation 172:15-28. 

Loring, P. H., Paton, P. W., Osenkowski, J. E., Gilliland, S. G., Savard, J.-P. L., and Mcwilliams, S. R. 

2014. Habitat use and selection of black scoters in southern New England and siting of offshore wind 

energy facilities. Journal of Wildlife Management 78:645-656. 

Loring, P. H. 2016. Evaluating digital VHF technology to monitor shorebird and seabird use of offshore 

wind energy areas in the western North Atlantic. Ph.D. thesis. University of Massachusetts, Amherst, 

MA. 172 pp. 

Mannocci, L., Boustany, A. M., Roberts, J. J., Palacios, D. M., Dunn, D. C., Halpin, P. N., Viehman, S., 

Moxley, J., Cleary, J., Bailey, H., Bograd, S. J., Becker, E. A., Gardner, B., Hartog, J. R., Hazen, E. L., 

Ferguson, M. C., Forney, K. A., Kinlan, B. P., Oliver, M. J., Perretti, C. T., Ridoux, V., Teo, S. L. H., and 

Winship, A. J. 2017. Temporal resolutions in species distribution models of highly mobile marine 

animals: Recommendations for ecologists and managers Diversity and Distributions 23:1098-1109. 

Mayr, A., Fenske, N., Hofner, B., Kneib, T., and Schmid, M. 2012. Generalized additive models for 

location, scale and shape for high dimensional data---a flexible approach based on boosting. Applied 

Statistics 61(3):403-427. 

Miller, P. I. and Christodoulou, S. 2014. Frequent locations of oceanic fronts as an indicator of pelagic 

diversity: Application to marine protected areas and renewables. Marine Policy 45:318-329. 

Mostello, C. S., Nisbet, I. C. T., Oswald, S. A., and Fox, J. W. 2014. Non-breeding season movements of 

six North American Roseate Terns Sterna dougallii tracked with geolocators. Seabird 27:1-21. 

Müller, K., Wickham, H., James, D. A., and Falcon, S. 2014. RSQLite: 'SQLite' interface for R. R 

package version 1.0.0. URL: https://CRAN.R-project.org/package=RSQLite (accessed March 1, 2018).

Nisbet, I. C. T., Mostello, C. S., Veit, R. R., Fox, J. W., and Afanasyev, V. 2011. Migrations and winter 

quarters of five Common Terns tracked using geolocators. Waterbirds 34:32-39. 
28 

Nisbet, I. C., Veit, R. R., Auer, S. A., and White, T. P. 2013. Marine birds of the eastern United States 

and the Bay of Fundy: Distribution, numbers, trends, threats, and management. Nuttall Ornithological 

Club, Cambridge, MA. 

Nisbet, I. C. T. and Veit, R. R. 2015. An explanation for the population crash of Red-necked Phalaropes 

Phalaropus lobatus staging in the Bay of Fundy in the 1980s. Marine Ornithology 43:119-121. 

Nychka, D., Furrer, R., Paige, J., and Sain, S. 2015. fields: Tools for spatial data. R package version 8.10. 

URL: https://CRAN.R-project.org/package=fields (accessed March 1, 2018).

O’Connell, A. F., Gardner, B., Gilbert, A. T., and Laurent, K. 2009. Compendium of Avian Occurrence 

Information for the Continental Shelf Waters along the Atlantic Coast of the United States, Final report

(Database Section – Seabirds). U.S. Department of the Interior, Bureau of Ocean Energy Management, 

OCS Study BOEM 2012-076. Herndon, VA. 344 pp. 

Oedekoven, C. S., Ainley, D. G., and Spear, L. B. 2001. Variable responses of seabirds to change in 

marine climate: California Current, 1985-1994. Marine Ecology Progress Series 212:265-281. 

Pebesma, E. J., and Bivand, R. S. 2005. Classes and methods for spatial data in R. R News 5(2). URL: 

https://CRAN.R-project.org/doc/Rnews/ (accessed March 1, 2018). 

Perrow, M. R., Harwood, A. J., Skeate, E. R., Praca, E., and Eglington, S. M. 2015. Use of multiple data 

sources and analytical approaches to derive a marine protected area for a breeding seabird. Biological 

Conservation 191:729 – 738. 

Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., and R Core Team. 2017. nlme: Linear and nonlinear 

mixed effects models. R package version 3.1-131. URL: http://CRAN.R-project.org/package=nlme 

(accessed March 1, 2018). 

Poli, C. 2015. Variability in movement patterns and habitat use of two species of Pelecaniformes. M.Sc. 

thesis. Clemson University, Clemson, SC. 104 pp. 

Powers, K. D., Wiley, D. N., Allyn, A. J., Welch, L. J., and Ronconi, R. A. 2017. Movements and 

foraging habitats of great shearwaters Puffinus gravis in the Gulf of Maine. Marine Ecology Progress 

Series 574:211-226. 

R Core Team. 2017. R: A language and environment for statistical computing. R Foundation for 

Statistical Computing, Vienna, Austria. URL: http://www.R-project.org/ (accessed March 1, 2018). 

R Special Interest Group on Databases (R-SIG-DB), Wickham, H., and Müller, K. 2017. DBI: R 

Database Interface. R package version 0.6-1. URL: https://CRAN.R-project.org/package=DBI (accessed 

March 1, 2018). 

Renner, M., Parrish, J. K., Piatt, J. F., Kuletz, K. J., Edwards, A. E., and Hunt Jr, G. L. 2013. Modeled 

distribution and abundance of a pelagic seabird reveal trends in relation to fisheries. Marine Ecology 

Progress Series 484:259-277. 

Robin, X., Turck, N., Hainard, A., Tiberti, N., Lisacek, F., Sanchez, J.-C., and Müller, M. 2011. pROC: 

an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics 12:77.
29 

Robinson Willmott, J. C., G. Forcey, and A. Kent. 2013. The relative vulnerability of migratory bird 

species to offshore wind energy projects on the Atlantic Outer Continental Shelf: An assessment method 

and database. Final Report to the U.S. Department of the Interior, Bureau of Ocean Energy Management, 

Office of Renewable Energy Programs. OCS Study BOEM 2013-207. 275 pp. 

Ryan, W. B. F., Carbotte, S. M., Coplan, J. O., O'Hara, S., Melkonian, A., Arko, R., Weissel, R.A., 

Ferrini, V., Goodwillie, A., Nitsche, F., Bonczkowski, J., and Zemsky, R. 2009. Global Multi-Resolution 

Topography synthesis. Geochemistry, Geophysics, Geosystems 10:Q03014. 

Sarkar, D. 2008. Lattice: Multivariate data visualization with R. Springer, New York, NY.

Schmid, M., Potapov, S., Pfahlberg, A., and Hothorn, T. 2008. Estimation and regularization techniques 

for regression models with multidimensional prediction functions. Technical Report 42, Department of 

Statistics, University of Munich, Germany.

Scott-Hayward, L. A. S., Mackenzie, M. L., Donovan, C. R., Walker, C. G., and Ashe, E. 2014. Complex 

Region Spatial Smoother (CReSS). Journal of Computational and Graphical Statistics 23:340-360. 

Sollmann, R., Gardner, B., Williams, K. A., Gilbert, A. T., and Veit, R. R. 2016. A hierarchical distance 

sampling model to estimate abundance and covariate associations of species and communities. Methods 

in Ecology and Evolution 7:529-537. 

Spear, L.B., D.G. Ainley, and P. Pyle. 1999. Seabirds in southeastern Hawaiian waters. Western Birds 

30:1‐32.

Spiegel, C. S., Berlin, A. M., Gilbert, A. T., Gray, C. O., Montevecchi, W. A., Stenhouse, I. J., Ford, S. 

L., Olsen, G. H., Fiely, J. L., Savoy, L., Goodale, M. W., and Burke, C. M. 2017. Determining fine-scale 

use and movement patterns of diving bird species in federal waters of the mid-Atlantic United States 

using satellite telemetry. U.S. Department of the Interior, Bureau of Ocean Energy Management, OCS 

Study BOEM 2017-069. Sterling, VA. 260 pp. 

Suryan, R. M., Kuletz, K. J., Parker-Stetter, S. L., Ressler, P. H., Renner, M., Horne, J. K., Farley, E. V., 

and Labunski, E. A. 2016. Temporal shifts in seabird populations and spatial coherence with prey in the 

southeastern Bering Sea. Marine Ecology Progress Series 549:199-215. 

Tasker, M. L., Jones, P. H., Dixon, T., and Blake, B. F. 1984. Counting seabirds at sea from ships: a 

review of methods employed and a suggestion for a standardized approach. Auk 101:567-577. 

Tremblay, Y., Bertrand, S., Henry, R. W., Kappes, M. A., Costa, D. P., and Shaffer, S. A. 2009. 

Analytical approaches to investigating seabird-environment interactions: a review. Marine Ecology 

Progress Series 391:153-163. 

Veit, R. R., White, T. P., Perkins, S. A., and Curely, S. 2016. Abundance and distribution of seabirds off 

southeastern Massachusetts, 2011-2015. U.S. Department of the Interior, Bureau of Ocean Energy 

Management, OCS Study BOEM 2016-067. Sterling, VA. 82 pp.

Venables, W. N., and Ripley, B. D. 2002. Modern Applied Statistics with S. Fourth Edition. Springer, 

New York, NY.
30 

Wakefield, E. D., Phillips, R. A., and Matthiopoulos, J. 2009. Quantifying habitat use and preferences of 

pelagic seabirds using individual movement data: a review. Marine Ecology Progress Series 391:165-182. 

Wakefield, E. D., Owen, E., Baer, J., Carroll, M. J., Daunt, F., Dodd, S. G., Green, J. A., Guilford, T., 

Mavor, R. A., Miller, P. I., Newell, M. A., Newton, S. F., Robertson, G. S., Shoji, A., Soanes, L. M., 

Votier, S. C., Wanless, S., and Bolton, M. 2017. Breeding density, fine-scale tracking, and large-scale 

modeling reveal the regional distribution of four seabird species. Ecological Applications 27:2074-2091. 

Weatherall, P., Marks, K. M., Jakobsson, M., Schmitt, T., Tani, S., Arndt, J. E., Rovere, M., Chayes, D., 

Ferrini, V., and Wigley, R. 2015. A new digital bathymetric model of the world's oceans. Earth and Space 

Science 2:331-345. 

White, T. P., Veit, R. R., and Perry, M. C. 2009. Feeding ecology of Long-tailed Ducks Clangula 

hyemalis wintering on the Nantucket Shoals. Waterbirds 32:293-299. 

White, T. P. and Veit. R. R. 2018. Spatial ecology of long-tailed ducks (Clangula hyemalis) and whitewinged scoters (Melanitta fusca) wintering on Nantucket Shoals: persistent foraging areas stabilized by 

local enhancement. Unpublished manuscript. 

Wickham, H. 2007. Reshaping data with the reshape package. Journal of Statistical Software 21(12).

Winiarski, K. J., Miller, D. L., Paton, P. W. C., and McWilliams, S. R. 2013. Spatially explicit model of 

wintering common loons: conservation implications. Marine Ecology Progress Series 492:273-283. 

Winiarski, K. J., Burt, M. L., Rexstad, E., Miller, D. L., Trocki, C. L., Paton, P. W. C., and McWilliams, 

S. R. 2014a. Integrating aerial and ship surveys of marine birds into a combined density surface model: A 

case study of wintering Common Loons. Condor 116:149-161. 

Winiarski, K. J., Miller, D. L., Paton, P. W., and McWilliams, S. R. 2014b. A spatial conservation 

prioritization approach for protecting marine birds given proposed offshore wind energy development. 

Biological Conservation 169:79-88. 

Wood, S.N. 2006. Generalized additive models: An introduction with R. Chapman and Hall/CRC. Boca 

Raton, FL. 392 pp. 

Yee, T. W. 2015. Vector generalized linear and additive models: With an implementation in R. Springer, 

New York, NY.

Zeileis, A., Leisch, F., Hornik, K., and Kleiber, C. 2002. strucchange: An R package for testing for 

structural change in linear regression models. Journal of Statistical Software 7(2).

Zeileis, A. and Grothendieck, G. 2005. zoo: S3 infrastructure for regular and irregular time series. Journal 

of Statistical Software 14(6).

Zeileis, A. 2006. Object-oriented computation of sandwich estimators. Journal of Statistical Software 

16(9). 

Zeileis, A., Kleiber, C., and Jackman, S. 2008. Regression models for count data in R. Journal of 

Statistical Software 27(8).
31 

Zipkin, E. F., Gardner, B., Gilbert, A. T., O'Connell, A. F., Royle, J. A., and Silverman, E. D. 2010. 

Distribution patterns of wintering sea ducks in relation to the North Atlantic Oscillation and local 

environmental characteristics. Oecologia 163:893-902. 

Zipkin, E. F., Leirness, J. B., Kinlan, B. P., O’Connell, A. F., and Silverman, E. D. 2014. Fitting 

statistical distributions to sea duck count data: Implications for survey design and abundance estimation. 

Statistical Methodology 17:67-81. 

Zipkin, E. F., Kinlan, B. P., Sussman, A., Rypkema, D., Wimer, M., and O'Connell, A. F. 2015. Statistical 

guidelines for assessing marine avian hotspots and coldspots: A case study on wind energy development 

in the U.S. Atlantic Ocean. Biological Conservation 191:216-223. 
32 

Table 1. Survey datasets, sample sizes, and surveyed area analyzed. All datasets (and their ID codes) were extracted from the USFWS Northwest Atlantic 

Seabird Catalog except ECSAS which was provided by CWS-ECCC. Continuous data were segmented as described in Section 2.2. A description of each 

individual survey is presented in Appendix A. Area surveyed does not necessarily represent unique area surveyed; some surveys covered similar areas 

over multiple days, months, or years. Seasons were defined as spring (March-May), summer (June-August), fall (September-November), and winter 

(December-February). Sample sizes and surveyed area analyzed for corrected MassCEC datasets (see Section 2.2) are shown in parentheses.

Source Dataset ID Platform Method

Year Number of transect segments Area surveyed (km2

) 

Start End Spring Summer Fall Winter Total Spring Summer Fall Winter Total

AMAPPS_FWS_Aerial_Fall2012 aerial continuous 2012 2012 0 0 2,986 0 2,986 0 0 4,765 0 4,765

AMAPPS_FWS_Aerial_Fall2013 aerial continuous 2013 2013 0 0 4,629 0 4,629 0 0 7,395 0 7,395

AMAPPS_FWS_Aerial_Fall2014 aerial continuous 2014 2014 0 0 2,876 0 2,876 0 0 4,608 0 4,608

AMAPPS_FWS_Aerial_Preliminary_

Summer2010 aerial continuous 2010 2010 0 1,131 0 0 1,131 0 1,802 0 0 1,802

AMAPPS_FWS_Aerial_Spring2012 aerial continuous 2012 2012 2,962 0 0 0 2,962 4,739 0 0 0 4,739

AMAPPS_FWS_Aerial_Summer2011 aerial continuous 2011 2011 0 3,442 0 0 3,442 0 5,502 0 0 5,502

AMAPPS_FWS_Aerial_Winter2010-2011 aerial continuous 2010 2011 0 0 0 513 513 0 0 0 823 823

AMAPPS_FWS_Aerial_Winter2014 aerial continuous 2014 2014 0 0 0 3,073 3,073 0 0 0 4,914 4,914

AMAPPS_NOAA/NMFS_NEFSCBoat2011 boat continuous 2011 2011 0 1,537 0 0 1,537 0 1,794 0 0 1,794

AMAPPS_NOAA/NMFS_NEFSCBoat2013 boat continuous 2013 2013 0 1,577 0 0 1,577 0 1,853 0 0 1,853

AMAPPS_NOAA/NMFS_NEFSCBoat2014 boat continuous 2014 2014 1,023 0 0 0 1,023 1,219 0 0 0 1,219

AMAPPS_NOAA/NMFS_NEFSCBoat2015 boat continuous 2015 2015 0 261 0 0 261 0 308 0 0 308

AMAPPS_NOAA/NMFS_SEFSCBoat2011 boat continuous 2011 2011 0 982 0 0 982 0 1,155 0 0 1,155

AMAPPS_NOAA/NMFS_SEFSCBoat2013 boat continuous 2013 2013 0 701 277 0 978 0 823 326 0 1,149

BarHarborWW05 boat continuous 2005 2005 0 876 181 0 1,057 0 1,048 217 0 1,265

BarHarborWW06 boat continuous 2006 2006 0 821 331 0 1,152 0 994 399 0 1,393

CapeHatteras0405 boat continuous 2004 2005 0 195 0 168 363 0 213 0 161 374

CapeWindAerial aerial continuous 2002 2004 1,117 1,112 1,176 1,271 4,676 1,788 1,782 1,881 2,041 7,492

CapeWindBoat boat continuous 2002 2003 117 109 29 0 255 748 712 184 0 1,644

CDASMidAtlantic aerial continuous 2001 2003 241 0 0 1,363 1,604 113 0 0 653 766

CSAP boat discrete 1980 1988 7,612 6,975 7,325 4,213 26,125 9,712 8,967 9,466 5,400 33,545

DOEBRIAerial2012 camera continuous 2012 2012 1,486 832 1,515 763 4,596 1,190 661 1,208 610 3,669

DOEBRIAerial2013 camera continuous 2013 2013 941 1,615 1,479 1,265 5,300 765 1,297 1,178 1,010 4,250

DOEBRIAerial2014 camera continuous 2014 2014 1,029 0 0 1,341 2,370 822 0 0 1,074 1,896
33 

Source Dataset ID Platform Method

Year Number of transect segments Area surveyed (km2

) 

Start End Spring Summer Fall Winter Total Spring Summer Fall Winter Total

DOEBRIBoatApr2014 boat continuous 2014 2014 164 0 0 0 164 195 0 0 0 195

DOEBRIBoatApril2012 boat continuous 2012 2012 165 0 0 0 165 197 0 0 0 197

DOEBRIBoatAug2012 boat continuous 2012 2012 0 164 0 0 164 0 197 0 0 197

DOEBRIBoatAug2013 boat continuous 2013 2013 0 166 0 0 166 0 199 0 0 199

DOEBRIBoatDec2012 boat continuous 2012 2013 0 0 0 162 162 0 0 0 194 194

DOEBRIBoatDec2013 boat continuous 2013 2013 0 0 0 170 170 0 0 0 202 202

DOEBRIBoatJan2013 boat continuous 2013 2013 0 0 0 164 164 0 0 0 198 198

DOEBRIBoatJan2014 boat continuous 2014 2014 0 0 0 164 164 0 0 0 197 197

DOEBRIBoatJune2012 boat continuous 2012 2012 0 166 0 0 166 0 200 0 0 200

DOEBRIBoatJune2013 boat continuous 2013 2013 0 168 0 0 168 0 200 0 0 200

DOEBRIBoatMar2013 boat continuous 2013 2013 166 0 0 0 166 201 0 0 0 201

DOEBRIBoatMay2013 boat continuous 2013 2013 168 0 0 0 168 201 0 0 0 201

DOEBRIBoatNov2012 boat continuous 2012 2012 0 0 165 0 165 0 0 197 0 197

DOEBRIBoatOct2013 boat continuous 2013 2013 0 0 170 0 170 0 0 201 0 201

DOEBRIBoatSep2012 boat continuous 2012 2012 0 0 168 0 168 0 0 201 0 201

DOEBRIBoatSep2013 boat continuous 2013 2013 0 0 168 0 168 0 0 201 0 201

DominionVirginia_VOWTAP boat continuous 2013 2014 24 18 12 24 78 77 58 38 77 250

EcoMonAug08 boat continuous 2008 2008 0 480 0 0 480 0 575 0 0 575

EcoMonAug09 boat continuous 2009 2009 0 458 0 0 458 0 547 0 0 547

EcoMonAug10 boat continuous 2010 2010 0 480 12 0 492 0 573 15 0 588

EcoMonAug2012 boat continuous 2012 2012 0 656 0 0 656 0 782 0 0 782

EcoMonFeb10 boat continuous 2010 2010 0 0 0 334 334 0 0 0 398 398

EcoMonFeb2012 boat continuous 2012 2012 0 0 0 549 549 0 0 0 661 661

EcoMonFeb2013 boat continuous 2013 2013 0 0 0 521 521 0 0 0 620 620

EcoMonJan09 boat continuous 2009 2009 0 0 0 391 391 0 0 0 474 474

EcoMonJun2012 boat continuous 2012 2012 31 513 0 0 544 38 613 0 0 651

EcoMonMay07 boat continuous 2007 2007 433 72 0 0 505 520 86 0 0 606

EcoMonMay09 boat continuous 2009 2009 199 422 0 0 621 241 505 0 0 746

EcoMonMay10 boat continuous 2010 2010 274 370 0 0 644 328 442 0 0 770
34 

Source Dataset ID Platform Method

Year Number of transect segments Area surveyed (km2

) 

Start End Spring Summer Fall Winter Total Spring Summer Fall Winter Total

EcoMonNov09 boat continuous 2009 2009 0 0 441 0 441 0 0 528 0 528

EcoMonNov10 boat continuous 2010 2010 0 0 418 0 418 0 0 500 0 500

EcoMonNov2011 boat continuous 2011 2011 0 0 454 0 454 0 0 542 0 542

EcoMonOct2012 boat continuous 2012 2012 0 0 498 0 498 0 0 598 0 598

ECSAS boat discrete 2006 2016 5,993 2,810 3,700 513 13,016 3,068 1,425 2,004 230 6,727

FLPowerLongIsland_Aerial aerial continuous 2004 2006 104 49 158 0 311 154 75 237 0 466

FLPowerLongIsland_Boat boat continuous 2004 2006 497 159 346 211 1,213 566 188 385 235 1,374

FWS_MidAtlanticDetection_Spring2012 aerial continuous 2012 2012 177 0 0 0 177 283 0 0 0 283

FWS_SouthernBLSC_Winter2012 aerial continuous 2012 2012 0 0 0 904 904 0 0 0 1,500 1,500

FWSAtlanticWinterSeaduck2008 aerial continuous 2008 2011 78 0 0 8,311 8,389 124 0 0 13,295 13,419

GeorgiaPelagic boat discrete 1982 1985 681 698 576 231 2,186 844 727 705 293 2,569

HatterasEddyCruise2004 boat continuous 2004 2004 0 131 0 0 131 0 117 0 0 117

HerringAcoustic06 boat continuous 2006 2006 0 0 287 0 287 0 0 341 0 341

HerringAcoustic07 boat continuous 2007 2007 0 0 334 0 334 0 0 395 0 395

HerringAcoustic08 boat continuous 2008 2008 0 0 822 0 822 0 0 990 0 990

HerringAcoustic09Leg1 boat continuous 2009 2009 0 0 127 0 127 0 0 151 0 151

HerringAcoustic09Leg2 boat continuous 2009 2009 0 0 289 0 289 0 0 341 0 341

HerringAcoustic09Leg3 boat continuous 2009 2009 0 0 263 0 263 0 0 315 0 315

HerringAcoustic2010 boat continuous 2010 2010 0 0 555 0 555 0 0 670 0 670

HerringAcoustic2011 boat continuous 2011 2011 0 0 808 0 808 0 0 950 0 950

HerringAcoustic2012 boat continuous 2012 2012 0 0 772 0 772 0 0 917 0 917

MassAudNanAerial aerial continuous 2002 2006 860 988 1,374 2,004 5,226 631 715 1,004 1,464 3,814

MassCEC2011-2012 aerial continuous 2011 2012 668

(670)

503

(503)

845

(846)

495

(495)

2,511

(2,514)

1,068

(1,072)

804

(804)

1,354

(1,355)

790

(790)

4,016

(4,021)

MassCEC2013 aerial continuous 2013 2013 498

(498)

499

(499)

921

(921)

330

(330)

2,248

(2,248)

797

(797)

798

(798)

1,473

(1,473)

528

(528)

3,596

(3,596)

MassCEC2014 aerial continuous 2014 2015 167

(167)

500

(500)

171

(171)

674

(674)

1,512

(1512)

267

(267)

801

(801)

273

(273)

1,080

(1,080)

2,421

(2,421)

NewEnglandSeamount06 boat discrete 2007 2007 61 4 0 0 65 34 2 0 0 36

NJDEP2009 boat continuous 2008 2009 1,224 1,350 1,430 967 4,971 1,460 1,618 1,727 1,162 5,967

NOAA/NMFS_NEFSCBoat2004 boat continuous 2004 2004 0 1,207 0 0 1,207 0 1,422 0 0 1,422
35 

Source Dataset ID Platform Method

Year Number of transect segments Area surveyed (km2

) 

Start End Spring Summer Fall Winter Total Spring Summer Fall Winter Total

NOAA/NMFS_NEFSCBoat2007 boat continuous 2007 2007 0 633 0 0 633 0 746 0 0 746

NOAAMBO7880 boat discrete 1978 1979 1,677 2,458 2,043 787 6,965 1,597 2,360 1,779 681 6,417

PlattsBankAerial aerial continuous 2005 2005 0 869 0 0 869 0 1,178 0 0 1,178

RISAMPAerial aerial continuous 2009 2010 985 767 0 714 2,466 1,180 918 0 855 2,953

RISAMPBoat boat continuous 2009 2010 202 288 128 163 781 276 354 173 219 1,022

SEFSC1992 boat continuous 1992 1992 0 0 0 783 783 0 0 0 938 938

SEFSC1998 boat continuous 1998 1998 0 1,365 0 0 1,365 0 1,596 0 0 1,596

SEFSC1999 boat continuous 1999 1999 0 730 524 0 1,254 0 852 623 0 1,475

StatoilMaine boat continuous 2012 2013 40 150 140 70 400 48 180 168 84 480

WHOIJuly2010 boat continuous 2010 2010 0 86 0 0 86 0 102 0 0 102

WHOISept2010 boat continuous 2010 2010 0 0 85 0 85 0 0 99 0 99

Total 1978 2016 32,064 42,543 42,008 33,606 150,221 35,491 50,866 51,722 43,061 181,140
36 

Table 2. List of species-season combinations modeled with sample sizes and counts. Only combinations with ≥50 sightings were modeled. Seasons were 

defined as spring (March-May), summer (June-August), fall (September-November), and winter (December-February). 

Species 

code Common name Scientific name Family

Number of segments Count

Spring

Summer

Fall

Winter

Total

Spring

Summer

Fall

Winter

Total

ARTE Arctic Tern Sterna paradisaea Sternidae 170 170 524 524

ATPU Atlantic Puffin Fratercula arctica Alcidae 362 287 124 342 1,115 723 639 174 589 2,125

AUSH Audubon's Shearwater Puffinus lherminieri Procellariidae 134 916 297 170 1,517 449 2,709 982 340 4,480

BCPE Black-capped Petrel Pterodroma hasitata Procellariidae 159 371 93 90 713 315 1,001 246 212 1,774

BLGU Black Guillemot Cepphus grylle Alcidae 90 90 167 167

BLKI Black-legged Kittiwake Rissa tridactyla Laridae 741 2,239 4,066 7,046 3,426 15,106 34,357 52,889

BLSC1 Black Scoter Melanitta americana Anatidae 516 443 1,330 2,289 19,582 10,224 64,882 94,688

BOGU Bonaparte's Gull Chroicocephalus philadelphia Laridae 467 329 1,585 2,381 5,808 2,297 14,072 22,177

BRPE Brown Pelican Pelecanus occidentalis Pelecanidae 66 127 164 76 433 307 567 815 316 2,005

BRSP2 Band-rumped Storm-Petrel Oceanodroma castro Hydrobatidae 276 276 558 558

BRTE Bridled Tern Sterna anaethetus Sternidae 101 65 166 224 148 372

COEI1 Common Eider Somateria mollissima Anatidae 906 146 640 2,172 3,864 208,330 23,280 63,016 536,264 830,890

COLO1 Common Loon Gavia immer Gaviidae 2,932 211 1,472 3,825 8,440 7,086 283 3,199 9,770 20,338

COMU Common Murre Uria aalge Alcidae 212 268 480 520 759 1,279

COSH1 Cory's Shearwater Calonectris diomedea Procellariidae 137 3,383 1,944 5,464 272 15,084 9,193 24,549

COTE1 Common Tern Sterna hirundo Sternidae 642 1,806 777 3,225 2,934 8,315 7,753 19,002

DCCO Double-crested Cormorant Phalacrocorax auritus Phalacrocoracidae 158 187 278 157 780 4,282 1,430 5,710 2,134 13,556

DOVE Dovekie Alle alle Alcidae 664 61 468 1,252 2,445 3,774 126 4,058 9,627 17,585

GBBG Great Black-backed Gull Larus marinus Laridae 3,882 3,513 6,155 3,902 17,452 27,567 12,000 37,287 32,058 108,912

GRSH Great Shearwater Puffinus gravis Procellariidae 682 7,351 7,531 140 15,704 6,047 186,133 102,183 591 294,954

GRSK Great Skua Stercorarius skua Stercorariidae 196 196 223 223

HERG Herring Gull Larus argentatus Laridae 6,384 3,202 8,612 5,300 23,498 55,681 11,147 63,122 32,945 162,895

HOGR Horned Grebe Podiceps auritus Podicipedidae 103 103 200 200

LAGU Laughing Gull Leucophaeus atricilla Laridae 742 1,750 1,871 134 4,497 1,767 5,769 10,339 354 18,229
37 

Species 

code Common name Scientific name Family

Number of segments Count

Spring

Summer

Fall

Winter

Total

Spring

Summer

Fall

Winter

Total

LESP Leach's Storm-Petrel Oceanodroma leucorhoa Hydrobatidae 279 2,540 591 3,410 920 10,336 1,496 12,752

LETE Least Tern Sterna antillarum Sternidae 126 98 224 457 1,202 1,659

LTDU1 Long-tailed Duck Clangula hyemalis Anatidae 1,335 539 3,564 5,438 90,219 18,805 151,955 260,979

MASH Manx Shearwater Puffinus puffinus Procellariidae 107 353 344 804 165 944 595 1,704

NOFU Northern Fulmar Fulmarus glacialis Procellariidae 2,678 956 2,171 2,061 7,866 31,524 13,376 10,666 24,805 80,371

NOGA Northern Gannet Morus bassanus Sulidae 6,729 1,358 5,532 7,932 21,551 41,213 2,951 24,453 61,036 129,653

PAJA Parasitic Jaeger Stercorarius parasiticus Stercorariidae 53 77 191 321 59 103 236 398

POJA Pomarine Jaeger Stercorarius pomarinus Stercorariidae 112 155 830 1,097 144 186 1,169 1,499

RAZO1 Razorbill Alca torda Alcidae 1,063 87 194 1,990 3,334 6,979 230 1,309 16,772 25,290

RBGU Ring-billed Gull Larus delawarensis Laridae 211 53 414 745 1,423 451 145 1,532 3,357 5,485

RBME Red-breasted Merganser Mergus serrator Anatidae 69 112 181 379 950 1,329

REPH Red Phalarope Phalaropus fulicarius Scolopacidae 480 250 338 1,068 86,165 26,657 2,597 115,419

RNPH Red-necked Phalarope Phalaropus lobatus Scolopacidae 143 182 201 526 2,452 2,354 1,372 6,178

ROST1 Roseate Tern Sterna dougallii Sternidae 59 212 83 354 203 805 533 1,541

ROYT Royal Tern Sterna maxima Sternidae 270 289 352 911 752 648 937 2,337

RTLO1 Red-throated Loon Gavia stellata Gaviidae 1,746 412 2,607 4,765 5,219 1,608 8,484 15,311

SOSH Sooty Shearwater Puffinus griseus Procellariidae 916 1,812 119 2,847 6,671 36,282 299 43,252

SOTE Sooty Tern Sterna fuscata Sternidae 60 119 179 617 595 1,212

SPSK South Polar Skua Stercorarius maccormicki Stercorariidae 92 142 234 108 191 299

SUSC1 Surf Scoter Melanitta perspicillata Anatidae 846 848 1,918 3,612 17,211 29,378 51,116 97,705

TBMU Thick-billed Murre Uria lomvia Alcidae 315 151 466 1,324 422 1,746

WISP1 Wilson's Storm-Petrel Oceanites oceanicus Hydrobatidae 1,750 9,269 1,481 12,500 15,059 100,317 7,820 123,196

WWSC1 White-winged Scoter Melanitta fusca Anatidae 568 667 1,536 2,771 21,676 12,241 30,142 64,059

Total 678,272 466,450 454,514 1,088,509 2,687,745

1 Models were re-run with the corrected MassCEC datasets (see Section 2.2)

2 The American Ornithological Union four-letter species code is BSTP
38 

Table 3. Predictor variables used in models with native resolutions, year ranges, and sources.

Predictor variable Native spatial 

resolution

Native temporal 

resolution Years Source

Survey

survey platform survey data

survey ID survey data

transect ID survey data

Temporal

year yearly

day of year daily

Atlantic Multidecadal 

Oscillation (AMO) index1 monthly 1978-2016

NOAA ESRL 

http://www.esrl.noaa.gov/psd/dat

a/correlation/amon.us.data

Multivariate El NiñoSouthern Oscillation Index 

(MEI)1 

monthly 1978-2016

NOAA ESRL 

http://www.esrl.noaa.gov/psd/dat

a/correlation/mei.data

North Atlantic Oscillation 

(NAO) index1 monthly 1978-2016

NOAA ESRL 

http://www.esrl.noaa.gov/psd/dat

a/correlation/nao.data

Trans-Niño Index (TNI)1 monthly 1978-2016

NOAA ESRL 

http://www.esrl.noaa.gov/psd/dat

a/correlation/tni.data

Geographic

longitude projected2

latitude projected2

distance to land3 50 m

derived from Prototype Global 

Shoreline Data

(http://msi.nga.mil/NGAPortal/DN

C.portal?_nfpb=true&_pageLabel

=dnc_portal_page_72)

Bathymetric
39 

Predictor variable Native spatial 

resolution

Native temporal 

resolution Years Source

depth4 

Okeanos Explorer: 25 m;

ACUMEN: approx. 25 m;

CCOM: approx. 100 m; 

GMRT: approx. 400 m;

CRM: approx. 100 m;

GEBCO: approx. 1 km

Okeanos Explorer 

(http://oceanexplorer.noaa.gov/o

keanos/welcome.html) 

obtainable from NOAA National 

Centers for Environmental 

Information Bathymetry Data 

Viewer 

(https://maps.ngdc.noaa.gov/vie

wers/bathymetry/); 

NOAA Atlantic Canyons 

Undersea Mapping Expeditions 

(ACUMEN) 

(http://oceanexplorer.noaa.gov/o

keanos/explorations/acumen12/s

ummary/welcome.html); 

University of New Hampshire 

Center for Coastal and Ocean 

Mapping (CCOM) / Joint 

Hydrographic Center

(http://ccom.unh.edu/theme/lawsea/law-of-the-sea-data/atlantic

and 

https://www.ngdc.noaa.gov/ships

/marcus_g_langseth/MGL1512_

mb.html); 

Global Multi-Resolution 

Topography Data Synthesis 

(GMRT) (http://www.marinegeo.org/portals/gmrt/; Ryan et al. 

2009);

NOAA National Geophysical 

Data Center U.S. Coastal Relief 

Model (CRM) 

(http://www.ngdc.noaa.gov/mgg/c

oastal/crm.html); 

General Bathymetric Chart of the 

Oceans (GEBCO) GEBCO_2014

grid version 20150318

(http://www.gebco.net/data_and_

products/gridded_bathymetry_da

ta/gebco_30_second_grid/; 

Weatherall et al. 2015)

slope (2 km) derived from depth5,6

slope (10 km) derived from depth5,6

slope of slope (10 km) derived from depth5,6

planform curvature (10 km) derived from depth6,7

profile curvature (10 km) derived from depth6,7

Oceanographic (seasonal 

climatologies)

surface chlorophyll-a 

(reflectance inputs between 

440 and 670 nm) 

approx. 4 km daily 1997-2017

NASA OceanColor 

(https://oceancolor.gsfc.nasa.gov

/); SeaWiFS, MODIS Aqua, and

VIIRS sensors/missions

turbidity (reflectance at 547

nm) approx. 4 km daily 2000-2017

NASA OceanColor 

(https://oceancolor.gsfc.nasa.gov

/); MODIS Aqua and MODIS 

Terra sensors/missions
40 

Predictor variable Native spatial 

resolution

Native temporal 

resolution Years Source

surface current velocity (u 

direction) approx. 9 km daily 1992-2012

Hybrid vertical coordinate system 

ocean model (HYCOM) global 

reanalysis (GLBu0.08 grid, 

experiments 19.0 and 19.1), 

1200 h

(http://hycom.org/dataserver)

surface current velocity (v 

direction) approx. 9 km daily 1992-2012

Hybrid vertical coordinate system 

ocean model (HYCOM) global 

reanalysis (GLBu0.08 grid, 

experiments 19.0 and 19.1), 

1200 h 

(http://hycom.org/dataserver)

surface current divergence approx. 9 km daily 1992-2012 derived from surface current 

velocity

surface current vorticity approx. 9 km daily 1992-2012 derived from surface current 

velocity

sea surface height approx. 28 km daily 1993-2016

AVISO Global 'DT all sat' MADT 

daily sea surface altimetry 

product

(https://www.aviso.altimetry.fr/en/

data/products/sea-surfaceheight-products/global/madt-huv.html) 

sea surface height 

standard deviation approx. 28 km daily 1993-2016 derived from sea surface height

probability of anticyclonic 

eddy ring approx. 28 km daily 1993-2016 derived from sea surface height8 

probability of cyclonic eddy 

ring approx. 28 km daily 1993-2016 derived from sea surface height8 

sea surface temperature approx. 1 km daily 1985-2010

AVHRR Pathfinder (Peter 

Cornillon, University of Rhode 

Island; 

http://www.sstfronts.org/opendap

/)

sea surface temperature 

standard deviation approx. 1 km daily 1985-2010 derived from sea surface temperature

probability of sea surface 

temperature front9 approx. 1 km monthly 2002-2013

derived from monthly front 

presence/absence in turn derived 

from daily NASA MUR SST data

(Peter Miller, Remote Sensing 

Group, Plymouth Marine 

Laboratory, UK; Miller and 

Christodoulou 2014)

upwelling index approx. 28 km monthly 1999-2009

NOAA CoastWatch grids derived 

from NASA Quick Scatterometer 

(QuikSCAT) data

(http://coastwatch.pfeg.noaa.gov/

erddap/griddap/erdQSstressmda

y.html)

Atmospheric (seasonal 

climatologies)

wind stress (x direction) approx. 28 km monthly 1999-2009

NOAA CoastWatch grids derived 

from NASA Quick Scatterometer 

(QuikSCAT) data

(http://coastwatch.pfeg.noaa.gov/

erddap/griddap/erdQSstressmda

y.html)
41 

Predictor variable Native spatial 

resolution

Native temporal 

resolution Years Source

wind stress (y direction) approx. 28 km monthly 1999-2009

NOAA CoastWatch grids derived 

from NASA Quick Scatterometer 

(QuikSCAT) data

(http://coastwatch.pfeg.noaa.gov/

erddap/griddap/erdQSstressmda

y.html)

wind divergence approx. 28 km monthly 1999-2009

NOAA CoastWatch grids derived 

from NASA Quick Scatterometer 

(QuikSCAT) data

(http://coastwatch.pfeg.noaa.gov/

erddap/griddap/erdQSstressmda

y.html)

1 Current and one-year lag

2 Oblique Mercator projection: azimuth = 40, origin = 75°W 35°N, k = 0.9996, ellipse = GRS80, earth 

datum = NAD83

3 Derived using ArcGIS 10.2 Spatial Analyst Tools

4 Datasets were converted to a common 100-m grid then mosaicked in order of decreasing priority: 

Okeanos Explorer, ACUMEN, CCOM, GMRT, CRM, and GEBCO

5 Derived using ArcGIS 10.4 Spatial Analyst Tools

6 Calculated from mosaicked 100-m depth grids that were smoothed using a Gaussian low-pass filter for 

each spatial scale

7 Derived using the ArcGIS 10.4 extension DEM Surface Tools (Jenness 2013) 

8 Derived using Duke University’s Marine Geospatial Ecology Lab's Marine Geospatial Ecology Tools 

(MGET) 0.8a64 for ArcGIS 10.4.1 (Okubo-Weiss algorithm)

9 A low-pass filter was applied to the climatologies to remove banding artifacts
42 

Table 4a. Pairwise Spearman rank correlation coefficients for spatial predictor variables (spring – March-May). High correlations are highlighted in 

yellow (>0.7), orange (>0.8), and red (>0.9). projected latitude distance to land depth slope (2 km) slope (10 km) slope of slope (10 km)

planform curvature (10 km)profile curvature(10 km)chlorophyll-a turbiditysurface current velocity (u)surface current velocity (v)surface current divergencesurface current vorticitysea surface heightsea surface height SDanticyclonic eddy probabilitycyclonic eddy probabilitysea surface temperaturesea surface temperature SDSST front probabilityupwelling indexwind stress (x)

wind stress (y)

wind divergence

projected longitude 0.84 0.30 -0.26 0.08 0.04 -0.02 0.02 0.12 0.49 0.02 -0.16 -0.36 0.12 0.00 -0.38 0.10 0.16 -0.27 -0.73 0.44 0.03 0.26 0.53 -0.50 0.00

projected latitude -0.18 0.21 0.01 -0.03 -0.03 0.02 0.05 0.82 0.32 -0.22 -0.38 0.18 0.02 -0.65 -0.20 -0.17 -0.56 -0.90 0.66 0.10 0.35 0.30 -0.65 0.12

distance to land -0.90 0.18 0.17 0.05 -0.03 0.15 -0.53 -0.41 0.17 -0.03 -0.14 -0.08 0.41 0.67 0.63 0.52 0.24 -0.24 0.01 -0.14 0.61 0.20 -0.19

depth -0.27 -0.25 -0.07 0.07 -0.29 0.59 0.46 -0.21 -0.02 0.12 0.06 -0.46 -0.70 -0.67 -0.54 -0.29 0.29 0.00 0.17 -0.63 -0.21 0.26

slope (2 km) 0.90 0.64 0.01 0.19 -0.21 -0.21 0.17 0.16 -0.08 0.19 -0.08 0.22 0.29 0.10 0.17 0.05 0.44 0.00 0.24 -0.17 -0.05

slope (10 km) 0.75 0.02 0.17 -0.22 -0.19 0.17 0.17 -0.08 0.21 -0.09 0.19 0.28 0.10 0.21 0.08 0.54 -0.01 0.23 -0.19 -0.01

slope of slope

(10 km) 0.01 0.10 -0.16 -0.17 -0.01 0.06 -0.10 0.22 -0.13 -0.02 0.13 0.04 0.14 -0.02 0.35 0.00 0.00 -0.11 0.05

planform curvature 

(10 km) -0.22 0.03 0.03 0.01 0.00 0.03 -0.01 -0.03 -0.02 -0.03 -0.02 -0.01 0.00 0.01 0.05 0.01 -0.06 0.01

profile curvature 

(10 km) -0.07 -0.10 0.11 0.05 0.07 -0.06 0.00 0.16 0.23 0.02 0.01 0.05 -0.01 0.06 0.15 -0.06 -0.06

chlorophyll-a 0.65 -0.27 -0.36 0.22 0.13 -0.73 -0.42 -0.43 -0.60 -0.86 0.77 0.08 0.30 -0.02 -0.64 0.27

turbidity -0.25 -0.28 0.01 0.22 -0.47 -0.24 -0.28 -0.20 -0.49 0.59 0.16 0.01 -0.11 -0.37 0.32

surface current 

velocity (u) 0.62 0.04 0.02 0.18 0.49 0.14 0.09 0.45 -0.01 0.30 0.16 0.34 -0.10 -0.26

surface current 

velocity (v) 0.07 0.04 0.22 0.27 0.01 0.01 0.58 -0.24 0.21 0.11 0.05 0.04 -0.29

surface current 

divergence -0.01 -0.21 -0.05 -0.06 -0.23 -0.14 0.17 -0.06 0.29 0.01 -0.27 -0.02

surface current 

vorticity -0.39 0.02 -0.06 0.12 -0.08 0.28 0.35 -0.05 -0.08 -0.37 0.21

sea surface height 0.41 0.24 0.54 0.67 -0.76 -0.33 -0.36 0.13 0.82 -0.43

sea surface height 

SD 0.56 0.55 0.37 -0.06 0.16 -0.04 0.65 0.08 -0.37

anticyclonic eddy 

probability 0.38 0.22 -0.09 0.09 -0.05 0.38 0.05 -0.10

cyclonic eddy 

probability 0.49 -0.36 -0.02 -0.38 0.15 0.45 -0.11

sea surface 

temperature -0.66 0.07 -0.23 -0.07 0.55 -0.29

sea surface 

temperature SD 0.47 0.32 0.18 -0.78 0.34

SST front 

probability 0.13 0.29 -0.47 0.16

upwelling index 0.13 -0.43 0.05

wind stress (x) -0.20 -0.23

wind stress (y) -0.29
43 

Table 4b. Pairwise Spearman rank correlation coefficients for spatial predictor variables (summer – June-August). High correlations are 

highlighted in yellow (>0.7), orange (>0.8), and red (>0.9). projected latitude distance to land depth slope (2 km) slope (10 km) slope of slope (10 km) planform curvature (10 km) profile curvature (10 km)

chlorophyll-a turbiditysurface current velocity (u)surface current velocity (v)surface current divergencesurface current vorticitysea surface heightsea surface height SDanticyclonic eddy probabilitycyclonic eddy probabilitysea surface temperaturesea surface temperature SDSST front probabilityupwelling indexwind stress (x)

wind stress (y)

wind divergence

projected longitude 0.84 0.30 -0.26 0.08 0.04 -0.02 0.02 0.12 0.33 0.31 -0.15 -0.30 0.08 0.00 -0.34 0.11 0.16 -0.22 -0.84 0.81 0.50 0.38 0.55 -0.22 -0.13

projected latitude -0.18 0.21 0.01 -0.03 -0.03 0.02 0.05 0.72 0.65 -0.22 -0.34 0.22 0.01 -0.65 -0.25 -0.15 -0.52 -0.94 0.90 0.70 0.51 0.33 -0.51 -0.03

distance to land -0.90 0.18 0.17 0.05 -0.03 0.15 -0.70 -0.58 0.17 -0.02 -0.27 -0.07 0.51 0.75 0.57 0.57 0.11 -0.07 -0.28 -0.16 0.56 0.58 -0.23

depth -0.27 -0.25 -0.07 0.07 -0.29 0.74 0.65 -0.18 -0.05 0.26 0.05 -0.54 -0.77 -0.65 -0.57 -0.16 0.13 0.32 0.18 -0.58 -0.60 0.28

slope (2 km) 0.90 0.64 0.01 0.19 -0.20 -0.23 0.21 0.20 -0.15 0.16 -0.01 0.26 0.32 0.10 0.07 -0.03 0.17 0.09 0.22 -0.05 -0.14

slope (10 km) 0.75 0.02 0.17 -0.20 -0.21 0.23 0.20 -0.18 0.18 -0.02 0.24 0.31 0.10 0.11 -0.05 0.21 0.07 0.21 -0.08 -0.09

slope of slope

(10 km) 0.01 0.10 -0.13 -0.13 0.03 0.10 -0.17 0.19 -0.11 0.03 0.12 0.01 0.07 -0.04 0.13 0.05 -0.01 -0.15 0.00

planform curvature 

(10 km) -0.22 0.05 0.05 0.01 0.00 -0.01 0.04 -0.03 -0.02 -0.03 0.00 -0.02 0.02 0.03 -0.01 -0.01 -0.07 0.02

profile curvature 

(10 km) -0.10 -0.11 0.12 0.10 0.04 -0.09 0.03 0.18 0.25 0.02 -0.03 0.03 -0.01 0.11 0.16 0.06 -0.14

chlorophyll-a 0.91 -0.22 -0.21 0.36 0.14 -0.79 -0.58 -0.51 -0.66 -0.66 0.64 0.67 0.40 -0.14 -0.72 0.21

turbidity -0.23 -0.28 0.31 0.16 -0.70 -0.49 -0.44 -0.54 -0.63 0.63 0.64 0.29 -0.08 -0.63 0.20

surface current 

velocity (u) 0.66 0.00 0.02 0.22 0.43 0.15 0.06 0.42 -0.39 0.09 0.10 0.27 0.16 -0.24

surface current 

velocity (v) 0.10 0.05 0.21 0.24 0.09 -0.02 0.51 -0.50 0.01 0.08 -0.01 0.05 -0.21

surface current 

divergence 0.00 -0.29 -0.18 -0.18 -0.28 -0.15 0.17 0.19 0.27 0.01 -0.27 0.04

surface current 

vorticity -0.32 0.06 -0.11 0.14 -0.02 0.07 0.22 -0.10 -0.08 -0.32 0.15

sea surface height 0.54 0.27 0.56 0.64 -0.67 -0.61 -0.44 0.20 0.85 -0.25

sea surface height 

SD 0.53 0.63 0.30 -0.30 -0.15 -0.13 0.63 0.59 -0.33

anticyclonic eddy 

probability 0.29 0.12 -0.08 -0.20 -0.03 0.36 0.28 -0.17

cyclonic eddy 

probability 0.43 -0.42 -0.45 -0.43 0.20 0.59 -0.06

sea surface 

temperature -0.96 -0.55 -0.39 -0.23 0.49 -0.07

sea surface 

temperature SD 0.60 0.41 0.24 -0.54 0.11

SST front 

probability 0.45 0.24 -0.58 -0.06

upwelling index 0.16 -0.37 -0.09

wind stress (x) 0.39 -0.25

wind stress (y) -0.27
44 

Table 4c. Pairwise Spearman rank correlation coefficients for spatial predictor variables (fall – September-November). High correlations are 

highlighted in yellow (>0.7), orange (>0.8), and red (>0.9). projected latitude distance to land depth slope (2 km) slope (10 km) slope of slope (10 km) planform curvature (10 km) profile curvature (10 km)

chlorophyll-a turbiditysurface current velocity (u)surface current velocity (v)surface current divergencesurface current vorticitysea surface heightsea surface height SDanticyclonic eddy probabilitycyclonic eddy probabilitysea surface temperaturesea surface temperature SDSST front probabilityupwelling indexwind stress (x)

wind stress (y)

wind divergence

projected longitude 0.84 0.30 -0.26 0.08 0.04 -0.02 0.02 0.12 0.29 -0.25 0.02 -0.19 0.10 0.00 -0.35 0.11 0.14 -0.20 -0.75 0.47 0.29 0.25 0.91 0.23 0.16

projected latitude -0.18 0.21 0.01 -0.03 -0.03 0.02 0.05 0.70 0.16 0.00 -0.30 0.18 0.02 -0.65 -0.23 -0.15 -0.52 -0.93 0.67 0.52 0.29 0.94 -0.07 0.26

distance to land -0.90 0.18 0.17 0.05 -0.03 0.15 -0.71 -0.67 0.10 0.14 -0.17 -0.09 0.50 0.71 0.54 0.57 0.24 -0.25 -0.35 -0.04 0.05 0.42 -0.15

depth -0.27 -0.25 -0.07 0.07 -0.29 0.77 0.71 -0.16 -0.22 0.14 0.07 -0.53 -0.73 -0.57 -0.59 -0.29 0.31 0.40 0.05 -0.01 -0.45 0.18

slope (2 km) 0.90 0.64 0.01 0.19 -0.20 -0.33 0.20 0.27 -0.14 0.19 -0.07 0.21 0.25 0.12 0.15 -0.13 0.19 0.05 0.07 -0.13 -0.08

slope (10 km) 0.75 0.02 0.17 -0.20 -0.32 0.20 0.27 -0.15 0.21 -0.08 0.19 0.23 0.12 0.18 -0.12 0.26 0.03 0.04 -0.20 -0.07

slope of slope

(10 km) 0.01 0.10 -0.12 -0.19 0.01 0.15 -0.13 0.23 -0.15 -0.02 0.08 0.05 0.13 -0.10 0.20 0.01 -0.01 -0.13 0.00

planform curvature 

(10 km) -0.22 0.04 0.06 0.03 0.00 0.01 0.00 -0.03 -0.01 0.00 -0.01 -0.01 -0.01 0.04 0.02 0.02 -0.05 0.01

profile curvature 

(10 km) -0.11 -0.24 0.13 0.12 0.05 -0.08 0.02 0.16 0.19 0.03 0.00 0.00 -0.06 0.09 0.09 0.04 -0.03

chlorophyll-a 0.65 -0.03 -0.32 0.23 0.14 -0.77 -0.54 -0.47 -0.68 -0.74 0.68 0.62 0.15 0.53 -0.44 0.25

turbidity -0.11 -0.24 0.10 0.19 -0.41 -0.45 -0.49 -0.33 -0.29 0.32 0.32 -0.13 -0.04 -0.43 0.09

surface current 

velocity (u) 0.59 0.07 0.03 0.07 0.40 0.08 -0.03 0.23 -0.20 0.26 0.11 0.02 -0.29 -0.26

surface current 

velocity (v) 0.01 0.03 0.25 0.34 0.14 0.08 0.53 -0.54 0.03 0.09 -0.26 -0.20 -0.34

surface current 

divergence 0.01 -0.17 -0.06 -0.04 -0.18 -0.14 0.15 0.15 0.10 0.11 -0.16 0.01

surface current 

vorticity -0.39 0.04 -0.02 0.16 -0.04 0.20 0.32 -0.31 -0.02 -0.32 -0.01

sea surface height 0.48 0.24 0.55 0.67 -0.72 -0.67 -0.11 -0.53 0.49 -0.28

sea surface height 

SD 0.52 0.64 0.35 -0.31 -0.20 -0.13 -0.05 0.19 -0.32

anticyclonic eddy 

probability 0.34 0.19 -0.12 -0.20 -0.02 0.00 0.15 -0.07

cyclonic eddy 

probability 0.48 -0.39 -0.42 -0.37 -0.38 0.27 -0.17

sea surface 

temperature -0.81 -0.40 -0.18 -0.85 0.01 -0.37

sea surface 

temperature SD 0.48 0.06 0.60 -0.26 0.39

SST front 

probability 0.14 0.45 -0.59 0.05

upwelling index 0.28 -0.04 0.04

wind stress (x) 0.01 0.21

wind stress (y) -0.04
45 

Table 4d. Pairwise Spearman rank correlation coefficients for spatial predictor variables (winter – December-February). High correlations are 

highlighted in yellow (>0.7), orange (>0.8), and red (>0.9). projected latitude distance to land depth slope (2 km) slope (10 km) slope of slope (10 km) planform curvature (10 km) profile curvature (10 km)

chlorophyll-a turbiditysurface current velocity (u)surface current velocity (v)surface current divergencesurface current vorticitysea surface heightsea surface height SDanticyclonic eddy probabilitycyclonic eddy probabilitysea surface temperaturesea surface temperature SDSST front probabilityupwelling indexwind stress (x)

wind stress (y)

wind divergence

projected longitude 0.84 0.30 -0.26 0.08 0.04 -0.02 0.02 0.12 0.16 -0.51 -0.05 -0.39 0.11 -0.02 -0.31 0.12 0.11 -0.24 -0.69 0.53 -0.20 0.27 0.84 -0.53 -0.03

projected latitude -0.18 0.21 0.01 -0.03 -0.03 0.02 0.05 0.61 -0.18 -0.10 -0.47 0.15 0.00 -0.61 -0.18 -0.19 -0.52 -0.89 0.73 0.07 0.41 0.71 -0.76 0.12

distance to land -0.90 0.18 0.17 0.05 -0.03 0.15 -0.78 -0.62 0.16 0.08 -0.12 -0.08 0.47 0.65 0.58 0.51 0.30 -0.23 -0.43 -0.22 0.42 0.32 -0.22

depth -0.27 -0.25 -0.07 0.07 -0.29 0.81 0.63 -0.22 -0.17 0.12 0.05 -0.53 -0.68 -0.61 -0.53 -0.35 0.29 0.43 0.26 -0.37 -0.32 0.29

slope (2 km) 0.90 0.64 0.01 0.19 -0.22 -0.40 0.16 0.16 -0.09 0.25 -0.07 0.18 0.22 0.11 0.20 0.07 0.22 0.07 0.21 -0.17 -0.10

slope (10 km) 0.75 0.02 0.17 -0.21 -0.39 0.16 0.16 -0.08 0.26 -0.09 0.16 0.21 0.11 0.23 0.09 0.32 0.07 0.20 -0.18 -0.06

slope of slope

(10 km) 0.01 0.10 -0.12 -0.22 -0.03 0.03 -0.06 0.26 -0.14 -0.06 0.06 0.00 0.16 0.00 0.22 0.06 0.04 -0.12 -0.04

planform curvature 

(10 km) -0.22 0.04 0.02 0.02 -0.01 0.00 -0.01 -0.03 -0.01 -0.03 -0.01 -0.01 0.02 0.03 0.03 0.03 -0.05 0.00

profile curvature 

(10 km) -0.13 -0.26 0.13 0.09 0.06 -0.04 0.01 0.15 0.21 0.02 0.02 0.05 -0.10 0.07 0.14 -0.04 -0.05

chlorophyll-a 0.44 -0.19 -0.32 0.19 0.15 -0.76 -0.52 -0.52 -0.59 -0.71 0.67 0.44 0.35 0.05 -0.62 0.38

turbidity -0.26 -0.10 0.00 -0.01 -0.06 -0.41 -0.44 -0.10 -0.06 -0.12 0.13 -0.17 -0.63 0.19 0.19

surface current 

velocity (u) 0.67 0.09 0.02 0.14 0.48 0.13 0.02 0.34 0.03 0.15 0.10 0.17 -0.08 -0.16

surface current 

velocity (v) 0.09 0.03 0.31 0.34 0.11 0.07 0.66 -0.29 0.13 0.04 -0.23 0.19 -0.26

surface current 

divergence -0.03 -0.15 0.02 -0.01 -0.19 -0.12 0.20 0.02 0.19 0.04 -0.16 0.06

surface current 

vorticity -0.36 0.06 -0.01 0.14 -0.04 0.22 0.34 -0.10 -0.03 -0.25 0.21

sea surface height 0.42 0.23 0.53 0.65 -0.73 -0.52 -0.35 -0.23 0.80 -0.43

sea surface height 

SD 0.49 0.53 0.33 -0.05 -0.17 -0.15 0.37 0.14 -0.21

anticyclonic eddy 

probability 0.30 0.25 -0.09 -0.18 -0.12 0.17 0.16 -0.08

cyclonic eddy 

probability 0.46 -0.39 -0.20 -0.46 -0.13 0.53 -0.05

sea surface 

temperature -0.71 -0.04 -0.28 -0.48 0.64 -0.31

sea surface 

temperature SD 0.48 0.40 0.54 -0.77 0.39

SST front 

probability 0.21 -0.05 -0.44 0.33

upwelling index 0.27 -0.45 0.05

wind stress (x) -0.58 -0.06

wind stress (y) -0.20
46 

Table 5. Base-learners employed in the boosted generalized additive modeling framework. Base-learner 

names are from the ‘mboost’ package for R (Hothorn et al. 2015; R Core Team. 2017). 

Name Description Predictor variables Model 

component

bols linear intercept p, μ, θ

bols linear survey platform p, μ, θ

brandom random effect survey ID θ

brandom random effect transect ID p, μ

bbs penalized regression spline1 year p, μ

bbs penalized regression spline1 day of year p, μ

btree tree2 all climate indices p, μ

bspatial penalized tensor product1 projected longitude

projected latitude p, μ

brad penalized radial basis3 projected longitude

projected latitude p, μ

btree tree4 distance to land and all bathymetric, 

oceanographic, and atmospheric variables p, μ

1 P-spline basis

2 Maximum depth = 1

3 Matern correlation function

4 Maximum depth = 5
47 

Table 6. Best models with model performance metrics.

Species Season Distribution Start 

values1 

Number of 

boosting 

iterations

Percent 

deviance 

explained2

AUC3 

Spearman 

rank 

correlation

Gaussian 

rank 

correlation4 

Median 

absolute 

error5

Mean 

absolute 

error5

Median 

bias5 Mean bias5 Root mean 

square error

ARTE summer ZINB GAM 19,537 0.49 0.95 0.10 0.14 0.08 1.66 0.08 -0.11 0.36

ATPU spring ZINB GAM 18,622 0.41 0.92 0.15 0.19 0.14 1.69 0.14 -0.09 0.27

ATPU summer ZIP simple 19,991 0.49 0.98 0.14 0.18 0.08 1.58 0.08 0.01 0.24

ATPU fall ZIP GAM 19,998 0.46 0.96 0.09 0.13 0.15 1.76 0.15 0.07 0.08

ATPU winter ZINB GAM 19,928 0.41 0.95 0.15 0.20 0.12 1.59 0.11 -0.05 0.22

AUSH spring ZINB simple 17,476 0.71 0.98 0.11 0.17 0.09 1.06 0.09 -0.31 0.33

AUSH summer ZINB GAM 14,710 0.52 0.95 0.22 0.28 0.06 1.40 0.06 -0.11 1.39

AUSH fall ZINB GAM 17,610 0.64 0.98 0.14 0.20 0.16 1.36 0.16 -0.21 0.82

AUSH winter ZIP GAM 7,332 0.78 1.00 0.12 0.20 0.00 0.88 0.00 -0.10 0.11

BCPE spring ZIP GAM 14,891 0.74 1.00 0.12 0.19 0.00 0.75 0.00 -0.01 0.12

BCPE summer ZINB GAM 17,498 0.67 0.99 0.16 0.22 0.00 1.35 0.00 -0.19 0.59

BCPE fall ZIP GAM 6,766 0.83 1.00 0.08 0.14 0.00 0.49 0.00 -0.01 0.07

BCPE winter ZIP GAM 9,734 0.86 0.99 0.09 0.16 0.00 0.39 0.00 -0.01 0.05

BLGU summer ZIP simple 19,954 0.59 0.99 0.08 0.13 0.15 1.50 0.15 0.09 0.09

BLKI spring ZINB GAM 19,122 0.51 0.91 0.22 0.25 0.05 1.43 0.05 -0.32 2.82

BLKI fall ZINB simple 18,020 0.56 0.95 0.35 0.40 0.02 1.17 0.02 -0.28 14.91

BLKI winter ZIP GAM 17,973 0.58 0.93 0.48 0.51 0.10 1.07 0.06 -0.07 9.05

BLSC spring ZINB GAM 19,910 0.59 0.95 0.19 0.24 0.02 1.24 0.02 -0.53 24.71

BLSC fall ZINB GAM 19,858 0.48 0.96 0.16 0.21 0.13 1.63 0.13 -0.12 6.34

BLSC winter ZIP simple 9,897 0.58 0.91 0.25 0.26 0.16 1.67 0.14 -0.10 71.51

BOGU spring ZINB GAM 18,906 0.39 0.92 0.17 0.18 0.12 1.63 0.11 -0.27 8.79

BOGU fall ZINB simple 19,892 0.46 0.90 0.11 0.14 0.19 1.95 0.19 0.14 1.49

BOGU winter ZIP GAM 6,819 0.52 0.87 0.25 0.29 0.23 1.40 0.20 -0.15 5.08

BRPE spring ZIP simple 17,205 0.45 0.99 0.08 0.13 0.20 1.51 0.20 -0.01 0.45

BRPE summer ZIP GAM 18,954 0.28 0.99 0.09 0.14 0.00 1.11 0.00 -0.40 0.48

BRPE fall ZIP GAM 16,244 0.43 0.99 0.10 0.16 0.02 1.25 0.02 -0.24 0.67
48 

Species Season Distribution Start 

values1 

Number of 

boosting 

iterations

Percent 

deviance 

explained2

AUC3 

Spearman 

rank 

correlation

Gaussian 

rank 

correlation4 

Median 

absolute 

error5

Mean 

absolute 

error5

Median 

bias5 Mean bias5 Root mean 

square error

BRPE winter ZIP GAM 17,807 0.59 0.98 0.08 0.12 0.00 1.43 0.00 -0.29 0.51

BRSP summer ZIP GAM 19,249 0.53 0.96 0.13 0.19 0.08 1.42 0.08 -0.07 0.18

BRTE summer ZINB GAM 19,970 0.60 0.95 0.08 0.11 0.10 1.66 0.10 -0.21 0.19

BRTE fall ZIP simple 19,987 0.60 1.00 0.07 0.12 0.30 1.70 0.30 0.33 0.10

COEI spring ZIP simple 4,941 0.43 0.98 0.25 0.31 0.20 1.82 0.19 0.17 260.83

COEI summer ZINB simple 19,997 0.62 0.99 0.10 0.15 0.00 2.62 0.00 1.22 29.83

COEI fall ZIP GAM 8,979 0.64 0.95 0.16 0.21 0.00 1.00 0.00 -1.00 109.23

COEI winter ZINB simple 19,988 0.37 0.94 0.37 0.40 0.03 1.56 0.02 -0.20 425.86

COLO spring ZIP GAM 19,981 0.40 0.91 0.41 0.44 0.19 1.24 0.13 -0.07 1.09

COLO summer ZINB GAM 19,988 0.45 0.94 0.11 0.14 0.08 1.77 0.08 -0.08 0.10

COLO fall ZIP GAM 19,980 0.44 0.95 0.28 0.34 0.11 1.31 0.09 -0.05 0.53

COLO winter ZINB GAM 19,518 0.32 0.82 0.35 0.37 0.50 1.39 0.33 -0.16 3.80

COMU spring ZIP GAM 18,023 0.38 0.92 0.12 0.15 0.22 1.64 0.22 -0.22 0.38

COMU winter ZIP GAM 19,994 0.49 0.95 0.14 0.19 0.02 1.42 0.02 -0.09 0.37

COSH spring ZINB GAM 19,522 0.56 0.98 0.11 0.15 0.07 1.65 0.07 -0.13 0.18

COSH summer ZINB GAM 18,254 0.34 0.84 0.31 0.33 0.32 1.40 0.23 -0.29 5.42

COSH fall ZINB GAM 19,994 0.43 0.89 0.28 0.32 0.15 1.38 0.12 -0.21 3.01

COTE spring ZINB GAM 18,997 0.52 0.97 0.23 0.29 0.02 1.51 0.02 0.03 1.31

COTE summer ZIP GAM 16,732 0.35 0.94 0.30 0.34 0.14 1.47 0.13 -0.02 2.51

COTE fall ZIP simple 13,236 0.53 0.94 0.20 0.24 0.17 1.46 0.15 -0.02 3.60

DCCO spring ZINB GAM 19,171 0.46 0.95 0.11 0.13 0.06 1.48 0.06 -0.43 7.43

DCCO summer ZIP simple 13,180 0.53 0.91 0.09 0.11 0.24 1.79 0.24 -0.10 2.45

DCCO fall ZINB simple 19,999 0.37 0.89 0.10 0.13 0.23 1.54 0.22 -0.42 4.77

DCCO winter ZINB GAM 8,179 0.65 0.91 0.10 0.13 0.09 1.22 0.09 -0.73 3.57

DOVE spring ZIP GAM 16,947 0.53 0.94 0.21 0.28 0.08 1.21 0.07 -0.18 1.39

DOVE summer ZINB GAM 20,000 0.67 0.98 0.06 0.10 0.01 1.36 0.01 -0.30 0.10

DOVE fall ZINB simple 19,792 0.62 0.98 0.18 0.25 0.02 1.07 0.02 -0.20 1.85
49 

Species Season Distribution Start 

values1 

Number of 

boosting 

iterations

Percent 

deviance 

explained2

AUC3 

Spearman 

rank 

correlation

Gaussian 

rank 

correlation4 

Median 

absolute 

error5

Mean 

absolute 

error5

Median 

bias5 Mean bias5 Root mean 

square error

DOVE winter ZINB GAM 18,923 0.43 0.91 0.27 0.33 0.08 1.17 0.07 -0.34 4.18

GBBG spring ZINB GAM 19,978 0.48 0.86 0.41 0.43 0.13 1.24 0.10 -0.33 18.20

GBBG summer ZINB GAM 19,074 0.47 0.90 0.38 0.40 0.09 1.35 0.05 -0.12 4.36

GBBG fall ZINB GAM 16,719 0.40 0.88 0.46 0.47 0.15 1.30 0.09 -0.13 12.13

GBBG winter ZIP GAM 16,954 0.66 0.89 0.42 0.45 0.10 1.12 0.07 -0.04 11.69

GRSH spring ZINB GAM 19,995 0.72 0.98 0.24 0.32 0.00 1.00 0.00 -0.22 3.36

GRSH summer ZIP GAM 11,197 0.78 0.90 0.50 0.50 0.08 1.03 0.06 -0.04 129.86

GRSH fall ZINB simple 19,961 0.50 0.95 0.59 0.60 0.03 1.08 0.01 -0.13 22.02

GRSH winter ZINB GAM 19,432 0.59 0.98 0.11 0.16 0.02 1.21 0.02 -0.33 0.71

GRSK fall ZIP simple 19,994 0.33 0.95 0.11 0.13 0.20 1.92 0.20 0.01 0.08

HERG spring ZINB simple 17,970 0.42 0.86 0.48 0.51 0.20 1.20 0.13 -0.21 22.70

HERG summer ZINB GAM 17,602 0.50 0.90 0.37 0.40 0.07 1.35 0.04 -0.13 3.92

HERG fall ZINB GAM 17,790 0.42 0.87 0.51 0.52 0.20 1.21 0.10 -0.12 14.84

HERG winter ZIP simple 14,194 0.46 0.85 0.43 0.45 0.30 1.25 0.18 -0.04 10.69

HOGR winter ZIP simple 19,993 0.41 0.94 0.09 0.12 0.26 1.71 0.26 -0.08 0.15

LAGU spring ZINB GAM 17,678 0.52 0.95 0.23 0.29 0.08 1.43 0.07 -0.15 0.75

LAGU summer ZINB GAM 19,997 0.58 0.95 0.30 0.36 0.06 1.24 0.05 -0.13 1.38

LAGU fall ZIP GAM 15,380 0.55 0.94 0.31 0.36 0.06 1.34 0.05 0.00 3.41

LAGU winter ZINB GAM 19,978 0.58 0.97 0.10 0.16 0.09 1.42 0.09 -0.17 0.24

LESP spring ZIP simple 18,370 0.52 0.96 0.15 0.19 0.14 1.47 0.13 -0.06 0.81

LESP summer ZINB GAM 19,925 0.50 0.94 0.36 0.40 0.05 1.21 0.04 -0.18 2.92

LESP fall ZINB GAM 19,989 0.60 0.97 0.19 0.26 0.05 1.29 0.05 -0.23 0.59

LETE summer ZIP GAM 13,589 0.36 0.94 0.08 0.11 0.18 1.76 0.18 -0.02 0.48

LETE fall ZINB GAM 19,977 0.62 0.97 0.08 0.11 0.00 1.99 0.00 0.12 1.22

LTDU spring ZINB GAM 18,504 0.57 0.98 0.33 0.40 0.00 1.11 0.00 -0.51 224.59

LTDU fall ZINB GAM 19,850 0.66 0.99 0.19 0.27 0.03 1.27 0.03 -0.14 15.84

LTDU winter ZINB GAM 14,859 0.61 0.97 0.49 0.53 0.02 1.18 0.01 -0.24 98.27
50 

Species Season Distribution Start 

values1 

Number of 

boosting 

iterations

Percent 

deviance 

explained2

AUC3 

Spearman 

rank 

correlation

Gaussian 

rank 

correlation4 

Median 

absolute 

error5

Mean 

absolute 

error5

Median 

bias5 Mean bias5 Root mean 

square error

MASH spring ZIP GAM 19,089 0.36 0.93 0.09 0.12 0.32 1.67 0.32 -0.10 0.11

MASH summer ZIP simple 16,998 0.56 0.84 0.10 0.11 0.39 1.68 0.38 -0.26 1.97

MASH fall ZINB simple 19,995 0.29 0.89 0.12 0.15 0.26 1.73 0.25 -0.18 0.27

NOFU spring ZINB GAM 19,994 0.66 0.94 0.43 0.48 0.01 1.08 0.01 -0.26 24.71

NOFU summer ZINB simple 19,619 0.68 0.98 0.24 0.32 0.01 0.94 0.00 -0.45 21.93

NOFU fall ZINB GAM 19,836 0.62 0.96 0.35 0.42 0.00 1.18 0.00 0.13 9.88

NOFU winter ZINB GAM 19,175 0.67 0.98 0.39 0.46 0.01 0.98 0.00 -0.28 15.23

NOGA spring ZINB GAM 19,547 0.44 0.85 0.50 0.52 0.24 1.13 0.12 -0.24 14.42

NOGA summer ZIP GAM 19,934 0.45 0.94 0.27 0.32 0.07 1.39 0.06 -0.07 0.58

NOGA fall ZINB GAM 19,962 0.51 0.91 0.48 0.51 0.12 1.09 0.05 -0.13 4.10

NOGA winter ZINB GAM 18,330 0.41 0.81 0.46 0.47 0.29 1.19 0.16 -0.33 20.58

PAJA spring ZIP GAM 18,341 0.14 0.86 0.05 0.06 0.51 2.01 0.50 0.02 0.05

PAJA summer ZINB GAM 18,311 0.23 0.87 0.05 0.06 0.56 1.97 0.56 -0.01 0.06

PAJA fall ZINB GAM 19,997 0.25 0.83 0.08 0.08 0.51 1.97 0.50 -0.01 0.11

POJA spring ZINB GAM 18,347 0.53 0.96 0.09 0.13 0.27 1.68 0.27 -0.15 0.08

POJA summer ZINB simple 19,991 0.23 0.84 0.07 0.08 0.53 1.96 0.53 -0.02 0.08

POJA fall ZINB GAM 19,998 0.37 0.90 0.19 0.21 0.13 1.77 0.11 -0.06 0.22

RAZO spring ZIP GAM 17,136 0.43 0.94 0.27 0.32 0.08 1.40 0.07 -0.04 1.72

RAZO summer ZIP simple 19,812 0.45 0.98 0.07 0.11 0.16 1.78 0.16 0.14 0.15

RAZO fall ZINB GAM 19,999 0.58 0.97 0.11 0.15 0.09 1.91 0.09 0.29 1.05

RAZO winter ZIP GAM 17,825 0.46 0.91 0.33 0.36 0.12 1.39 0.10 -0.06 4.74

RBGU spring ZINB GAM 19,449 0.46 0.94 0.12 0.16 0.21 1.70 0.21 -0.08 0.36

RBGU summer ZINB GAM 19,580 0.54 0.96 0.06 0.08 0.20 1.44 0.20 -0.54 0.21

RBGU fall ZINB GAM 18,878 0.41 0.90 0.14 0.17 0.15 1.49 0.15 -0.34 1.20

RBGU winter ZIP simple 6,421 0.41 0.88 0.19 0.23 0.30 1.71 0.29 -0.02 2.24

RBME spring ZINB GAM 19,998 0.51 0.93 0.07 0.10 0.11 1.59 0.11 -0.33 0.43

RBME winter ZINB GAM 19,649 0.31 0.89 0.07 0.07 0.53 1.68 0.53 -0.31 1.48
51 

Species Season Distribution Start 

values1 

Number of 

boosting 

iterations

Percent 

deviance 

explained2

AUC3 

Spearman 

rank 

correlation

Gaussian 

rank 

correlation4 

Median 

absolute 

error5

Mean 

absolute 

error5

Median 

bias5 Mean bias5 Root mean 

square error

REPH spring ZIP GAM 15,470 0.74 0.96 0.19 0.23 0.01 1.23 0.01 -0.27 90.82

REPH summer ZINB simple 18,496 0.31 0.94 0.11 0.16 0.22 1.29 0.22 -0.63 96.77

REPH fall ZINB GAM 19,241 0.38 0.95 0.14 0.17 0.19 1.69 0.19 -0.16 1.45

RNPH spring ZIP simple 17,426 0.31 0.93 0.10 0.13 0.21 1.86 0.20 -0.04 3.18

RNPH summer ZINB GAM 19,999 0.36 0.94 0.10 0.13 0.15 1.56 0.15 -0.27 1.98

RNPH fall ZINB GAM 19,998 0.32 0.89 0.09 0.11 0.48 1.93 0.48 -0.02 0.87

ROST spring ZINB GAM 19,999 0.57 0.98 0.07 0.10 0.12 1.59 0.12 -0.32 0.25

ROST summer ZINB GAM 19,139 0.59 0.97 0.11 0.16 0.00 1.65 0.00 0.00 0.69

ROST fall ZIP simple 12,530 0.46 0.96 0.07 0.10 0.36 2.38 0.36 0.89 0.53

ROYT spring ZINB GAM 19,317 0.68 0.98 0.15 0.22 0.00 1.21 0.00 -0.28 0.58

ROYT summer ZINB GAM 19,995 0.65 0.98 0.14 0.20 0.01 1.32 0.01 -0.15 0.23

ROYT fall ZINB simple 19,985 0.53 0.97 0.15 0.20 0.05 1.61 0.05 -0.08 0.36

RTLO spring ZINB GAM 19,564 0.48 0.93 0.34 0.39 0.08 1.25 0.06 -0.18 1.87

RTLO fall ZIP GAM 17,216 0.45 0.96 0.16 0.21 0.05 1.47 0.05 -0.07 0.61

RTLO winter ZIP GAM 17,073 0.42 0.87 0.33 0.37 0.27 1.42 0.19 -0.03 2.40

SOSH spring ZINB GAM 18,420 0.53 0.95 0.26 0.32 0.02 1.22 0.02 -0.32 5.06

SOSH summer ZINB simple 19,999 0.71 0.93 0.30 0.34 0.01 1.17 0.01 -0.52 57.21

SOSH fall ZINB simple 17,594 0.39 0.92 0.08 0.09 0.27 1.85 0.27 -0.12 0.50

SOTE spring ZIP GAM 11,468 0.63 1.00 0.07 0.13 0.00 1.14 0.00 -0.06 0.84

SOTE summer ZIP GAM 18,403 0.56 0.98 0.09 0.12 0.03 1.67 0.03 -0.07 0.50

SPSK summer ZIP simple 19,994 0.35 0.92 0.07 0.10 0.32 1.85 0.32 -0.02 0.06

SPSK fall ZIP GAM 19,990 0.38 0.96 0.09 0.13 0.16 1.78 0.16 -0.07 0.09

SUSC spring ZINB GAM 17,711 0.57 0.98 0.26 0.33 0.02 1.62 0.02 0.27 9.04

SUSC fall ZINB GAM 19,811 0.55 0.97 0.23 0.31 0.01 1.04 0.01 -0.17 10.73

SUSC winter ZINB simple 19,787 0.61 0.97 0.37 0.44 0.05 1.07 0.04 -0.18 16.31

TBMU spring ZINB GAM 19,993 0.44 0.95 0.15 0.20 0.06 1.43 0.06 -0.20 0.90

TBMU winter ZIP GAM 19,736 0.50 0.97 0.11 0.15 0.16 1.52 0.16 0.03 0.24
52 

Species Season Distribution Start 

values1 

Number of 

boosting 

iterations

Percent 

deviance 

explained2

AUC3 

Spearman 

rank 

correlation

Gaussian 

rank 

correlation4 

Median 

absolute 

error5

Mean 

absolute 

error5

Median 

bias5 Mean bias5 Root mean 

square error

WISP spring ZINB GAM 18,134 0.58 0.97 0.36 0.41 0.02 1.33 0.01 -0.10 10.65

WISP summer ZINB GAM 19,033 0.38 0.86 0.50 0.52 0.18 1.18 0.10 -0.31 28.20

WISP fall ZIP simple 18,090 0.25 0.96 0.29 0.35 0.08 1.46 0.07 -0.01 3.76

WWSC spring ZIP GAM 9,376 0.49 0.94 0.19 0.24 0.20 1.34 0.20 -0.28 27.23

WWSC fall ZIP GAM 11,397 0.56 0.97 0.20 0.27 0.02 1.23 0.02 -0.09 5.90

WWSC winter ZINB GAM 19,906 0.53 0.95 0.32 0.37 0.05 1.25 0.05 -0.34 20.66

1 See Section 2.5.5

2 To calculate percent deviance explained, the saturated likelihood was assumed to be the maximum possible likelihood value, and the null 

likelihood was calculated from an intercepts-only zero-inflated model fit to the data (unpublished)

3 Area under the receiver operating characteristic curve for classifying observed counts as 0 or >0

4 Boudt et al. (2012) and Bodenhofer et al. (2013)

5 Predicted values minus observed values expressed as a proportion of the mean observed value (i.e., mean count)
53 

Figure 1. Study area with BOEM Wind Energy Lease and Planning Areas overlaid (approximate 

boundaries current as of 2017-04-10).
54 

Figure 2. Area surveyed by year. 
55 

Figure 3. Total area surveyed each season within the study area, binned into 10 x 10 km cells. 

White areas represent no survey effort. See Appendix A for maps of individual survey effort.
56

Step 1. Data preparation

segmented survey data

+

predictor variables

intersect data combined data

Tuning

Boosting iteration

Step 2. Model fitting

boosting parameters

• learning rate

• number of iterations

identify optimal values that

minimize prediction error

(using cross-validation)

model

fitting

fitted

model

Step 4. Model performance

Step 3. Prediction across space

Step 6. Precision of predictions

Step 5. Model selection

predictor variables fitted functional

relationships

sea surface temperature SD

log(mean count)

+

prediction

predicted relative density

predicted relative density segmented survey data

+

statistical fit performance metrics

multiple models select model with

best performance

best model

best model re-sample data & re-fit 

model 200 times

multiple predictions coefficient of variation

Negative log-likelihood (out-of-bag data)

Figure 4. Schematic overview of statistical modeling process. See Section 2 Methods for details.
57

Figure 5. Example maps of predicted relative density and its coefficient of variation (CV) for one 

species (Atlantic Puffin) and season (spring).
58

Figure 6a. Relative importance of predictor variables for p (red) and μ (blue) components of best spring (March-May) models. Areas of 

circles are proportional to importance.
59 

Figure 6b. Relative importance of predictor variables for p (red) and μ (blue) components of best summer (June-August) models. 

Areas of circles are proportional to importance.
60 

Figure 6c. Relative importance of predictor variables for p (red) and μ (blue) components of best fall (September-November) models. 

Areas of circles are proportional to importance.
61 

Figure 6d. Relative importance of predictor variables for p (red) and μ (blue) components of best winter (December-February) models. 

Areas of circles are proportional to importance.
62 

Figure 7. Example annual maps of predicted relative density and its coefficient of variation (CV) 

for one species (Atlantic Puffin).
63 

Figure 8. Example map of predicted relative density for Coastal Waterfowl species group during 

winter. Coastal Waterfowl included Black Scoter, Common Eider, Common Loon, Long-tailed 

Duck, Red-throated Loon, Surf Scoter, and White-winged Scoter.
64 

Figure 9. Maps of predicted relative density for phalarope multi-species models for spring and 

fall. These models were fit to pooled sighting data for Red Phalarope, Red-necked Phalarope, 

and unidentified phalaropes.
65 

Glossary

Note: many of these definitions are specific to the context of this report. 

Area under the receiver operating characteristic (ROC) curve (AUC) – An ROC curve is a graphical 

representation of how well a model can discriminate between (or predict) two categories of data (e.g., 

presence/absence), and the AUC is the integral of this curve. AUC values range from 0-1 where a value of 

0.5 indicates model performance equivalent to random; a value >0.5 indicates performance better than 

random; and a value <0.5 indicates performance worse than random. Thus, higher AUC values indicate 

better model performance.

Autocorrelation – Correlation between data points or residual errors that are close in space and/or time. 

Spatial data often exhibit autocorrelation, and not accounting for it in predictive models can bias model 

predictions and artificially inflate statistical precision and significance.

Base-learner – In a boosted generalized additive modeling framework, a relatively simple model relating 

the response variable to a predictor variable(s). One base-learner is selected in each boosting iteration, 

and the final model is essentially the sum of modeled relationships across the selected base-learners.

Boosting – Iterative model fitting technique. In each iteration a single base-learner is selected. Each 

selected base learner’s contribution to the final model is controlled by the learning rate.

Bootstrap (non-parametric) – A data re-sampling technique for estimating the statistical uncertainty in 

model predictions. A dataset of size n is re-sampled with replacement x times to derive x new datasets of 

size n. The model is fit to each new dataset to derive x predictions. The variability across these x

predictions can then be used to evaluate their precision (e.g., coefficient of variation).

Climatology – Long-term spatial pattern in an environmental variable. For example, average values 

across years at different locations in space during a given annual time period (e.g., monthly, seasonal).

Coefficient of variation (CV) – Measure of dispersion for a distribution, representing the standard 

deviation (SD) as a proportion of the mean. In this report, CV is calculated from the mean and SD of 

distributions of bootstrapped model predictions, so the CV actually reflects the standard error (SE) of the 

prediction relative to the mean prediction. A larger CV indicates more variation (uncertainty) in the 

prediction relative to the mean prediction and thus lower precision and higher uncertainty. 

Cross-validation – A technique for evaluating the predictive ability of a fitted model. The data are divided 

into in-bag data and out-of-bag data, the model is fit to the in-bag data, and then the fitted model’s ability 

to predict the out-of-bag data is measured.

Ensemble model – A model created by combining multiple models into a single model. In the context of 

boosting, models are fit iteratively (i.e., after a model is fit, the remaining variation is used to fit the next 

model) and then combined.

Generalized additive model – A model whose response variable is the sum of multiple, potentially nonlinear (e.g., smooth) functional relationships with predictor variables. A link function is employed for the 

response variable.
66 

In-bag – In a cross-validation context, the subset of data, drawn at random without replacement from the 

full dataset, that the model is fit to in order to estimate the model parameters.

Learning rate – In a boosting context, the degree to which each base learner contributes to the final model. 

The optimal learning rate is one that results in well-defined model convergence sooner than later.

Offset (boosting) – In a boosting context, the initial values of each model component that are used to 

initialize the boosting algorithm.

Offset (effort) – In count models, a model term that accounts for survey effort by enforcing a proportional 

relationship between the expected count and effort (e.g., area surveyed).

Out-of-bag – In a cross-validation context, the subset of data, drawn at random without replacement from 

the full dataset, from which the predictive performance of the fitted model is assessed.

Percent deviance explained (PDE) – Measure of the percentage of variation in the data explained by a 

model beyond that explained by the simplest model without predictor variables. Values normally range 

from 0-100%, although negative values are possible. Higher values indicate better model performance. 

PDE is a generalized model analogue of the coefficient of determination (R2

).

Predictive performance – The ability of a model to explain variation in data that the model was not fit to 

(e.g., out-of-bag data in a cross-validation framework). 

Predictor – An independent variable in a model that is used to explain variation in the response.

Rank correlation coefficient – Measure of the correspondence between observed and predicted values. 

Values range between -1 and 1 where a value of 0 indicates no correspondence between observed and 

predicted values; a value >0 indicates positive correspondence between observed and predicted values; 

and a value <0 indicates negative correspondence between observed and predicted values. Thus, positive 

values closer to 1 indicate better model performance.

Relative density – Model predicted values that are proportional to the expected number of birds per unit 

area, but that do not represent the actual expected number of birds per unit area. An index of density. 

Relative variable importance – Measure of the importance of a predictor variable in terms of the 

frequency with which that predictor occurred in the selected base-learners across boosting iterations and 

that predictor’s ability to explain variation in the data when it was selected. Relative variable importance 

was scaled so that it summed to 1 across predictors.

Re-sampling – A method of using randomly drawn subsets of data to estimate statistical precision (e.g., 

bootstrapping model predictions) or to perform model validation (e.g., cross-validation).

Residual error – Difference between observed data and corresponding model predictions. In this report, 

residual error is expressed as a percentage of the mean of the data.

Response – The dependent variable in a model representing the quantity of interest for which predictions 

are to be made.
67 

Spatial predictive modeling – Modeling technique whereby relationships between environmental 

predictors and a response variable are estimated for areas with survey data, and then these relationships 

are used to predict the response as a function of the same environmental predictors in areas without 

survey data.

Stochastic gradient boosting – A type of boosting whereby the data are sub-sampled in each iteration 

before the base-learner(s) is fit to the gradient.

Tree depth – In a tree base-learner, a parameter that controls the number of allowable nodes in the tree, 

thereby controlling the number of possible interactions between predictor variables. Greater tree depth 

allows for more interactions.

Tuning – Procedure by which model fitting parameters are adjusted to maximize the predictive

performance of a model. For example, in boosting the learning rate and number of boosting iterations are 

adjusted during tuning. 

Zero-inflated negative binomial (ZINB) distribution – A statistical distribution used to model count data 

that accounts for a large number of zeroes and an overdispersed count distribution (e.g., because of 

aggregative behavior of animals).

Zero-inflated Poisson (ZIP) distribution – A statistical distribution used to model count data that accounts 

for a large number of zeroes.

The Department of the Interior Mission

As the Nation's principal conservation agency, the Department of the Interior has 

responsibility for most of our nationally owned public lands and natural resources. This 

includes fostering sound use of our land and water resources; protecting our fish, 

wildlife, and biological diversity; preserving the environmental and cultural values of 

our national parks and historical places; and providing for the enjoyment of life through 

outdoor recreation. The Department assesses our energy and mineral resources and 

works to ensure that their development is in the best interests of all our people by 

encouraging stewardship and citizen participation in their care. The Department also 

has a major responsibility for American Indian reservation communities and for people 

who live in island territories under US administration.

The Bureau of Ocean Energy Management

As a bureau of the Department of the Interior, the Bureau of Ocean Energy (BOEM) 

primary responsibilities are to manage the mineral resources located on the Nation's 

Outer Continental Shelf (OCS) in an environmentally sound and safe manner.

The BOEM Environmental Studies Program

The mission of the Environmental Studies Program (ESP) is to provide the information 

needed to predict, assess, and manage impacts from offshore energy and marine 

mineral exploration, development, and production activities on human, marine, and 

coastal environments.