Methodology

INTRODUCTION

The North American Breeding Bird Survey (BBS) is a major source of information on the changes in North American bird populations at broad spatial scales. Since 1966, the BBS has provided trend information at broad geographic scales (range-wide, national, and regional) across much of Canada and the United States for up to 500 species of birds (Hudson et al. 2017, Sauer et al. 2017). BBS data are collected annually by expert volunteers conducting 50, 3-minute point-counts along a roughly 40-km long roadside route, and approximately 5000 routes are surveyed each year (Hudson et al. 2017). The field methods are designed to estimate changes in relative abundance through time by controlling for the effects of survey location, weather, time of day, and season, as well as variations among observers (Sauer et al. 2003). Although the BBS is an excellent source of information on trends, its road-side surveys may not capture all species. Similarly, although the roadside locations are generally representative of the surrounding landscape in most regions (Veech et al. 2017), in some regions it may not represent changes in the landscape far from roads (Van Wilgenburg et al. 2015). It also does not include the information needed to directly model variations in detectability, instead relying on the strictly controlled field methods and statistical adjustments for the most likely sources of variation. These potential sources of bias in BBS trends have been studied and reviewed (e.g., Sauer et al. 2017), and with some care in interpreting the results, BBS trend estimates provide important conservation information (Rosenberg et al. 2017). The BBS was designed to monitor changes in species’ populations over time across broad regions such as the intersection of states/provinces with Bird Conservation Regions (BCRs; Sauer et al. 2003, Link et al. 2020, Smith and Edwards 2021). BCRs are regions of North America, similar to ecoregions (CEC 1997) that share similar ecological characteristics as well as similar bird communities (Bird Studies Canada and North American Bird Conservation Initiative 2014). These regional-scale summaries are critical for identifying and prioritizing species in peril (Government of Canada 2010, IUCN 2012, Rosenberg et al. 2017) and understanding broad-scale patterns of change in North American birds (North American Bird Conservation Initiative Canada 2019, Rosenberg et al. 2019, North American Bird Conservation Initiative 2022).

The BBS dataset can also be analyzed at a finer spatial resolution to complement the regional estimates, and to inform different ecological questions and conservation efforts. Fine-scale estimates of trends should benefit from including the spatial relationships among individual survey locations (e.g., Smith et al. 2024), and allow for visualizing spatial patterns in trends to better understand the local context of range-wide trends and to generate hypotheses on the drivers of population trends. Compared to regional estimates, fine-scale estimates of population trends may provide more useful information for local conservation efforts, as covariates with local effects such as local land cover change, and agricultural practices can be modeled (Thogmartin et al. 2004, Paton et al. 2019, Mirochnitchenko et al. 2021). Many factors influence the relative abundance and trends in bird populations, and they act and interact to induce spatial patterns across a range of spatial scales (Morrison et al. 2010). Factors such as habitat change (Stanton et al. 2018, Betts et al. 2022), biotic factors like prey availability (Drever et al. 2018), or broad-scale patterns in abiotic factors like precipitation, temperature, and phenology (Renfrew et al. 2013, Wilson et al. 2018) can induce spatial patterns in trends or relative abundance and can act across different periods in the species’ annual cycles (Morrison et al. 2010, Wilson et al. 2011). Likewise, conservation actions occur at many scales, from the broad scales of international conventions to the fine scales of an individual wetland (Prairie Habitat Joint Venture 2020).

Including both relative abundance and trend as parameters can improve our understanding of population change, by separately modeling the pattern and covariates that affect the variation in relative abundance in space (e.g., distribution, range dynamics, habitat availability) from those that affect changes in abundance through time (e.g., habitat change, weather, and climate). Earlier fine-scale models for the BBS data did not include the rate of population change (trend) as a parameter in the model (Bled et al. 2013). However, including both relative abundance and trends as separate parameters allows the model to include covariates on each, such as a recent analysis of the effects of forest change on species’ trends on BBS routes (Betts et al. 2022). A spatially explicit hierarchical regression can model both spatial patterns in relative abundance and trend (Ver Hoef et al. 2018, Wright et al. 2021). Here we use the term “hierarchical” model in a general sense to describe models with layered structures where parameters at one level are drawn from distributions and the parameters of those distributions are modeled at higher levels (Cressie et al. 2009, Gelman et al. 2013, Kruschke 2015). Separating these parameters in the model allows for the inclusion of a broader range of covariates (Meehan et al. 2019) to better understand the processes affecting relative abundance (e.g., mean habitat amount or distance to core of a species’ range) and trends (habitat change through time, or differences in climate change effects at northern or southern range limits).

One goal of our work here is to compare two conceptually different approaches to modeling spatial relationships because the spatial locations of BBS observations are not perfectly represented by either of the most common approaches. Spatially explicit models treat sample locations as either discrete areas with neighborhood relationships (Ver Hoef et al. 2018), or points within continuous space (Golding and Purse 2016). It is not obvious whether an area-based or point-based approach better reflects reality for BBS data (Pebesma and Bivand 2023), because the observations from a given BBS route are collected along a transect that is approximately 40 km long. Intrinsic Conditional Autoregressive (iCAR) structures are area-based and model space as a network of polygons with binary neighborhood relationships, e.g., only polygons sharing a border are considered neighbors (Besag and Kooperberg 1995). This area-based approach has been used to model the relatively fine-scale population trends in Christmas Bird Count data (Meehan et al. 2019) and the annual relative abundance of birds using BBS data (Bled et al. 2013). Gaussian Process (GP) models use a continuous measure of distance between points to estimate the correlation of parameters (e.g., trends) between pairs of points and the rate of decrease in the correlation with increasing distance (Golding and Purse 2016). Both approaches are necessary simplifications of the true spatial processes underlying variation in relative abundance and trends among BBS routes. The iCAR approach simplifies the spatial structure by assuming each route represents a discrete area of space (i.e., a polygon surrounding the route; Fig. 1), but the neighboring routes may be separated by a wide range of distances depending on the spatial distribution and spatial density of those routes. The GP approach simplifies spatial relationships by assuming each route represents a point in space, but the measure of intervening distance among routes only applies to the distance between the start points, not to the full transect.

Here we describe and compare four fine-scale BBS models that expand on the published broad-scale BBS models in three ways: (1) they estimate bird population trends and relative abundance at a fine-scale (individual BBS routes); (2) they allow for route-level covariates on the trends and relative abundances; and (3) their output visualizes spatial patterns in both trend and relative abundance. Most of the models share information on relative abundance and trend in a spatially explicit way, and we have included a non-spatial model for comparison. We describe two spatial models that rely on an iCAR structure for the spatial relationships: the first is the iCAR model; and the second is a version of the BYM model, named for Besag, York, and Mollié (Besag et al. 1991), which is identical to the iCAR model but includes an additional random effect to allow extra non-spatial variation in trends. The third spatial model is an isotropic Gaussian Process (GP) model that models spatial relationships using the Euclidean distances among routes. Finally, the fourth model is a non-spatial version that estimates route-level variation in trends and relative abundances as a log-normally distributed random effect. We fit all four models to 71 species and fit one of the spatial models (iCAR) and the non-spatial model to an additional 216 species (details below). We compare the predictive accuracy of models in a leave-future-out cross-validation to assess the benefits of including spatial information and to compare the various approaches to modeling space for the BBS data. Finally, we provide two examples of route-level covariate analyses, to demonstrate the application of these models to conservation and research into understanding the drivers of population change and how those drivers may vary in space and time.

METHODS

Data

We limited most of our analyses to a 15-year period, which we considered short enough that the log-linear slope that represents population trends in these models can be a meaningful summary of the population change (Buckland et al. 2004, Thompson and La Sorte 2008). In effect, 15 years is likely long enough to estimate a rate of change on each route, but also short enough to reduce the likelihood of complex non-linear population patterns. The only exception is the Horned Grebe (Podiceps auritus) covariate example, where we used a 43-year period because the covariate was designed to adjust for annual fluctuations and non-linear patterns in regional moisture/drought cycles (details below). This 15-year period is somewhat arbitrary and for many species or ecological questions, it may be very informative to fit these models (or modifications of these models) to longer or shorter periods of time.

