Analysis of the North American Breeding Bird Survey Using Hierarchical Models

Abstract

We analyzed population change for 420 bird species from the North American Breeding Bird Survey (BBS) using a hierarchical log-linear model and compared the results with those obtained through route-regression analysis. Survey-wide trend estimates based on the hierarchical model were generally more precise than estimates from the earlier analysis. No consistent pattern of differences existed in the magnitude of trends between the analysis methods. Survey-wide trend estimates changed substantially for 15 species between route-regression and hierarchical-model analyses. We compared regional estimates for states, provinces, and Bird Conservation Regions; differences observed in these regional analyses are likely a consequence of the route-regression procedure's inadequate accommodation of temporal differences in survey effort. We used species-specific hierarchical-model results to estimate composite change for groups of birds associated with major habitats and migration types. Grassland, aridland, and eastern-forest-obligate bird species declined, whereas urban—suburban species increased over the interval 1968–2008. No migration status group experienced significant changes, although Nearctic—Neotropical migrant species showed intervals of decline and permanent resident species increased almost 20% during the interval. Hierarchical-model results better portrayed patterns of population change over time than route-regression results. We recommend use of hierarchical models for BBS analyses.

Resumen

Analizamos los cambios poblacionales de 420 especies con base en los censos de aves reproductivas de Norte América (BBS, por sus siglas en inglés) mediante un modelo jerárquico log-lineal y comparamos los resultados con los obtenidos mediante análisis de regresión por rutas. Los estimados de tendencias a nivel de todos los censos basados en modelos jerárquicos fueron más precisos que los estimados del análisis previo. No existió un patrón consistente de diferencias en la magnitud de las tendencias entre los métodos de análisis. Los estimados de tendencias a nivel de todos los censos cambiaron sustancialmente entre los análisis de regresión y de modelos jerárquicos para 15 especies. Comparamos los estimados regionales para estados, provincias y regiones de conservación de aves; las diferencias observadas en esos análisis regionales probablemente son consecuencia de la consideración inadecuada de las diferencias temporales en el esfuerzo de muestreo que hace el procedimiento de regresión. Utilizamos los resultados de modelos jerárquicos específicos de cada especie para estimar el cambio conjunto de grupos de aves asociadas con los principales ambientes y tipos de migración. Las poblaciones de especies restringidas a pastizales, zonas áridas y a bosques del oriente disminuyeron, mientras que las de aves de ambientes urbanos-suburbanos aumentaron entre 1968 y 2008. Ninguna agrupación basada en el tipo de migración experimentó cambios significativos, aunque las poblaciones de especies migratorias neártico-neotropicales disminuyeron en algunos intervalos y las poblaciones de especies residentes permanentes aumentaron en casi un 20% durante el intervalo. Los resultados de los modelos jerárquicos ilustraron los patrones de cambio poblacional de mejor manera que los resultados de los análisis de regresión. Recomendamos el uso de modelos jerárquicos para analizar datos del BBS.

THE NORTH AMERICAN Breeding Bird Survey (BBS) provides information regarding population change for >420 bird species in the United States and Canada (Sauer et al. 2008). The BBS is the primary data source for modeling consequences of changes in land uses, climate, or other possible stressors on most North American bird populations (U.S. NABCI Committee 2009); for many species, the BBS provides the only information for modeling population dynamics. Its extensive use by conservation planners (e.g., Rich et al. 2004) and managers (e.g., Butcher et al. 2007) illustrates its importance to bird conservation, science, and management.

From 1994 to 2007, analysis of the BBS was conducted using a route-regression procedure (Geissler and Sauer 1990, Link and Sauer 1994). In route regression, a species' trend is defined, for each survey route, as the slope of a Poisson regression with covariates to accommodate observer differences (Link and Sauer 1994). Regional trend is estimated as a weighted average of route-specific estimates. After estimating regional trends, annual indices are estimated as residuals from a regionally smoothed population trajectory (Sauer and Geissler 1990). Route regression can be applied to any subset of BBS data, and even extremely unbalanced data sets with many missing count years can be summarized using route regression to provide an estimate of regional trend. The analysis can also be implemented with limited computer resources. These attributes made route regression an important tool for summary analyses of BBS data (e.g., Sauer et al. 2008). However, it has limitations as an analytical technique. Although conceptually complete, implementation of the method is based more on computational convenience than on a formal model. Regional trend estimates are weighted averages in which ad hoc weights are used to combine route data, and variances of the estimate must be estimated by bootstrapping. The two-step process of first estimating trend and then estimating annual indices also makes it difficult to estimate the precision of the residual indices, leading to estimates of the annual indices without confidence intervals. Finally, the method's focus on estimating trends can obscure changes in survey efficiency over time (Link and Sauer 2002).

Link and Sauer (2002) proposed a log-linear hierarchical model for analysis of population change in the BBS. This model has features similar to the Sauer and Geissler (1990) model, in that year effects vary around a long-term trend, but in the hierarchical analysis the model is fit regionally and year effects are included in the model as random effects. Annual indices and trend are estimated directly from year effects and other model components; interval estimates are calculated for annual indices of abundance by region. Although structural similarities exist between the hierarchical-model analysis and route regression, the hierarchical model has the advantages that all components are estimated in the analysis, avoiding ad hoc precision weightings, and that variance can be directly estimated for indices and trends. The hierarchical model can be conveniently implemented using Bayesian methods (Link and Sauer 2002).