We used 71 species (Table S1) to compare the model predictions and predictive accuracy for all four models and used the Baird’s Sparrow (Centronyx bairdii) to demonstrate model fit and convergence. We chose these 71 species because they have small ranges with relatively few BBS routes, which improved computational efficiency, and yet are also commonly observed during surveys and so provide high-quality data on any given route. We chose the Baird’s Sparrow to demonstrate model fit and convergence because it has some interesting spatial variation in abundance and trends and a very restricted distribution confined to the northern Great Plains region (Fig. 1), which reduces model run-time. Species with large ranges that appear on many routes will increase the computational power required to run the models, increasing the model run-time. Specifically, from 2006 to 2021, these small-range species were observed on 125–400 BBS routes, with at least 600 total observations of the species and an average of at least four observations per route. These thresholds on numbers of routes and observations are effectively arbitrary but balanced the need to have sufficient count-data to estimate parameters well and few enough routes that models would fit relatively quickly. We only compared the fit and predictive accuracy of all four models for these species with fewer than 400 routes, because fitting the GP model requires days or even weeks for species with many routes.

To better understand the benefits of including space in a model for more species, we compared the non-spatial model and one of the spatial models for an additional 216 species that have larger ranges. For these species, we only compared the predictions and predictive accuracy of the non-spatial model to the iCAR model, to reduce our computations. These 216 species were observed on 400 or more BBS routes during 2006–2021 (i.e., species with more routes than the small-range species). We fit these two models to these additional species to assess the more general benefits of including spatial information for more species and for species with populations spread across large ranges that may include many different factors influencing trends and relative abundance (Doser et al. 2024).

Model structure

The four models are hierarchical log-link negative binomial regressions broadly similar to other models commonly applied to BBS data (Sauer and Link 2011, Smith et al. 2014), but modeling trend and relative abundance as route-specific, spatially varying coefficients (Barnett et al. 2021, Thorson et al. 2023, Doser et al. 2024). In all four models, each route has a separate slope (trend) and intercept (relative abundance), but there are no parameters to model yearly fluctuations or non-linear temporal patterns. Therefore, the interpretation of “trend” in these models is limited to this log-linear slope parameter (i.e., a single mean rate of change over the entire modeled time-series).

All of the models have the same basic structure (Fig. 2), varying only in the way the intercept and slope terms were estimated (Fig. 3). In all models, we treated BBS counts as being drawn from a negative binomial distribution (equation 1 and 2, Fig. 2). We included the same observer effects commonly included in BBS trend models (Smith et al. 2024), keeping the observer effects (equation 5, Fig. 2), the inverse dispersion parameter (equation 3, Fig. 2), and the first-year observer effects (equation 4, Fig. 2) in all models. We included the first-year BBS observer parameter because first-year observers have distinct variations in their counts of some species when compared to more experienced observers (Kendall et al. 1996).

Route-level intercepts (alpha terms highlighted in lighter yellow, Fig. 3) and route-level slope parameters (beta terms highlighted in darker green, Fig. 3) were estimated as hierarchical effects, sharing information among routes. Specifically, both the intercepts and the slopes were estimated as an additive combination of a species-mean and a random route-level term (equations 4 and 5, Fig. 3). Three of the models used spatial information to estimate the route-level variation in the intercepts and slopes (i.e., effectively shrinking toward a local mean of neighboring routes), while the fourth non-spatial model estimated them as exchangeable random effects (i.e., shrinking toward a global mean of all routes). To encourage convergence, we constrained many of the random effects in the model, including the spatial route-level parameters, to sum to zero. These constraints often improved model sampling efficiency, but because they are centered on a mean across all routes, they do not affect the interpretation of the final route-level slopes or intercepts (Morris et al. 2019).

To estimate route-level relative abundance, while accounting for variation among observers, we modeled separate intercepts for routes and observers. Using separate observer and route effects has not been commonly included in hierarchical Bayesian models for the BBS (Sauer and Link 2011, Smith et al. 2014, Link et al. 2020, Edwards and Smith 2021), until recently (Betts et al. 2022, Smith et al. 2024). In general, observers and routes are partly associated in the BBS dataset by design as an experimental control for variation among observers (Kendall et al. 1996). However, observers and routes vary in the number of surveys conducted; from 2006 to 2021 more than 69% of surveys were conducted on routes with more than one observer during those 15 years, and 55% of surveys were conducted by observers who have surveyed more than one route.

Finally, we also used an informative prior on the standard deviation of the observer effects (equation 6, Fig. 2), and we ensured that all parameters had converged when fitting the models (details below). We used a half-normal prior on the standard deviation among observers, scaled so that variation among observers is unlikely to result in a six-fold increase, or reduction, in a given species count (equation 3, Fig. 2), and that variation among observers is less than variation among routes. We suggest this prior is reasonable given that BBS observers are highly skilled and familiar with the local bird community (Link and Sauer 1997).

Spatial structures

We fit models with two different approaches to account for spatially explicit relationships among routes: (1) an iCAR structure that treats spatial relationships as a series of discrete neighbors, producing a sparse matrix of adjacencies between pairs of routes; and (2) an isotropic GP model that treats space as the continuous distance between routes, creating a matrix of Euclidean distances between the start locations of each BBS route. To illustrate one difference between the approaches, the GP may consider the relative abundance or trends of two distant routes as effectively independent if the distance is large enough. In contrast, the iCAR structure considers any routes whose polygons share a border as having a very close connection, regardless of polygon size or distance. In some cases, treating two relatively distant routes as close neighbors may be useful if their relative proximity could inform the parameter estimates, but may also introduce error into the estimate of spatial variance (Pebesma and Bivand 2023).

We used a Voronoi tessellation to generate the discrete neighborhood relationships required to support the iCAR model (Ver Hoef et al. 2018, Pebesma and Bivand 2023). iCAR models are often applied to contiguous area-based stratifications, such as regular grids, census regions, or political jurisdictions, which have natural neighborhood relationships defined by their adjacencies (Ver Hoef et al. 2018, Meehan et al. 2019). To generate contiguous discrete spatial units without imposing a regular grid structure, we used a Voronoi tessellation to create contiguous polygons, centered on the start point of each BBS route (Pebesma 2018). We further limited the adjacency matrix to the approximate boundaries of the species’ range by clipping the tessellated surface using the BBS strata where the species occurs (Sauer and Link 2011) and a concave polygon surrounding start locations of all routes with data for that species (Gombin 2023). This clipping ensured that adjacency relationships did not extend beyond the borders of the species’ range and allowed the adjacency matrix to respect large-scale, complex range boundaries (e.g., gaps in forest bird ranges created by the Great Plains).

We adapted functions and code in the Stan probabilistic programming language from the “rethinking” R-package for inclusion in our GP model (McElreath 2023). Like the iCAR approach, we used independent GPs to model the covariance of the intercept parameters and the slope parameters. We estimated the full matrix for between-route distances using functions in the “sf” package for R (Pebesma 2018).

Intrinsic Conditional Autoregressive model: iCAR

We estimated the route-level intercepts and slopes using an iCAR structure, where the parameter for route-r is drawn from a normal distribution, centered on the mean of that parameter’s values in all neighboring routes, with an estimated standard deviation that is proportional to the inverse of the number of neighbors for that route (Morris et al. 2019). Specifically, the route-level variation in the intercept term a random route-level term drawn from a normal distribution centered on the mean of the intercepts for all neighboring routes (equation 6, Fig. 3). The route-level variation in the slope was estimated similarly as a random draw from a normal distribution centered on the mean of the slopes for all neighboring routes (equation 7, Fig. 3).

Besag, York, Mollié iCAR model: BYM

We also used an implementation of the BYM spatial iCAR model (Besag et al. 1991) to estimate route-level slopes. The route intercepts were estimated in the same way as for the iCAR model (equation 6, Fig. 3). The route-level variation in slopes used the same structure as for the iCAR model (equations 8 and 9, Fig. 3), but we added a non-spatial component, estimated as a random effect drawn from a normal distribution with an estimated standard deviation (equation 10, Fig. 3).The additional random effect in the BYM model allows the route-level trend estimates to vary more among neighboring routes, if supported by the data (Besag et al. 1991).

Gaussian Process model: GP

In the GP model, the route-level variation in slope and intercept random terms were estimated as multivariate normal distributions (equations 11 and 12, Fig. 3), with covariance matrices estimated using a squared exponential kernel function (Gelman et al. 2013). The covariance terms for the intercept and slope parameters were separately estimated and are each an exponentially decreasing function of the distance between the routes (equations 13 and 14, Fig. 3).

Non-spatial model

To assess the benefits of the spatially explicit models, we compared the predictions and predictive accuracy of the spatial models to an otherwise identical model that lacked spatial information. This non-spatial model had all the same parameters as the spatial models, except that the route-level variation in the intercepts and slopes were estimated as random effects (equations 15 and 16, Fig. 3).

Convergence

We fit all models using 1000–2000 warmup iterations and an equal number of sampling iterations for each of the four independent chains (or three independent chains for each iteration of cross-validation). We assessed convergence by monitoring for divergent transitions and estimating split-Rhat values and bulk effective sample sizes for all parameters. We considered convergence to have failed if any Rhat was > 1.03 or if any parameter’s effective sample size was < 100 (although the vast majority of parameters had effective sample sizes > 1000 and Rhat < 1.01). For a few GP models that failed to converge, we re-fit the models with the alternative priors described in the supplemental methods.

Model assessment

To assess the benefits of adding spatial information into the model, we compared the 1-year-ahead, leave-future-out (LFO) predictive success of the four models for the 71 species with relatively small ranges (Roberts et al. 2017, Bürkner et al. 2020). In this application of LFO, we fit the model to the first eight years of data (2006–2013; the minimum length of time we considered sufficient for prediction), and used the parameter estimates from this model to predict the counts in the following year (2014). Then we iterated this approach making predictions for the remaining years (2015–2019, and 2021), predicting the held-out data in that year using data for all years prior to fit the model. We could not assess predictive accuracy for the year 2020 because the BBS survey season was canceled because of the public health travel restrictions of spring 2020. We also compared the iCAR spatial model with the non-spatial version of the model using a LFO assessment for an additional 216 species (Table S1). We used the LFO approach to directly test the accuracy of predictions for the next year’s observations.

The cross-validation process generated predictions for every count in the dataset after 2013 and an estimate of the log pointwise predictive density (lppd) of the observed count, given the model and the data in all previous years (Gelman et al. 2014). The lppd is a metric of predictive accuracy that measures how likely it would be to observe the held-out data point (i.e., the observed count in the next year) given the parameter estimates and the model. More specifically, it is the log of the posterior probability (or probability density) of each held-out data point (pointwise), given the model and parameter estimates. Values of lppd are negative (log of probability values, which are < 1.0), and values closer to 0 indicate predictions that are similar to the observed data. For interpretation and visualization, we calculated pairwise differences in lppd between pairs of models for each count and transformed summaries of these lppd differences across all counts into approximate z-scores (mean divided by the standard error of the point-wise differences in lppd). These z-scores summarize the support in the data for each model, accounting for the variation across all observations and providing a consistent scale to summarize pair-wise model comparisons across datasets with different numbers of observations (Link and Sauer 2016). They are an approximation of the test statistic in a paired t-test; e.g., absolute values greater than approximately 2 could be interpreted as a “significant difference” in predictive success, although we have avoided emphasizing arbitrary thresholds.

Route-level covariate examples

Modeling covariates of fine-scale trends and relative abundances is a major benefit of modeling BBS data at the route level. Including covariates, we can investigate conservation and management concerns such as estimating the effects of local habitat change on population trends, or the effects of moisture and climate on local abundance. Further, a fine-scale allows for the use of finer-scale, more localized covariates. To demonstrate this, we present two examples, each including route-level predictors to inform estimates of relative abundance and trend. The first example uses data on habitat suitability from an independent study (Jefferys et al. 2024) on the Rufous Hummingbird (Selasphorus rufus) and models the effect of habitat suitability on relative abundance and change in habitat suitability on trends in BBS data (see Jefferys et al. 2024 and supplemental methods). The second example looks at the effects of annual variation in available habitat—the number of ponds surrounding a BBS route each year in the Prairie Pothole Region (PPR)—on the expected counts of a water bird, the Horned Grebe.

Rufous Hummingbird covariate example

This example model estimates the effect of the amount of suitable habitat around a BBS route on the mean number of birds observed, and the effect of the change in suitable habitat on the change in the number of birds observed through time. This example application is an elaboration of the iCAR route-level trend model, where the route-level intercepts and slopes are additive combinations of two components: (1) one that is a function of a route-level predictor, and (2) one that is a residual component, estimated as a spatially varying component using the iCAR structure (Ver Hoef et al. 2018). As with our previous models, this model used data from the BBS to estimate relative abundance and trend. The route-level predictors are derived from an independent study on Rufous Hummingbirds (Jefferys et al. 2024). In that study, habitat suitability was estimated with a species distribution model using an independent dataset of Rufous Hummingbird observations and data on weather, climate, and landcover and changes in suitability were driven by the loss of early successional forest and warming temperatures in the northeastern regions of the breeding range (Jefferys et al. 2024). In our model, we used the mean habitat suitability from that study across a 15-year period (2006–2021) in a 200- m buffer surrounding each BBS route as a predictor on the intercept (i.e., mean habitat suitability as a predictor on the mean relative abundance on a given route). We used the rate of change in habitat suitability over the same 15-year period within the same buffer as a predictor on the slope (i.e., change in habitat suitability as a predictor on the trend).

We estimated the route-level variation in intercepts and slopes by extending equations 4 and 5 (Fig. 3), to include a component that was a function of the mean habitat suitability for the intercept and the rate of change in habitat suitability for the slope (equations 17 and 18, Fig. 4). The intercepts and slopes for each route were an additive combination of a mean species-level intercept or slope spatially varying residual component (equations 6 and 7, Fig. 4), and a component that was a function of the mean habitat suitability on the route (equations 19 and 20, Fig. 4) or rate of change in habitat suitability on the slope (equations 21 and 22, Fig. 4).

This partitioning of the intercept and slope parameters allows the model to generate two alternative estimates of the mean relative abundance and trend on each route. The full trend (equation 18, Fig. 4) represents the estimated trend on a given route, including the effects of habitat change. The residual trend (i.e., removing the final term from equation 18, Fig. 4) represents an alternate trend if habitat suitability stayed constant on a given route. A similar partitioning of the residual and full estimates of the intercepts is possible, although we did not explore that here.

Horned Grebe covariate example

This example application was an elaboration of the iCAR route-level trend model, where trends and relative abundances are estimated while accounting for the annual variation in climatically dependent habitat. The route-level predictors are derived from a study of the effects of moisture/drought patterns on Horned Grebe (more detail in the supplemental methods), a waterbird species that breeds in small to moderately sized shallow, freshwater ponds (Stedman 2020). To represent annual variation in available habitat for the Horned Grebe in the Canadian Prairie Pothole Region (PPR), we used a dataset collected collaboratively by the U.S. Fish and Wildlife Service (USFWS) and the Canadian Wildlife Service (CWS) across the entire PPR. From this dataset, we used only the Canadian data on the number of ponds (natural or artificial ponds that are flooded seasonally, semi-permanently, and permanently) during the Waterfowl Breeding Population and Habitat Survey (Smith 1995). Annual fluctuations in moisture affect the number of wetlands available, which in turn has a strong influence on waterbird populations that depend on wetlands (Sorenson et al. 1998, Johnson et al. 2005, Roy 2015, Steen et al. 2016). The model was based on the iCAR model and added an additional iCAR component to create a varying-coefficient model on the effects of available habitat on the expected counts during a given survey on a given route. We also fit the same species data to a simpler iCAR model with no covariates to compare the difference in estimated trends with and without accounting for the annual variations in available habitat.