We estimated population change for 420 bird species from BBS data using the hierarchical model. We compared results from the hierarchical models with those from a BBS analysis using the route-regression method. We focused on estimates of long-term trend from the analyses as overall summaries of the results; an appendix containing a more comprehensive analysis and summary is available online (Appendix S1; see below). We estimated trends by species at the scale of the entire survey, by states and provinces, and by Bird Conservation Regions (BCRs; Sauer et al. 2003) and tested for consistent differences in estimated trends and relative precision between the analysis methods. We also compared hierarchical-model results of two alternative definitions of trend at the survey-wide scale. To summarize patterns of population change for groups of species, we used a hierarchical model to estimate composite trends for species groups used in the “State of the Birds” report (U.S. NABCI Committee 2009).

A comprehensive summary of population trends and annual indices for species is posted online as a supplement to the present study (see Acknowledgments). Species information includes graphs of population change, along with a comparative summary of estimates of change derived from the hierarchical model and route-regression analysis for each of the 420 species summarized in online Appendix S2 (see Acknowledgments). Survey-wide and regional (state, province, and BCR) population change for 1966– 2008 and 1999–2008 was calculated using hierarchical models and is presented with 95% credible (or confidence) intervals. We also present the long-term change estimate based on route regression, along with sample sizes (number of standard Breeding Bird Survey [BBS] routes used in the route-regression analysis) and relative abundance (mean count of birds of that species on BBS routes). We indicate whether estimated trends differed between the long-term hierarchical-model and route-regression estimates. We also estimated the proportion of the population in each BCR, and we note whether the species is of management concern in the region (as determined by Partners in Flight species-prioritization efforts; Panjabi et al. 2005). Time-series graphs show annual indices based on hierarchical-model results.

Methods

North American Breeding Bird Survey

Breeding Bird Surveys (Sauer et al. 2008) are conducted by volunteer bird watchers. Each observer selects a morning to conduct a survey, primarily during June, although late May dates are allowed in southern states and early July dates are allowed in northern provinces. Surveys are point counts conducted at 50 stops placed along randomly preselected roadside routes. Counts start 30 min before local sunrise, and stops are 0.8 km apart. At each stop, the observer conducts a 3-min count of all birds heard or seen within 400 m (0.25 miles). Of the >5,100 survey routes that have been established, ∼2,500 are surveyed each year. The density of survey routes varies greatly across North America, with lower route densities in the western United States and very few routes in northern Canada (Fig. 1).

Fig. 1.

Map of North America, showing locations of North American Breeding Bird Survey (BBS) routes. State and provincial boundaries are shown, along with boundaries of Bird Conservation Regions. The heavy line indicates the northern edge of the area covered in our analysis.

Route Regression

Route regression has been the primary method of summary of BBS data for >20 years (e.g., Holmes and Sherry 1988, lgl and Johnson 1997), and the method is described in detail by Geissler and Sauer (1990) and Link and Sauer (1994). In the present version of BBS route-regression analysis, trend is estimated for each survey route using estimating equations to directly fit a multiplicative model in which counts are regressed on year with observer identity included as a categorical predictor (Link and Sauer 1994). Trend is defined as the slope of the regression on year. Although estimates of trends from the estimating equations are identical to those from a Poisson regression with log links, estimating equations do not require or use the distributional assumptions needed for the Poisson regression. Composite regional trends are estimated as weighted averages of route slopes for the time intervals of interest. Weighting by mean abundance of routes and by a variance weight (Sauer et al. 2003) limits the influence of poorly surveyed or low-abundance routes on the estimates (Geissler and Sauer 1990). For regional (multistratum) estimates, an additional area weight is included to accommodate differences in population sizes among strata. Bootstrapping is used to reduce bias and estimate variances of regional trend estimates. The bootstrap procedure consists of repeatedly subsampling collections of n routes from the n routes in a region. A trend estimate is computed for each subsample; the mean of these estimates is used as a regional trend estimate, and their variance is used to estimate its precision. Annual indices of abundance are estimated as yearly deviations from a long-term trend line (Sauer and Geissler 1990).

Hierarchical Model

Link and Sauer (2002) and Sauer and Link (2002) described a hierarchical model for analysis of BBS data. This model allows year, stratum, and observer effects to be described by parameters that are random variables and directly accommodates the multiple scales at which data are collected in the survey. Regional summaries are defined as functions of the underlying regional abundance and population change parameters. Bayesian methods, in which inference is based on distributions of parameters conditional on the data (i.e., posterior distributions), provide both a conceptual basis for interpretation of results and a computational approach for fitting the hierarchical model to BBS data. Bayesian analyses require specification of prior distributions of parameters and sampling distributions of the data conditioned on the parameters; from these distributions the posterior distribution of the parameters conditional on the data can be computed. The posterior distribution is used to make probability statements about the parameters. For instance, the posterior mean (or median) might be used as a point estimator, and percentiles of the posterior distribution can be used to create credible intervals (Bayesian confidence intervals). Unlike classical (frequentist) statistical inference, statements based on the posterior distributions of parameters are direct statements about probabilities associated with the parameters. Posterior distributions of parameters are often difficult to compute analytically but can be approximated using simulation-based Markov chain Monte Carlo (MCMC) methods (Lunn et al. 2000). Markov chain samples can be used to calculate means, medians, and credible intervals for the posterior distributions of the parameters of interest. The MCMC approach also has the convenient property that the posterior distribution of derived parameters (functions of the parameters sampled in the MCMC) can be directly calculated from the Markov chain samples (Link and Barker 2010).