RESULTS

In general, there are clear spatial patterns in the estimated trends and relative abundances from the spatial models, with similar patterns among the three types of spatial models. Those patterns are obscured or completely lacking from the non-spatial version of the model (e.g., the results for Baird’s Sparrow in Figs. 5 and 6). The GP model tended to smooth the spatial pattern in trends more than the iCAR model, which in turn smoothed more than the BYM model (Fig. 5). The spatial smoothing in relative abundance was stronger in both the iCAR and BYM models than the GP model for Baird’s Sparrow (Fig. 5). The covariance in relative abundance of Baird’s Sparrow among routes was effectively 0 at distances of only 100 km, whereas the covariance in trend was relatively strong even at distances > 1000 km (Fig. S1). Predictions of route-level trends had smaller standard errors when including spatial information, and trend precision generally increased with the degree of spatial smoothing (Figs. 6, S2, and S3). For Baird’s Sparrow, all three spatial models had better predictive accuracy than the non-spatial model, with z-scores of pairwise differences between one of the spatial models and the non-spatial model ranging from 2.7 to 3.3 (Fig. S4). The iCAR model had better predictive accuracy than the BYM model (z-score of the difference = 3.8; Fig. 7), and there was little difference in predictive accuracy between the iCAR and GP models (z-score difference = -0.51; Fig. 7).

The leave-future-out (LFO) cross-validation shows that the iCAR and GP models out-perform the non-spatial model (i.e., more accurately predicted the next year’s data) for almost all the 71 small-range species (Fig. 7 and Fig. S4). Out of the spatial models, the GP model had the highest predictive accuracy for the greatest number of species, followed by the iCAR model, and the BYM model had the lowest predictive accuracy. The BYM model was the only spatial model that had clearly lower predictive accuracy than the non-spatial model for any species (i.e., four species for which the z-score difference is < -2, Fig. 7 and Fig. S4). The iCAR model and the GP model had comparable predictive accuracy for many species (most values between -2 and +2 in the iCAR-GP comparison of Fig. 7); 69% (49 of 71 species) of the species were better predicted by the GP model (negative values Fig. 7) and the remaining species were better predicted by the iCAR model (positive values in Fig. 7). When including the additional 216 species for which fitting the GP model was prohibitively time-consuming (days or even weeks are required for convergence for a given species), the iCAR spatial model had higher predictive accuracy than the non-spatial model for 283 of 287 species, and predictive accuracy was very similar for the remaining four (Fig. 8).

The iCAR model generated trend prediction maps with clear spatial patterns that likely relate to spatially dependent variation in processes affecting populations (Fig. 9). These patterns are not evident in predictions from an identical model without spatial information (Fig. 9). The spatial patterns in route-level trends vary widely among these species and among the others we tested (Figs. S2 and S3), suggesting varied drivers of population change across the continent and among species.

In general, the iCAR and GP models were comparable in predictive accuracy for the 71 small-range species we analyzed (Fig. 7). In addition, the spatial patterns in predicted trends were very similar between these two models, even for species where the predictive accuracy differed between the models (Fig. 10). For example, the GP model had higher predictive accuracy than the iCAR model (z-score difference = -4.3, Fig. S4) for Canyon Towhee (Melozone fusca), but the opposite was true for Western Bluebird (Sialia mexicana; z-score difference = 2.3, Fig. S4). Regardless, the spatial pattern in predicted trends between the two models is quite similar for both species (Fig. 10). For both species, and in general, the GP model trend estimates had narrower credible intervals (higher estimated precision) than the iCAR model (Figs. S5 and S6). Precision of the iCAR trend estimates also showed a clear relationship to the number of neighbors for any given route, in that routes with few neighbors (on the edges of the species’ range) were much less precise than estimates in the core of the species’ range (Fig. S6).

Including habitat suitability in the Rufous Hummingbird population model affected estimates of route-level relative abundance and improved estimates of the spatial pattern in trends (Fig. 11). However, much of the overall decline was not related to covariates describing route-level habitat change, as the negative population trends across the species’ range remained after removing the effects of local habitat change covariates (right panel, Fig. 11). The effect of habitat suitability on mean relative abundance was strong and positive (ρ_α = 3 [95% CI 2.2:3.8]), such that routes with higher overall habitat suitability had higher mean counts. From 2006 to 2021, the Rufous Hummingbird’s overall population declined steeply, decreasing by approximately -43% (95% CI -52:-33). There was a positive effect of change in habitat suitability on trends, such that routes with habitat loss had more negative population trends ρ_β = 0.025 (95% CI 0.003:0.047). Trends were negative across the species’ range, but most negative in the coastal regions where the habitat has changed the most and where the species is also most abundant (left panel, Fig. 11, and Fig. S7). The change in habitat suitability affected the spatial patterns in trend (Fig. 11), with the greater loss of habitat in the coastal regions (Fig. S7) accounting for most of the increased rates of decline in the core of the species’ range. The residual trend component alone does not show the same coastal-decline pattern (right panel, Fig.11).

Annual variation in the number of ponds around BBS routes affected the overall rate of population change in Horned Grebes and showed a spatial relationship (Fig. 12). In a model including the annual pond variation, the Horned Grebe population declined overall at a rate of -1.9%/year from 1975 to 2017. After removing the effect of annual pond variation, the long-term rate of decline was -2.2%/year. The effect of annual fluctuations in the number of ponds was positive across the region: the mean value of ρʹ = 0.42 (95% CI 0.29:0.55), but there was also a spatial gradient in intensity. The effect of the number of ponds per year was strongest in the northwest part of the Prairies (Fig. 12).

DISCUSSION

Spatially explicit, route-level models are useful for visualizing fine spatial patterns at scales more relevant to local conservation, understanding the drivers of population change, and estimating the effects of covariates on relative abundance and trends (e.g., Betts et al. 2022). At this fine spatial scale, incorporating spatial information improves the models’ predictions of future data. This improvement is particularly clear for both the iCAR and the GP models, where the spatial models had higher accuracy for out-of-sample predictions than the non-spatial model for almost every species we compared. Fine spatial patterns in trend estimates across a species’ range are useful for generating hypotheses on the ecological drivers of population change. Route-level models also allow for the incorporation of local habitat covariates on relative abundance and trend at fine scales, which is important as some covariates affect bird populations at scales much smaller than the strata often used for broad-scale analyses, such as BCRs or states/provinces/territories (Thogmartin et al. 2004, Paton et al. 2019, Monroe et al. 2022). Route-level patterns are also useful in guiding conservation and/or further monitoring efforts, such as identifying small areas for conservation purposes or diverging population trends within management areas (i.e., strata or BCR).

These route-level, spatial models generate smoothed patterns of variation in population trends across a species’ range, which will greatly facilitate hypothesis generation and direct investigation to better understand the drivers of population change similar to Fink et al. (2023a). For example, the spatial models show relatively smooth patterns in Baird’s Sparrow trends across the species’ range (Fig. 4), which are not evident in the simpler, non-spatial model. In the spatial models, Baird’s Sparrow has increased in the west and decreased in the eastern portion of its range. This longitudinal pattern may suggest hypotheses related to spatial variation in factors such as climate, or habitat amount, which could then be directly tested by incorporating covariates representing these factors into a subsequent model. Similarly, the complex spatial patterns in the trends of Hairy Woodpecker (Dryobates villosus, Fig. 8) show some latitudinal variation in trends in the west that is not as clear in the east, suggesting that there may be distinct processes related to latitude driving trends in these two regions. Comparisons of these patterns among species may be particularly informative. For example, the similar southeast to northwest gradients in trends for Canyon Towhee and Western Bluebird may suggest some similarity in the underlying drivers of population change (Fig. 9). These observations illustrate the types of hypothesis-generating that these fine-scale, spatially explicit models can help generate.

All three of the spatial models (iCAR, GP, and BYM) generated broadly similar spatial patterns in route-level trends for the subset of species we compared (Fig. 5 and Fig. S2). For the species in this study, there is little support for the extra variation in route-level trends in the BYM model, given it had lower predictive accuracy than the simpler iCAR model in all cases. The iCAR structure outperformed the GP models for 31% of the species and is more computationally efficient. Overall, the GP model outperformed the iCAR model for most (69%) of the species we compared. The GP model also produces smoother spatial patterns in population trends than the other spatial models and for some, the difference is striking (e.g., Black-throated Gray Warbler Setophaga nigrescens, California Quail Callipepla californica, and the Golden-winged Warbler Vermivora chrysoptera in Fig. S3). For the first two species, the GP outperformed the iCAR for accuracy, while for the third species, the iCAR was better (Fig. S4).

Though the GP model had better predictive accuracy for most species, the best spatial structure to use will depend on the species and the goals of a given study. Until GP models become more efficient to implement, the iCAR structure may be preferable for larger datasets (e.g., broad-ranging species and or longer time-series). The iCAR structure may also provide more direct control to model discontinuities in the spatial relationships, such as complex range boundaries (Ver Hoef et al. 2018, Pebesma and Bivand 2023), because there are many ways to define neighborhood relationships (Freni-Sterrantino et al. 2018). A species with limited dispersal may be particularly sensitive to the Euclidean distance between points and therefore better modeled with the GP, but the simplification of space using the iCAR structure may be sufficient for most wide-ranging migratory birds. For example, for some species, there are routes that are separated from most other routes by relatively large distances. These “isolated” routes are treated very differently by the iCAR and GP models: they are considered close neighbors in the iCAR model regardless of the distance between them, whereas in the GP model, the large distance between routes reduces their correlation with their nearest neighbors. Interestingly, when we compared the predictive accuracy between the two models for these “isolated” routes (nearest neighboring route with data was greater than 200 km away), the iCAR tended to outperform the GP (Fig. S8). Therefore, a more accurate representation of the long distances separating these isolated routes in the GP model does not necessarily result in more accurate predictions, and when combined with the GP’s computational load, it may be more effective to treat space as a series of relative spatial relationships using the iCAR structure.

These route-level BBS models provide many opportunities for further comparisons, applications, and elaborations. Fine-scale estimates could be summarized across species and within regions, such as summaries of the spatial patterns in grassland bird trends or summaries for a given species within BCRs or states/provinces/territories and compared to estimates from models fit at those broader spatial scales. The spatial patterns in trend estimates also allow for comparison of BBS data to other fine-grained maps of species trend and relative abundance, such as eBird (Sullivan et al. 2014, Fink et al. 2023b) or the Integrated Monitoring in Bird Conservation Regions (IMBCR) program (Pavlacky et al. 2017). Comparison of trend estimates between the two programs for the same species and periods of time could provide useful validation of and or help understand differences between the two sources of information. Similarly, there are many possible avenues to integrate information across programs for a given period (e.g., recent trends) or through time (e.g., long-term information from the BBS with more recent information from eBird and/or IMBCR). Through data integration, we can overcome some of the limitations of the BBS, such as the lack of detectability data (Edwards et al. 2023a), and road-side survey, while retaining the program’s benefits of a long time-series with a structured spatial design and consistent sampling through time.

Separating the route-level intercepts from the observer-level intercepts allows us to better model patterns in relative abundance. It should also allow for improved modeling of the variation among observers. Although many previous BBS analyses have combined observer and route effects (Link et al. 2020, Smith and Edwards 2021), doing so contributes some of the variation in relative abundance among routes to observer variation, which is effectively sampling noise. The model will struggle to estimate separate intercepts for observers and routes in situations where there are few data to inform the estimates (e.g., cases where a route has only been surveyed by one observer). However, we suggest that if a model has sufficient data to estimate these parameters separately, however weakly, it is preferable to a model that confounds the variation in relative abundance among routes from the sampling noise of observer variation. This separation of the observer from route effects is improved by the hierarchical structure of the models, inclusion of spatial information, weakly informative priors, and the improved efficiency of HMC algorithms over the earlier Gibbs sampling algorithms. Although we were motivated by our desire to directly model route-level relative abundance, this approach is equally applicable to other BBS analyses (Smith et al. 2024), at any scale, and is included in the models in the R-package bbsBayes2 (Edwards et al. 2023b).

In both covariate examples, incorporating spatial covariates into the trend analyses tested hypotheses related to the drivers of population change and helped identify specific areas for further research and conservation action. For the Rufous Hummingbird, the model shows higher mean relative abundance on routes with higher habitat suitability and positive effects of the change in habitat suitability on the species’ trend (more negative trends on routes where habitat has decreased). These example findings coincide with a recent study that found that the survival rate of Rufous Hummingbirds was negatively affected by high human population density (English et al. 2024), where there is likely less habitat. Interestingly, it also shows that during this period, the variation among routes in habitat change does not account for all of the decline in the species’ population (Fig. 11, and Fig. S7), suggesting that factors other than local habitat or factors acting outside of the breeding range may be driving the overall decline. However, covariates other than habitat suitability could represent local habitat better for the Rufous Hummingbird and could result in a different relationship between local habitat and relative abundance. For the Horned Grebe, the effect of annual fluctuations in available wetland habitat (the number of ponds) is positive overall and varies in magnitude across the species’ range. The effect is strongest in the western Prairies where the effects of drought are often strongest (Johnson et al. 2005, Millett et al. 2009, Roy 2015). These results highlight the importance of continued investment in wetland conservation programs for waterbird populations breeding in the Prairie Potholes Region, and the vulnerability of these species to climate change because their breeding habitat is highly sensitive to climatic conditions.

The structure of the BBS, designed for monitoring temporal changes in bird populations, allows for the efficient estimation of fine-scale patterns in trends and the effects of local drivers of those trends, provided the survey design and the model sufficiently account for potential changes in detectability. Unmodeled changes in detectability of birds through time could explain some of the spatial patterns in trends from these models, if the changes in detectability were coincident across many BBS routes in a particular region. For example, changes in vegetation that affect detectability, such as forest loss adjacent to the roads where surveys occur, could bias estimates of trends from BBS observations or bias estimates of the effects of that changing vegetation on bird populations. The BBS field methods control many factors that are known to affect detectability, including weather, season, time of day, and among observer variation. Other potential sources of bias in BBS trends include changes in phenology, changes in traffic during surveys, among observer variation, and within observer variation, although so far each of these when studied has been shown to have minor effects on trend estimates (Kendall et al. 1996, Griffith et al. 2010, English et al. 2021), or can be statistically adjusted in the models (Sauer et al. 1994). There is ongoing work to further explore the potential bias in trends due to observer aging (Farmer et al. 2014, U.S. Geological Survey and Canadian Wildlife Service 2020). The models here could also be modified to integrate BBS observations with additional data that could support adjustments for possible changes in detectability (e.g., Edwards et al. 2023a), or additional data that could directly estimate changes in detectability through time (Pavlacky et al. 2017, Zhao et al. 2024).

Fine-scale models can also be used to inform different scales of decisions and communities. Decisions on land use for industries such as agriculture, forestry, and housing are often made at fine scales (Sodhi et al. 2011, Malek et al. 2019). Likewise, habitat protection and restoration by community organizations, municipal governments, and non-governmental organizations occur at fine scales (Sheppard 2005, Aronson et al. 2017). For example, the Horned Grebe covariate analysis confirms the vulnerability of waterbird species in the northwestern Prairie Potholes Region and supports a current initiative to protect critical shallow wetlands in the region (Prairie Habitat Joint Venture 2020). Community support is important for the success of conservation initiatives (Berkes 2004, Bennett and Dearden 2014), and so providing estimates at scales relevant to communities may increase community support for conservation and encourage a feeling of stewardship. Further, routes are a relevant scale for the volunteer observers dedicated to the BBS, with the average BBS volunteer participating for 12 years. Producing estimates at a route-level allows volunteers to see the direct results of their efforts over the years, a large motivator for many citizen science volunteers (Phillips et al. 2019).