In the hierarchical model for BBS analyses, population change is modeled with an overdispersed Poisson regression. Counts Y_i,j,t (i for stratum, j for unique combinations of route and observer, and t for year) are independent Poisson random variables with means λ_i,j,t that are log-linear functions of explanatory variables:  

$$$$

Explanatory variables in Equation 1 are stratum-specific intercepts (S) and slopes (β; t* is the baseline year), observer-route combinations (ω), year (γ), start-up (η, with I[j,t] an indicator that takes the value 1 for an observer's first year of survey on a route, 0 otherwise), and overdispersion effects (ε). Year effects are random; explicit modeling of this variance component is necessitated by variation in sample sizes (Link and Sauer 2002).

The hierarchical model requires distributional assumptions and specification of prior distributions. In the present analysis, parameters, S_i η, and β_i were given diffuse (essentially flat) normal distributions (mean 0, variance 10⁶), whereas observer-route effects (ω), year effects (γ), and overdispersion effects (ε) were treated as mean-zero normal random variables. The variances of the observer effects and overdispersion effects were set as constants and , respectively, while the variance of the year effects (γ) was allowed to vary among strata (). All variances were assigned diffuse inverse gamma prior distributions (scale and shape parameters set to 0.001).

Stratum-specific annual indices of abundance are exponentiated sums of year, stratum, and trend effects, scaled by the proportion of routes in the stratum on which the species was observed (z_i). Variance components are also added to accommodate asymmetries in estimating the mean from the log normal distribution and to scale indices to a level similar to mean counts on a route in a region:  

$$$$

n_i,t is an index to the number of birds per route in stratum i at year t. Link and Sauer (2002) defined annual indices without the variance components; we include them here to provide a scaling of the indices that is consistent with historical analyses. Indices for stratum totals are calculated as N_i,t = A_i n_i,t, where A_i is the area of the stratum. To obtain indices for larger areas (groups of strata; e.g., states, BCRs, countries), we sum the N_i,t over the relevant i. For presentation, the composite indices N_t are scaled by the total areas, obtaining a summary scaled to birds per route:  

$$$$

The β_i values from the model could be used as a measure of composite change within the model but are relevant only at the scale of the strata and for the full time interval for which the model is fit. To estimate composite population change at the scale of states, BCRs, and other regions of interest, and to permit change estimation for any time interval, composite summaries of population change are defined as functions of the annual indices n_it. Link and Sauer (1998, 2002) defined trend as an interval-specific geometric mean of yearly changes in population size, expressed as a percentage. Thus, the trend from year t_a to year t_b for stratum i is 100(B_i - 1)%, where  

$$$$

The composite trend is calculated as 100( - 1)%, using the composite indices  

$$$$

and

Link and Sauer's (2002) definition of trend is based on ratios of parameters associated with single years. Alternative definitions of trend abound. For instance, Thomas et al. (2004) defined trend as the ratio of endpoints from semiparametric models that describe a nonlinear curve through the data. One common definition of trend is the ratio of predicted values from a linear regression through the annual indices. To document the consequences of an alternative definition of trend, we also estimated trend at the survey-wide scale as the slope of a regression of log-transformed annual indices on year. These regression-based trend estimates were calculated as derived statistics during the MCMC, and we compared these estimates to the geometric mean estimator defined above.

We used the program WINBUGS to conduct the MCMC analysis for strata within states and provinces (Lunn et al. 2000). We discarded at least 20,000 initial samples from each Markov chain as burn-in, and then ran another 20,000 iterations to obtain results for estimating the posterior distributions. Some species required additional iterations before usable estimates were obtained, and for others, the large data sets proved difficult to manage and we could only summarize 10,000 replicates. To minimize storage needs we used results from every second iteration (i.e., “thinned” the Markov chains by a factor of 2) to calculate estimates and credible intervals. We also output the MCMC replicates for additional summaries at larger geographic scales.

The hierarchical model requires sufficient samples that span the period of interest to allow estimation of the year effects and annual indices. Regions with very few survey routes or with missing years of data in the interval of interest produced very imprecise results, and occasionally inclusion of these results led to extremely imprecise regional estimates. In those instances, we removed the region that produced the imprecision and reran the regional analysis. We have noted the eliminated regions in Appendix S2.

Scales of Summary

BBS trend and annual indices were estimated at the scale of states and provinces, physiographic regions, and for the entire surveyed area (excluding Alaska; Fig. 1). Conservation of migratory birds is increasingly focused on modeling and estimation of population change for BCRs, which are ecoregions developed to provide a common framework for bird conservation. Following a comparative analysis of BBS data, we now use BCRs as strata (Sauer et al. 2003). Because the BBS was originally stratified and coordinated within states and provinces, we retained the states and provinces as components of the stratification, estimated population change for BCRs within states or provinces as our fundamental strata, and aggregated these regions to estimate composite trends within states and provinces, BCRs, and larger regions.

Summary of Start-up Effects

Kendall et al. (1996) documented the presence of start-up effects for observers in the BBS. In their analysis, removing each observer's first year count reduced the trend estimates in 275 of 416 species. Because the effect size was small, removal of initial counts for each observer as an ad hoc approach to accommodate start-up effects was not implemented in the route-regression procedure. The hierarchical model includes a parameter to model start-up effect (η) and controls for it in the analysis. We estimated start-up effects for species, and we present estimates of the parameters and their credible intervals.

Comparison of Hierarchical-model with Route-regression Results

We conducted a route-regression analysis of trends and annual indices by state or province and management units for the period 1966–2008 and compared the results with hierarchical-model results. The BBS was implemented first in the eastern United States in 1966, in the central United States in 1967, and in the western United States and southern Canada by 1968. Trend estimates in these analyses reflect only years for which there were data (i.e., for species that occur only in the western United States, the trend estimates are based on the interval 1968–2008).

Evaluation of precision of results.—Route-regression analysis and other statistical summaries provide confidence intervals, but the hierarchical-model analysis provides credible intervals. Conceptually, these intervals are not equivalent, because they represent the outcomes of frequentist and Bayesian inference, respectively. In frequentist statistics, the confidence interval is a random range of values in which, over repeated sampling under replicate conditions, the unknown but fixed parameter will occur with a certain probability. In Bayesian statistics, probability distributions are used to quantify uncertainty about parameters and the credible interval is calculated from the posterior distribution of the parameter, providing the probability that the parameter occurs within an interval. Regardless of their derivation, we used the width of the interval to measure the precision of the results. We considered an estimate imprecise if the half-width of the credible interval was larger than 3% year^-1 for the species (Sauer et al. 2003). We summarized the number of imprecise estimates for species by region. We also evaluated whether consistent differences in precision occurred between methods. We calculated the mean difference in width of the 95% intervals between the estimates (hierarchical model - route regression) and estimated the 95% confidence intervals of these differences by region. We considered a difference significant if the confidence intervals did not overlap zero. In both of these analyses, we restricted the comparisons to species found on 14 or more survey routes in the region.

Comparison of estimated trend.—For each species within each region, we compared the hierarchical-model and route-regression trend estimates. Trend estimates were considered in disagreement if the 95% credible interval from the hierarchical-model analysis did not overlap the 95% confidence intervals from the route-regression analysis. We evaluated whether consistent differences occurred between methods by calculating the mean difference between the estimates (hierarchical model - route regression) and the 95% confidence intervals of the mean difference. We considered a difference to be significant if the confidence intervals did not contain zero.

Multispecies Summaries

Aside from documenting individual species patterns, BBS results are often used to evaluate population change for groups of species (Sauer and Link 2002). The “State of the Birds” report (U.S. NABCI Committee 2009) used data from the BBS, the Christmas Bird Count, and waterfowl surveys conducted by the U.S. Fish and Wildlife Service to describe change in population status for several groups of species in the United States. Collections of species were defined as those that occur in major biomes (Arctic; aridland; eastern, central, western, boreal, and subtropical forests; and grassland). Species that were limited to one of these biomes were considered obligates. A variety of secondary habitat groups were also defined, including birds that use urban and suburban habitats, wetland birds, and marsh-specialist birds. A generalist category described widespread bird species that occur in three or more major biomes. Birds were also grouped by migration status: permanent residents, temperate migrants, and Nearctic—Neotropical migrants. Finally, exotic (non-native) species were considered a group. See the report and the associated website (U.S. NABCI Committee 2009) for additional information about the species groups.

For each of these species groups, a composite summary of population change was constructed using BBS data, following the approach of Gregory et al. (2005). Survey-wide annual indices for each species in the group were scaled to set the first year as 1.0, and then indices from subsequent years for each species were averaged to form a composite time series. We used 1968 as the first (or base) year because no data were available from the western United States prior to 1968. Simple averages of the species indices do not provide a valid summary of population change, because estimates for species vary greatly in quality of information. Any simple summary will be dominated by extreme, yet imprecise, estimates (e.g., Sauer and Link 2002). Sauer and Link (2002) developed a hierarchical model that summarized a collection of trend estimates (e.g., population change from the base year to any future year). Population trend parameters (percent per year) for each species were assumed to be hierarchically governed by a common mean and variance. These hyperparameters were used to estimate the trend parameters for the collection of estimated trend. This model has the benefit that the mean is based on the parameters, not on the often imprecise estimates. The model was fit in a Bayesian analysis using MCMC methods (Sauer and Link 2002).

Gregory et al. (2005) estimated composite change using a geometric mean by year among the species indices. To be consistent with their approach (see rationale in Gregory et al. 2005), we used a model in which the logarithms of the change parameters () for s = 1,…, n species are assumed to be normally distributed (i.e., [ln()|β_s,]= N[ln(β_s),]) and distributions for the trend parameters (β_s) are lognormal ([ln(β_s)|µ,τ²] = N[µ,τ²]). Note that year indexing is suppressed here; the βs values are estimated and summarized for each year. Our Bayesian analysis of this model used diffuse priors for µ (normal distributions with means 0, variances 10⁶) and τ² (inverse gamma distributions with scale and shape parameters set to 0.001). We applied the model to data for each year for each collection of species, and estimated posterior distributions of β_s and µ. We present the median and percentile credible intervals of the posterior distribution of exp(µ) as composite change from the base year to the year of interest. We provide the summary results for the final year, showing the total percentage change from 1968 to 2008 for the groups, and graphs of change over time for selected groups.

Results

Survey-wide trends.—Trend estimates from the hierarchical-model and route-regression analyses for 420 species are presented in Appendix S2. Mean difference in trends (hierarchical-model estimate — route-regression estimate) among the species was very small (-0.1% year^-1; confidence intervals overlapped zero). The 95% credible intervals of the hierarchical-model estimates were 0.9% year^-1 shorter than 95% confidence intervals from the route regression, a significant result (Table 1). Seventy-three species had long (>6%) credible intervals.

Fifteen species (3.6%) had credible intervals in the hierarchicalmodel analysis that did not overlap the confidence intervals of the route-regression analysis. Of these, 10 species trends agreed in direction between methods but differed in magnitude. For example, Northern Bobwhite (Colinus virginianus) had an estimated trend of -3.7% year^-1 from the hierarchical-model analysis and -2.9% year^-1 from the route-regression analysis. We note also that Henslow's Sparrow (Ammodramus henslowii) had trend estimates that, while agreeing as to sign, were quite different between the analyses. Five species, Common Merganser (Mergus merganser), Mourning Dove (Zenaida macroura), Northern Waterthrush (Parkesia noveboracensis), Carolina Chickadee (Poecile carolinensis), and Blue-gray Gnatcatcher (Polioptila caerulea), had trends that disagreed in sign between the methods. We review the patterns for these species in more detail below.

Differences in results associated with trend definitions appeared to be minor. Regression-based estimates of trend from the hierarchical models produced results that were very similar to those from the geometric mean trend estimator, with a median difference of 0.025% year^-1 (range: -1.09 to 4.26, percentile confidence interval: -0.46 to 1.03) between the two estimates (Fig. 2). The slightly higher geometric mean trend estimate for some of the more dramatically increasing species (Fig. 2) is likely a reflection of consistently positive residual values for the final year in the interval. We note that none of those differences fall outside the credible intervals of the alternative estimates. Median absolute differences (disregarding sign) were 0.12% year^-1 (range: 0.001–4.26). The greatest difference in trend estimates was for Tricolored Blackbird, which appeared as an outlier from the general correspondence of trends (Fig. 2), but estimates for the species were very imprecise. We use the geometric mean trend estimator for regional comparisons with route-regression results. Detailed summaries of individual species results, including annual indices of abundance at the survey-wide scale, are presented in Appendix S1.

Bird Conservation Region trends.—Among BCRs, the percentage of species with inconsistent estimates (nonoverlapping confidence intervals) between methods averaged 2.2%, varying from zero in six BCRs to 6.1% (9 of 147 species) in the Appalachian Mountains (Table 1). On average, 16.5% of species (range: 4.7% [Appalachian Mountains] to 37.5% [Boreal Taiga Plains]) in a BCR had credible intervals with half-widths >3% year^-1; thus, an estimated trend of <3% year^-1 would not be significant. The average difference in trends (hierarchical model — route regression) among BCRs was 0.2%, with a range of -1.5 (Northern Rockies) to 1.7 (Gulf Coastal Prairie), and 40.5% (13 of 32) of these differences had confidence intervals that did not overlap zero. Among those BCRs with non-overlapping confidence intervals, the hierarchical models tended to have more positive trends in seven, the route regression in six. The credible intervals of trends estimated from hierarchical models averaged 3.6% year^-1 (range: 8.1–1.6) smaller than the confidence intervals of trends from route regressions.

Table 1.

Hierarchical-model and route-regression estimates of trend for Bird Conservation Regions. For each region, the number of species (n), the number of species for which the route-regression trends were greater (>) or less (<) than the hierarchical-model trends (as shown by nonoverlapping confidence and credible intervals), the number of species for which a 3% year^-1 change would not have been detectable (V), the mean difference in estimated trends (hierarchical model — route regression) and the confidence interval associated with the difference, and the difference in 95% intervals (hierarchical model credible interval width — route regression confidence interval width) and the confidence interval associated with the difference.

State and provincial trends.—Among states and provinces, the percentage of species with inconsistent estimates (nonoverlapping confidence intervals) between methods averaged 1.1%, varying from zero in 24 states or provinces to 4.7% in Maryland. On average, 15.3% of species (range: 0% [Delaware] to 48.8% [Nevada]) in a state or province had credible intervals large enough that a 3% year^-1 trend would not be significant. Average difference in trends (hierarchical model — route regression) among states and provinces was 0.04% (range: -2.0 [Washington] to 1.7 [Iowa]). As with BCRs, a large percentage of these differences in trends (23.6%, 13 of 55) had confidence intervals that did not overlap zero; of these, nine had hierarchical-model estimates that tended to be larger and three had hierarchical-model estimates that tended to be smaller. The credible intervals of the hierarchical models averaged 3.8% year^-1 (range: 8.5–1.4) smaller than the confidence intervals of the route regression.

Start-up effects.—Of the 420 species, the mean η was -0.039 (CI: -0.053 to -0.026), indicating that counts in the first year of the survey by an observer for an average species were ∼4% smaller than what would have been expected, having controlled for population change. First-year effect varied among species (Appendix S2): significant positive first-year effects were documented in 21 species, whereas significant negative first-year effects occurred in 127 species. The positive first-year effects did not appear to be associated with large differences in trends between the methods, with a mean difference between the route-regression and hierarchicalmodel trend estimates for positive effects of 0.03% year^-1 (-1.10 to 1.16); significant negative effects were associated with larger differences, of 0.2% year^-1 (-0.02 to 0.48).

Fig. 2.

Scatterplot by species of population change from 1966 to 2008, estimated from North American Breeding Bird Survey data using a logregression based estimator (Regression) and the geometric mean estimator based on the ratio of endpoints. A 1:1 line is added to show equality.

Group summaries.—Regional habitat-obligate bird species, migration-status groups, and the secondary groups showed widely varying patterns of population change. Composite changes from 1968 to 2008 showed that grassland, aridland, and eastern-forestobligate bird species collectively declined, while generalist and urban— suburban species collectively increased. None of the migration status categories showed significant changes, although Nearctic—Neotropical migrant species showed indications of decline and permanent resident species increased almost 20% over the interval (Table 2).

Table 2.

Total percentage change from 1968 to 2008 for 13 species groups, with credible intervals and number of species in the group (n).

Fig. 3.

Plots of composite population status from 1968 to 2008 for collections of species grouped by migration habits: (A) Nearctic—Neotropical migrant, (B) temperate migrant, and (C) permanent resident species. Hierarchical-model estimates of total percentage change from 1968 (Index) are shown as lines with year markers, and credible intervals are shown as lines.

Plots of composite population status over time for the species indicate that Nearctic—Neotropical and temperate migrants experienced recent increases after declines in the 1970s, whereas permanent resident species had stable populations in early years and increased in recent years (Fig. 3). Trajectories of habitat-obligate groupings also showed variation in changes over time, with grassland, eastern forest, and boreal forest birds declining in the early years of the survey and then varying in pattern in recent years, aridland birds declining over the entire survey period with few years of increase, general stability of western forest birds, and increases in populations of urban—suburban birds after 1980 (Fig. 4). Detailed information regarding species in these groups and additional summaries of species-group results are presented online on the BBS and State of the Birds websites (see Acknowledgments).

Fig. 4.

Plots of composite population status from 1968 to 2008 for breeding-habitat-obligate species groups: (A) grassland, (B) eastern forest, (C) boreal forest, (D) aridland, and (E) western forest birds. Composite population status is also presented for (F) urban—suburban bird species. Hierarchicalmodel estimates of total percentage change from 1968 (Index) are shown as lines with year markers, and credible intervals are shown as lines.

Discussion

The BBS samples the continental United States and southern Canada. It has been partially implemented in Alaska and Puerto Rico, and plans are underway for a Mexican BBS. Data quality and consistency vary greatly over this large geographic range, and it is a challenge to develop statistical models that accommodate this variation in the quality of the survey while controlling for factors that influence detectability of birds along BBS routes. Route regression provided a means to control for observer differences and to accommodate variation in the quality of information from survey routes but did not provide a comprehensive framework for analysis of BBS data that could fully accommodate the multiscale aspects of the survey, be sensitive to the temporal patterns of data quality, and permit a full evaluation of covariates that influence both detection of birds and actual population change. Hierarchical models provide this framework and provide clear advantages over the route-regression analysis in terms of efficiency of estimation. Hierarchical models and MCMC fitting methods are easily modified to include covariates and spatial modeling (e.g., Thogmartin et al. 2004, Nielson et al. 2008), and they permit composite analysis of the BBS and other surveys, such as the Christmas Bird Count (Link and Sauer 2007). We view implementation of the stratum-based hierarchical model as an important step in providing better analyses of BBS data. We focused here on estimation of long-term trend as a primary component of population change; a more comprehensive comparative analysis of population change and annual indices is available online in the supplementary materials (see Acknowledgments).

Alternative Definitions of Trend and Their Consequences for Analysis of BBS Data

Our comparative analysis of two alternative definitions of trend indicated that choice of metric of trend did not have major consequences for the analysis of population change. This result is intuitive, because the extreme year indices used in the geometric mean definition of trends are also influential values in determining regression-based slopes. We note that the geometric mean change is a derived parameter based on a model including a loglinear regression on time. Estimation in this model-based context reduces its sensitivity to poorly informed endpoints.

Definitions of trend that are based on the notion of a long-term, consistent change that is an underlying characteristic of a time series (e.g., Dagum and Dagum 1988) have great conceptual appeal, but we suggest that the geometric mean definition is preferable in two situations that are frequently encountered in the BBS. First, many conservation applications such as state-of-the-birds summaries (U.S. NABCI Committee 2009) require estimates of interval-specific change for many time intervals, and the geometric mean estimate is clearly the appropriate metric for those analyses. Second, populations fluctuate, and Thomas et al. (2004) suggested that trend estimates based on the regression through the indices are likely to be inadequate as a model for population change in longer time series because they will not be sensitive to possibly important population changes occurring at the end of the time series. Examples of species undergoing large population fluctuations are easy to find in BBS data because populations fluctuate in response to severe winters or disease outbreaks (e.g., LaDeau et al. 2007), and species with rapidly increasing populations in recent years tended to have larger differences in the alternative trend estimates. We note that only Eastern Bluebirds had alternative trend estimates for which the credible intervals did not overlap. Examination of annual indices suggests that Eastern Bluebirds, after increasing in population size from the late 1970s, have experienced population declines in the most recent 10 years, possibly related to severe winters in the eastern United States. Regardless of how trend is defined, estimation of trend is readily accomplished using the hierarchical models we applied to the BBS.

Comparative Analysis of Route-regression and Hierarchical-model Trends

Precision.—When applied to BBS data, hierarchical models provided long-term trend estimates that generally had narrower credible intervals than the corresponding confidence intervals from route-regression analysis. This is true among BCRs and states and provinces, and the differences in precision are larger at these lower geographic scales than at the survey-wide scale. The lower precision of route-regression-based trend estimates was likely associated with the inefficiency of estimating route-specific trends with incomplete time series. Although weighting was used to limit the influence of routes with limited data on the overall estimate, most regions did not contain consistent information from all periods (Sauer et al. 2003); averaging route-specific trends based on different periods increased the variance of the route-regression trend estimates.

Differences in trend estimates.—Survey-wide, there were no consistent differences in mean estimates of population change between hierarchical-model and route-regression-based trend estimates. Averaged trend differences among BCR estimates and among state and province estimates were also not significantly different. However, many BCRs and states and provinces had small but significant differences in trend estimates between the two analyses. These differences were not consistently positive or negative, but they indicate that some component of the survey differentially influenced the results of one of the analyses.

The differences may have resulted from changes in numbers of routes surveyed over time. Route regression is based on estimated trends for individual routes, averaged to estimate an overall trend. Missing data on routes could have made the overall trend estimate unrepresentative of the entire interval. Although routes with more data (larger temporal spans with data, fewer observer changes) were weighted more heavily in the analysis, data from routes with short periods contributed to the overall trend estimate even though their data reflected a limited period. If the preponderance of routes was from a subinterval of the survey period, the change during that interval was likely influential in the overall trend estimate; if that period experienced drought or some other environmental stressor that influenced many species, it could have had a common influence on the species trend estimates. This confounding of number of surveys and trend was controlled in the hierarchical-model analysis. Estimates of trend from the hierarchical model were based on ratios of annual indices; the year-effects models controlled for differences in data quality among the years. In the year-effects model, limited data from portions of the interval will influence the precision of the year effects and, hence, the precision of annual indices and trends, but the trend estimate will always represent the appropriate interval.

The number of surveyed routes has changed dramatically over time. In 1970, 1,316 routes were surveyed; in 2008, 2,670 were surveyed, and an average of 58 more routes were surveyed each year from 1975 to 1993. This large increase raised the concern that recent data would tend to have a disproportionate effect on trend estimates, and in recent years, trend analyses were conducted that both included and omitted recently initiated routes to determine whether systematic differences were detectable in the analysis (J. R. Sauer unpubl. data). Although those analyses did not detect a significant bias associated with new routes in the context of the route-regression analysis, the differences that we documented between hierarchical-model trends and route-regression trends suggest the possibility that the changes in route coverage over time influenced the route-regression estimates.

We calculated the difference in average number of routes surveyed during 1969–1973 and 2004–2008 for each BCR, and correlated them with the mean difference in trends (hierarchical model - route regression) among species in the BCR. The correlation was -0.46 (n = 32), indicating that changes in number of routes were associated with differences in trends between methods. The average difference in routes was 71.4, 23.6, and 34.9 for BCRs with route-regression trends that were significantly larger, significantly smaller, or not significantly different than hierarchical-model trends, respectively. These results suggest that increased influence of the recent years' data led to more positive trends in the route-regression analysis. The species-group analyses indicate that in recent years most bird species had increasing populations; hence, a greater influence of these years in the route-regression analysis may have influenced the trend estimate for the entire interval.

Species-specific differences in trend estimates.—Within BCRs and states and provinces, only 2.3% of species had confidence intervals of trends from the route-regression analysis that did not overlap the credible intervals from the hierarchical-model analysis. Among the 83 species that differed, 62 cases are a consequence of differences in magnitudes of estimates (i.e., the direction of the change was the same in both analyses; only the magnitude differed). The other 21 (0.6% of the 3,652 estimates of species trends in the analysis) differed in sign as well as significance. For species with conflicting trend results at a survey-wide scale, it is informative to examine consistency between trends among BCRs to determine whether the overall difference is a consequence of a few very different regions or a consistent difference within regions. Common Mergansers were infrequently observed on BBS routes because they occur on wetlands that are often not visible from roadsides, and the species occurred in several regions that were poorly sampled. Although none of the BCR-scale results were significantly different between methods, results from the route-regression analysis were larger in magnitude and less precise than the hierarchical-model results in the influential Northern Rockies, Northern Pacific Rainforest, and Atlantic Northern Forest BCRs. Also, Coastal California, on the basis of a small sample of 11 routes, had large estimated trend of 19.5% year^-1 from the route regression but only 8.2% from the hierarchical-model analysis. The Mourning Dove is a widespread and generally abundant species in the continental United States, and the significant difference in results associated with the small difference of 0.2% year^-1 reflected the precision of the estimates. The difference was likely a consequence of substantial differences in regional results, with lower estimated trends in the hierarchical model in the Great Basin, Eastern Tallgrass Prairies, and Western Gulf Coastal Plain. Northern Waterthrush trends were significantly more positive in the hierarchical-model analysis; although none of the BCR-scale estimates were significantly different, the Boreal Hardwood Transition hierarchical model estimate of -1.6% year^-1 was higher than the -4.0% year^-1 from the route regression. Carolina Chickadees showed a significantly more positive trend overall in the hierarchical-model estimates. Although only two BCRs were significantly different, larger estimated trends based on hierarchicalmodel results occurred in 12 of the 13 BCRs for which trends were estimated; only Peninsular Florida had smaller estimated trends with the hierarchical results. Carolina Chickadees had a positive start-up effect (0.02, credible interval = 0.00 to 0.05), which may have caused the route-regression trends of be slightly more negative. Blue-gray Gnatcatchers also showed a pattern of more positive trends in the hierarchical-model analyses. Although only the Appalachian BCR trend differences were significant, 17 of the 19 BCRs with >14 routes had higher estimated trends with the hierarchical-model analysis. The gnatcatchers did not have a positive start-up effect (-0.01, credible interval = -0.04 to 0.02).

Fig. 5.

Annual indices of abundance for Henslow's Sparrow from BBS data, showing the differences in results based on the alternative analyses. Hierarchical-model indices (line with year markers, credible intervals are lines) are scaled year effects from the log-linear model. The route-regression indices (dashed line) are based on average residuals from a composite regression scaled to the mean counts on BBS routes (Sauer and Geissler 1990).

Survey-wide results for Henslow's Sparrow differed dramatically between the analyses. Route-regression results showed extreme declines of -6.9% year^-1, whereas the hierarchical model indicated a nonsignificant -0.6% year^-1 change. Declines in the early years of the BBS are undisputed, but in recent years many investigators (e.g., Herkert 2007) have documented increases in Henslow's Sparrows in many regions. The annual indices associated with the hierarchical-model analysis show recent increases that somewhat mitigated a severe early decline, whereas the routeregression annual indices show a consistent decline (Fig. 5). We attribute the differences in these analyses to an important deficiency of the route-regression approach associated with populations undergoing rapid changes in number. Because the route regression weighted results from routes with more data, the consistently surveyed routes that covered the interval of early population decline had a large influence on the overall trend estimate. Regional weightings may exacerbate this effect, because eastern routes that tend to be well surveyed were likely in the area experiencing the largest declines. The large estimated change over the long term also dominated the residual indices; a positive pattern in the residuals had limited magnitude in relation to the estimated decline. The hierarchical model controlled for this unevenness of coverage because it was fitted using all data simultaneously.

Our definition of trend as a ratio of endpoints reflects our pragmatic view of trend as an interval-specific estimate of change rather than a consistent long-term change in a population. Many species in the BBS have experienced population fluctuations over the interval 1966–2008; estimates of trends for these species are likely to be highly interval-specific. We note that estimates based on alternative definitions of trend such as regressions through the indices can be easily formulated as derived statistics from the hierarchical model.

Species-group Analyses

In this summary, we followed the long-standing practice of grouping birds for summary analysis (e.g., Robbins et al. 1989, Sauer et al. 2008), but we changed species groupings from earlier BBS analyses (e.g., Sauer et al. 2008) to correspond to recently defined groupings (U.S. NABCI Committee 2009). Selection of species groups for summary analysis has generally been controversial, particularly when the groups are given the status of environmental indicators that presumably reflect overall environmental conditions (e.g., Gregory et al. 2005). Two of the intents of the “State of the Birds” report were to stimulate discussion of the use of birds as indicators and to show that well-documented patterns of population decline in some species groups may be indicators of our changing landscapes (U.S. NABCI Committee 2009). These results extend this analysis with additional data from the BBS, reinforce the results described in the report (U.S. NABCI Committee 2009), and form a model for future BBS summary analyses.

Acknowledgments

Supplementary materials for this article, with more data and more comprehensive analysis and summary, are available online. Appendix Sl is available at www.mbr-pwrc.usgs.gov, and Appendix S2 is available at dx.doi.org/10.1525/auk.2010.09220. We thank the many volunteer observers who contribute time, expertise, and energy to support the BBS. C. S. Robbins developed the survey, and the United States and Canadian BBS coordinators have managed and enhanced the survey over the past 45 years. K. Pardieck supplied us with BBS data and commented on the methods and manuscript. J. E. Fallon assisted with data analysis. P. Blank, D. H. Johnson, C. S. Robbins, and two anonymous reviewers provided constructive and sometimes humbling comments on the manuscript. The BBS and State of the Birds websites are at www.mbr-pwrc.usgs.gov/bbs/bbs.html and www.stateofthebirds.org, respectively.

Author notes

Associate Editor: D. H. Johnson