We hope that the models and examples here will facilitate greater use of the BBS data, providing new ways to explore the spatiotemporal patterns in relative abundance and trends, and new tools with which to better understand the drivers of those patterns. The long-term information from the BBS provides a priceless benchmark against which to calibrate an otherwise shifting ecological baseline. In addition to its monitoring role, the structured design of the BBS also provides excellent data to study the drivers and correlates of population change using tractable and interpretable models such as these.

LITERATURE CITED

Aronson, M. F., C. A. Lepczyk, K. L. Evans, M. A. Goddard, S. B. Lerman, J. S. MacIvor, C. H. Nilon, and T. Vargo. 2017. Biodiversity in the city: key challenges for urban green space management. Frontiers in Ecology and the Environment 15:189-196. https://doi.org/10.1002/fee.1480

Barnett, L. A. K., E. J. Ward, and S. C. Anderson. 2021. Improving estimates of species distribution change by incorporating local trends. Ecography 44:427-439. https://doi.org/10.1111/ecog.05176

Bennett, N. J., and P. Dearden. 2014. Why local people do not support conservation: community perceptions of marine protected area livelihood impacts, governance and management in Thailand. Marine Policy 44:107-116. https://doi.org/10.1016/j.marpol.2013.08.017

Berkes, F. 2004. Rethinking community-based conservation. Conservation Biology 18:621-630. https://doi.org/10.1111/j.1523-1739.2004.00077.x

Besag, J., and C. Kooperberg. 1995. On conditional and intrinsic autoregressions. Biometrika 82:733-746. https://doi.org/10.1093/biomet/82.4.733

Besag, J., J. York, and A. Mollié. 1991. Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics 43:1-20. https://doi.org/10.1007/BF00116466

Betts, M. G., Z. Yang, A. S. Hadley, A. C. Smith, J. S. Rousseau, J. M. Northrup, J. J. Nocera, N. Gorelick, and B. D. Gerber. 2022. Forest degradation drives widespread avian habitat and population declines. Nature Ecology & Evolution 6:709-719. https://doi.org/10.1038/s41559-022-01737-8

Bird Studies Canada, and North American Bird Conservation Initiative. 2014. Bird Conservation Regions. Published by Bird Studies Canada on behalf of the North American Bird Conservation Initiative. https://birdscanada.org/bird-science/nabci-bird-conservation-regions

Bled, F., J. Sauer, K. Pardieck, P. Doherty, and J. A. Royle. 2013. Modeling trends from North American Breeding Bird Survey data: a spatially explicit approach. PLoS ONE 8:e81867. https://doi.org/10.1371/journal.pone.0081867

Buckland, S. T., K. B. Newman, L. Thomas, and N. B. Koesters. 2004. State-space models for the dynamics of wild animal populations. Ecological Modelling 171:157-175. https://doi.org/10.1016/j.ecolmodel.2003.08.002

Bürkner, P.-C., J. Gabry, and A. Vehtari. 2020. Approximate leave-future-out cross-validation for Bayesian time series models. Journal of Statistical Computation and Simulation 90:2499-2523. https://doi.org/10.1080/00949655.2020.1783262

Commission for Environmental Cooperation (CEC). 1997. Ecological Regions of North America: Toward a Common Perspective. CEC, Montréal, Québec, Canada.

Cressie, N., C. A. Calder, J. S. Clark, J. M. V. Hoef, and C. K. Wikle. 2009. Accounting for uncertainty in ecological analysis: the strengths and limitations of hierarchical statistical modeling. Ecological Applications 19:553-570. https://doi.org/10.1890/07-0744.1

Doser, J. W., M. Kéry, S. P. Saunders, A. O. Finley, B. L. Bateman, J. Grand, S. Reault, A. S. Weed, and E. F. Zipkin. 2024. Guidelines for the use of spatially varying coefficients in species distribution models. Global Ecology and Biogeography 33:e13814. https://doi.org/10.1111/geb.13814

Drever, M. C., A. C. Smith, L. A. Venier, D. J. H. Sleep, and D. A. MacLean. 2018. Cross-scale effects of spruce budworm outbreaks on boreal warblers in eastern Canada. Ecology and Evolution 8:7334-7345. https://doi.org/10.1002/ece3.4244

Edwards, B. P. M., and A. C. Smith. 2021. bbsBayes: an R Package for Hierarchical Bayesian Analysis of North American Breeding Bird Survey Data. Journal of Open Research Software 9:19. https://doi.org/10.5334/jors.329

Edwards, B. P. M., A. C. Smith, T. D. S. Docherty, M. A. Gahbauer, C. R. Gillespie, A. R. Grinde, T. Harmer, D. T. Iles, S. M. Matsuoka, N. L. Michel, A. Murray, G. J. Niemi, J. Pasher, D. C. Pavlacky Jr, B. G. Robinson, T. B. Ryder, P. Sólymos, D. Stralberg, and E. J. Zlonis. 2023a. Point count offsets for estimating population sizes of north American landbirds. Ibis 165:482-503. https://doi.org/10.1111/ibi.13169

Edwards, B. P. M., A. C. Smith, and S. LaZerte. 2023b. bbsBayes2: Hierarchical Bayesian Analysis of North American BBS Data.

English, S. G., C. A. Bishop, S. Wilson, and A. C. Smith. 2021. Current contrasting population trends among North American hummingbirds. Scientific Reports 11:18369. https://doi.org/10.1038/s41598-021-97889-x

English, S. G., S. Wilson, Q. Zhao, C. A. Bishop, and A. J. Moran. 2024. Demographic mechanisms and anthropogenic drivers of contrasting population dynamics of hummingbirds. Biological Conservation 289:110415. https://doi.org/10.1016/j.biocon.2023.110415

Farmer, R. G., M. L. Leonard, J. E. Mills Flemming, and S. C. Anderson. 2014. Observer aging and long-term avian survey data quality. Ecology and Evolution 4:2563-2576. https://doi.org/10.1002/ece3.1101

Fink, D., T. Auer, A. Johnson, M. Strimas-Mackey, S. Ligocki, O. Robinson, W. Hochachka, L. Jaromczyk, C. Crowley, K. Dunham, A. Stillman, I. Davies, A. Rodewald, V. Ruiz-Gutierrez, and C. Wood. 2023b. eBird status and trends, data version: 2022; Released: 2023. https://doi.org/10.2173/ebirdst.2022

Fink, D., A. Johnston, M. Strimas-Mackey, T. Auer, W. M. Hochachka, S. Ligocki, L. Oldham Jaromczyk, O. Robinson, C. Wood, S. Kelling, and A. D. Rodewald. 2023a. A double machine learning trend model for citizen science data. Methods in Ecology and Evolution 14:2435-2448. https://doi.org/10.1111/2041-210X.14186

Freni-Sterrantino, A., M. Ventrucci, and H. Rue. 2018. A note on intrinsic conditional autoregressive models for disconnected graphs. Spatial and Spatio-temporal Epidemiology 26:25-34. https://doi.org/10.1016/j.sste.2018.04.002

Gelman, A., J. Carlin B., H. S. Stern, D. Dunson B., A. Vehtari, and D. B. Rubin. 2013. Bayesian data analysis. Third edition. Chapman and Hall/CRC, New York, New York, USA. https://doi.org/10.1201/b16018

Gelman, A., J. Hwang, and A. Vehtari. 2014. Understanding predictive information criteria for Bayesian models. Statistics and Computing 24:997-1016. https://doi.org/10.1007/s11222-013-9416-2

Golding, N., and B. V. Purse. 2016. Fast and flexible Bayesian species distribution modelling using Gaussian processes. Methods in Ecology and Evolution 7:598-608. https://doi.org/10.1111/2041-210X.12523

Gombin, J. 2023, May 27. Concaveman package. JavaScript.

Government of Canada. 2010. COSEWIC definitions associated with quantitative criteria. Canadian Wildlife Service, Gatineau, Québec, Canada. https://cosewic.ca/index.php/en/assessment-process/cosewic-assessment-process-categories-and-guidelines/quantitative-criteria-definitions.html

Griffith, E. H., J. R. Sauer, and J. A. Royle. 2010. Traffic effects on bird counts on North American Breeding Bird Survey routes. Auk 127:387-393. https://doi.org/10.1525/auk.2009.09056

Hudson, M.-A. R., C. M. Francis, K. J. Campbell, C. M. Downes, A. C. Smith, and K. L. Pardieck. 2017. The role of the North American Breeding Bird Survey in conservation. Condor 119:526-545. https://doi.org/10.1650/CONDOR-17-62.1

International Union for Conservation of Nature (IUCN). 2012. IUCN Red List categories and criteria: Version 3.1. Second edition. IUCN, Gland, Switzerland.

Jefferys, K. M., M. G. Betts, W. D. Robinson, J. R. F. Curtis, T. A. Hallman, A. C. Smith, C. Strevens, and J. Aguirre-Gutiérrez. 2024. Breeding habitat loss linked to declines in Rufous Hummingbirds. Avian Conservation and Ecology 19(2):2. https://doi.org/10.5751/ACE-02681-190202

Johnson, W. C., B. V. Millett, T. Gilmanov, R. A. Voldseth, G. R. Guntenspergen, and D. E. Naugle. 2005. Vulnerability of Northern Prairie wetlands to climate change. BioScience 55:863-872. https://doi.org/10.1641/0006-3568(2005)055[0863:VONPWT]2.0.CO;2

Kendall, W. L., B. G. Peterjohn, and J. R. Sauer. 1996. First-time observer effects in the North American Breeding Bird Survey. Auk 113:823-829. https://doi.org/10.2307/4088860

Kruschke, J. K. 2015. Doing Bayesian data analysis: a tutorial with R, JAGS, and Stan. Second edition. Academic, Boston, Massachusetts, USA.

Link, W. A., and J. R. Sauer. 1997. New approaches to the analysis of population trends in land birds: comment. Ecology 78:2632-2634. https://doi.org/10.1890/0012-9658(1997)078[2632:NATTAO]2.0.CO;2

Link, W. A., and J. R. Sauer. 2016. Bayesian cross-validation for model evaluation and selection, with application to the North American Breeding Bird Survey. Ecology 97:1746-1758. https://doi.org/10.1890/15-1286.1

Link, W. A., J. R. Sauer, and D. K. Niven. 2020. Model selection for the North American Breeding Bird Survey. Ecological Applications 30:e02137. https://doi.org/10.1002/eap.2137

Malek, Ž., B. Douw, J. V. Vliet, E. H. Van Der Zanden, and P. H. Verburg. 2019. Local land-use decision-making in a global context. Environmental Research Letters 14:083006. https://doi.org/10.1088/1748-9326/ab309e

McElreath, R. 2023. rethinking. R. https://github.com/rmcelreath/rethinking

Meehan, T. D., N. L. Michel, and H. Rue. 2019. Spatial modeling of Audubon Christmas Bird Counts reveals fine-scale patterns and drivers of relative abundance trends. Ecosphere 10:e02707. https://doi.org/10.1002/ecs2.2707

Millett, B., W. C. Johnson, and G. Guntenspergen. 2009. Climate trends of the North American prairie pothole region 1906-2000. Climatic Change 93:243-267. https://doi.org/10.1007/s10584-008-9543-5

Mirochnitchenko, N. A., E. F. Stuber, and J. J. Fontaine. 2021. Biodiversity scale-dependence and opposing multi-level correlations underlie differences among taxonomic, phylogenetic and functional diversity. Journal of Biogeography 48:2989-3003. https://doi.org/10.1111/jbi.14248

Monroe, A. P., J. A. Heinrichs, A. L. Whipple, M. S. O’Donnell, D. R. Edmunds, and C. L. Aldridge. 2022. Spatial scale selection for informing species conservation in a changing landscape. Ecosphere 13:e4320. https://doi.org/10.1002/ecs2.4320

Morris, M., K. Wheeler-Martin, D. Simpson, S. J. Mooney, A. Gelman, and C. DiMaggio. 2019. Bayesian hierarchical spatial models: implementing the Besag York Mollié model in stan. Spatial and Spatio-temporal Epidemiology 31:100301. https://doi.org/10.1016/j.sste.2019.100301

Morrison, C. A., R. A. Robinson, J. A. Clark, and J. A. Gill. 2010. Spatial and temporal variation in population trends in a long-distance migratory bird. Diversity and Distributions 16:620-627. https://doi.org/10.1111/j.1472-4642.2010.00663.x

North American Bird Conservation Initiative. 2022. State of the birds 2022. https://www.stateofthebirds.org/2022/

North American Bird Conservation Initiative Canada. 2019. The state of Canada’s birds, 2019. Environment and Climate Change Canada, Ottawa, Canada. https://publications.gc.ca/collections/collection_2020/eccc/cw66/CW66-312-2019-eng.pdf

Paton, G. D., A. V. Shoffner, A. M. Wilson, and S. A. Gagné. 2019. The traits that predict the magnitude and spatial scale of forest bird responses to urbanization intensity. PLOS ONE 14:e0220120. https://doi.org/10.1371/journal.pone.0220120

Pavlacky, D. C., P. M. Lukacs, J. A. Blakesley, R. C. Skorkowsky, D. S. Klute, B. A. Hahn, V. J. Dreitz, T. L. George, and D. J. Hanni. 2017. A statistically rigorous sampling design to integrate avian monitoring and management within Bird Conservation Regions. PLoS ONE 12:e0185924. https://doi.org/10.1371/journal.pone.0185924

Pebesma, E. 2018. Simple features for R: standardized support for spatial vector data. R Journal 10:439-446. https://doi.org/10.32614/RJ-2018-009

Pebesma, E., and R. Bivand. 2023. Spatial data science: with applications in R. First edition. Chapman and Hall/CRC, New York, New York, USA. https://doi.org/10.1201/9780429459016

Phillips, T. B., H. L. Ballard, B. V. Lewenstein, and R. Bonney. 2019. Engagement in science through citizen science: moving beyond data collection. Science Education 103:665-690. https://doi.org/10.1002/sce.21501

Prairie Habitat Joint Venture. 2020. Prairie habitat joint venture: the Prairie Parklands Implementation Plan 2013-2020.

Renfrew, R. B., D. Kim, N. Perlut, J. Smith, J. Fox, and P. P. Marra. 2013. Phenological matching across hemispheres in a long-distance migratory bird. Diversity and Distributions 19:1008-1019. https://doi.org/10.1111/ddi.12080

Roberts, D. R., V. Bahn, S. Ciuti, M. S. Boyce, J. Elith, G. Guillera-Arroita, S. Hauenstein, J. J. Lahoz-Monfort, B. Schröder, W. Thuiller, D. I. Warton, B. A. Wintle, F. Hartig, and C. F. Dormann. 2017. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography 40:913-929. https://doi.org/10.1111/ecog.02881

Rosenberg, K. V., P. J. Blancher, J. C. Stanton, and A. O. Panjabi. 2017. Use of North American Breeding Bird Survey data in avian conservation assessments. Condor 119:594-606. https://doi.org/10.1650/CONDOR-17-57.1

Rosenberg, K. V., A. M. Dokter, P. J. Blancher, J. R. Sauer, A. C. Smith, P. A. Smith, J. C. Stanton, A. Panjabi, L. Helft, M. Parr, and P. P. Marra. 2019. Decline of the North American avifauna. Science 366:120-124. https://doi.org/10.1126/science.aaw1313

Roy, C. 2015. Quantifying geographic variation in the climatic drivers of midcontinent wetlands with a spatially varying coefficient model. PLoS ONE 10:e0126961. https://doi.org/10.1371/journal.pone.0126961

Sauer, J. R., J. E. Fallon, and R. Johnson. 2003. Use of North American Breeding Bird Survey data to estimate population change for bird conservation regions. Journal of Wildlife Management 67:372-389. https://doi.org/10.2307/3802778

Sauer, J. R., and W. A. Link. 2011. Analysis of the North American Breeding Bird Survey using hierarchical models. Auk 128:87-98. https://doi.org/10.1525/auk.2010.09220

Sauer, J. R., K. L. Pardieck, D. J. Ziolkowski, A. C. Smith, M.-A. R. Hudson, V. Rodriguez, H. Berlanga, D. K. Niven, and W. A. Link. 2017. The first 50 years of the North American Breeding Bird Survey. Condor 119:576-593. https://doi.org/10.1650/CONDOR-17-83.1

Sauer, J. R., B. G. Peterjohn, and W. A. Link. 1994. Observer differences in the North American Breeding Bird Survey. Auk 111:50-62. https://doi.org/10.2307/4088504

Sheppard, S. R. 2005. Participatory decision support for sustainable forest management: a framework for planning with local communities at the landscape level in Canada. Canadian Journal of Forest Research 35:1515-1526. https://doi.org/10.1139/x05-084

Smith, A. C., A. D. Binley, L. Daly, B. P. M. Edwards, D. Ethier, B. Frei, D. Iles, T. D. Meehan, N. L. Michel, and P. A. Smith. 2024. Spatially explicit Bayesian hierarchical models improve estimates of avian population status and trends. Ornithological Applications 126:duad056. https://doi.org/10.1093/ornithapp/duad056

Smith, A. C., and B. P. M. Edwards. 2021. North American Breeding Bird Survey status and trend estimates to inform a wide range of conservation needs, using a flexible Bayesian hierarchical generalized additive model. Condor 123:duaa065. https://doi.org/10.1093/ornithapp/duaa065

Smith, A. C., M.-A. R. Hudson, C. Downes, and C. M. Francis. 2014. Estimating breeding bird survey trends and annual indices for Canada: how do the new hierarchical Bayesian estimates differ from previous estimates? Canadian Field-Naturalist 128:119-134. https://doi.org/10.22621/cfn.v128i2.1565

Smith, G. W. 1995. A critical review of the aerial and ground surveys of breeding waterfowl in North America. U.S. Department of the Interior, National Biological Service, Washington, D.C., USA.

Sodhi, N. S., R. Butler, W. F. Laurance, and L. Gibson. 2011. Conservation successes at micro-, meso- and macroscales. Trends in Ecology & Evolution 26:585-594. https://doi.org/10.1016/j.tree.2011.07.002

Sorenson, L. G., R. Goldberg, T. L. Root, and M. G. Anderson. 1998. Potential effects of global warming on waterfowl populations breeding in the Northern Great Plains. Climatic Change 40:343-369. https://doi.org/10.1023/A:1005441608819

Stanton, R. L., C. A. Morrissey, and R. G. Clark. 2018. Analysis of trends and agricultural drivers of farmland bird declines in North America: a review. Agriculture, Ecosystems & Environment 254:244-254. https://doi.org/10.1016/j.agee.2017.11.028

Stedman, S. J. 2020. Horned Grebe (Podiceps auritus), version 1.0. In S. M. Billerman, editor. Birds of the world. Cornell Lab of Ornithology, Ithaca, New York, USA. https://doi.org/10.2173/bow.horgre.01

Steen, V. A., S. K. Skagen, and C. P. Melcher. 2016. Implications of climate change for wetland-dependent birds in the Prairie Pothole Region. Wetlands 36:445-459. https://doi.org/10.1007/s13157-016-0791-2

Sullivan, B. L., J. L. Aycrigg, J. H. Barry, R. E. Bonney, N. Bruns, C. B. Cooper, T. Damoulas, A. A. Dhondt, T. Dietterich, A. Farnsworth, D. Fink, J. W. Fitzpatrick, T. Fredericks, J. Gerbracht, C. Gomes, W. M. Hochachka, M. J. Iliff, C. Lagoze, F. A. La Sorte, M. Merrifield, W. Morris, T. B. Phillips, M. Reynolds, A. D. Rodewald, K. V. Rosenberg, N. M. Trautmann, A. Wiggins, D. W. Winkler, W. Weng-Keen, C. L. Wood, J. Yu, and S. Kelling. 2014. The eBird enterprise: an integrated approach to development and application of citizen science. Biological Conservation 169:31-40. https://doi.org/10.1016/j.biocon.2013.11.003

Thogmartin, W. E., J. R. Sauer, and M. G. Knutson. 2004. A hierarchical spatial model of avian abundance with application to Cerulean Warblers. Ecological Applications 14:1766-1779. https://doi.org/10.1890/03-5247

Thompson, F. R., and F. A. La Sorte. 2008. Comparison of methods for estimating bird abundance and trends from historical count data. Journal of Wildlife Management 72:1674-1682. https://doi.org/10.2193/2008-135

Thorson, J. T., C. L. Barnes, S. T. Friedman, J. L. Morano, and M. C. Siple. 2023. Spatially varying coefficients can improve parsimony and descriptive power for species distribution models. Ecography 2023:e06510. https://doi.org/10.1111/ecog.06510

U.S. Geological Survey and Canadian Wildlife Service. 2020. Strategic plan for the North American Breeding Bird Survey, 2020-30. Circular 1466. U.S. Geological Survey, Laurel, Maryland, USA. https://doi.org/10.3133/cir1466

Van Wilgenburg, S., E. Beck, B. Obermayer, T. Joyce, and B. Weddle. 2015. Biased representation of disturbance rates in the roadside sampling frame in boreal forests: implications for monitoring design. Avian Conservation and Ecology 10(2):5. https://doi.org/10.5751/ACE-00777-100205

Veech, J. A., K. L. Pardieck, and D. J. Ziolkowski Jr. 2017. How well do route survey areas represent landscapes at larger spatial extents? An analysis of land cover composition along Breeding Bird Survey routes. Condor 119:607-615. https://doi.org/10.1650/CONDOR-17-15.1

Ver Hoef, J. M., E. E. Peterson, M. B. Hooten, E. M. Hanks, and M.-J. Fortin. 2018. Spatial autoregressive models for statistical inference from ecological data. Ecological Monographs 88:36-59. https://doi.org/10.1002/ecm.1283

Wilson, S., S. L. LaDeau, A. P. Tøttrup, and P. P. Marra. 2011. Range-wide effects of breeding- and nonbreeding-season climate on the abundance of a Neotropical migrant songbird. Ecology 92:1789-1798. https://doi.org/10.1890/10-1757.1

Wilson, S., A. C. Smith, and I. Naujokaitis-Lewis. 2018. Opposing responses to drought shape spatial population dynamics of declining grassland birds. Diversity and Distributions 24:1687-1698. https://doi.org/10.1111/ddi.12811

Wright, W. J., K. M. Irvine, T. J. Rodhouse, and A. R. Litt. 2021. Spatial Gaussian processes improve multi-species occupancy models when range boundaries are uncertain and nonoverlapping. Ecology and Evolution 11:8516-8527. https://doi.org/10.1002/ece3.7629

Zhao, Q., Q. S. Latif, B. L. Nuse, D. C. Pavlacky Jr., C. L. Kilner, T. B. Ryder, and C. E. Latimer. 2024. Integrating counts from rigorous surveys and participatory science to better understand spatiotemporal variation in population processes. Methods in Ecology and Evolution 15:1380-1393. https://doi.org/10.1111/2041-210X.14368