Overview

We designed a Bayesian hierarchical model for estimating status and trends from the North American BBS that uses a GAM smooth to estimate the medium- and long-term temporal components of a species population trajectory (i.e. changes occurring over time periods ranging from 3 to 53 yr). In the model, the parameters of the GAM smooths are treated as random effects within the geographic strata (the spatial units of the predictions, intersections of Bird Conservation Regions, and province/state/territory boundaries), so that information is shared on the shape of the population trajectory across the species' range. In comparison to the non-Bayesian HGAMs in the work of Pedersen et al. 2019, our model is most similar to the “GS” model, which has a global smooth in addition to group-level smooths with a similar degree of flexibility. We applied 2 versions of the GAM: one in which the GAM smooth was the only component modeling changes in abundance over time (GAM) and another in which random year effects were also estimated to allow for single-year departures from the GAM smooth (GAMYE, which is conceptually similar to the model described in the work of Knape 2016).

For a selection of species, we compared estimates and predictive accuracy of our 2 models using the GAM smooth, against 2 alternative models that have been used to analyze the BBS data. We chose the main comparison species (Barn Swallow, Hirundo rustica) because of the striking differences between trajectories from the SLOPE model and a number of nonlinear models (Smith et al. 2015, Sauer and Link 2017). We added a selection of other species to represent a range of anticipated patterns of population change, including species with known change points in their population trajectories (Chimney Swift, Chaetura pelagica; Smith et al. 2015), and species with relatively more data and known large and long-term trends (Wood Thrush, Hylocichla mustelina and Ruby-throated Hummingbird, Archilochus colubris) and species with relatively fewer data and long-term changes (Canada Warbler, Cardellina canadensis, Cooper's Hawk, Accipiter cooperii, and Chestnut-collared Longspur, Calcarius ornatus). Finally, we also added a few species with strong annual fluctuations and/or abrupt step-changes in abundance (Pine Siskin, Spinus pinus and Carolina Wren, Thryothorus ludovicianus).

The BBS data are collected along roadside survey routes that include 50 stops at which a 3-min point count is conducted, once annually, during the peak of the breeding season (Robbins et al. 1986, Hudson et al. 2017, Sauer et al. 2017). All of the models here use the count of individual birds observed on each BBS route (summed across all 50 stops) in a given year by a particular observer. The 4 statistical models differed only in the parameters used to model changes in species relative abundance over time. We used 15-fold cross-validation (Burman 1989) to estimate the observation-level, out-of-sample predictive accuracy of all 4 models (Vehtari et al. 2017, Link et al. 2020). We compared the overall predictive accuracy among the models, and we explored the spatial and temporal variation in predictive accuracy in depth.

We compared 4 alternative BBS models, all of which have the same basic structure:

The models treat the observed BBS counts as overdispersed Poisson random variables, with mean λ_s,j,t (i.e. geographic stratum s, observer and route combination j, and year t). The means are log-linear functions of stratum-specific intercepts (θ_s, estimated as fixed effects and with the same priors following Smith et al. 2014), observer-route effects (ω_j, estimated as random effects and with the same priors following Sauer and Link 2011), first-year startup effects for an observer(η, estimated as fixed effects and with the same priors following Sauer and Link 2011), a count-level random effect to model overdispersion (η_s,j,t, estimated using heavy-tailed, t-distribution and with the same priors following Link et al. 2020), and a temporal component estimated using a function of year, which varies across the 4 models (Δ_s (t)). The models here only varied in their temporal components (Δ_s (t)).

Bayesian Hierarchical GAMs

Generalized additive model. The main temporal component Δ_s (t) in the GAM was modeled with a semi-parametric smooth, estimated following Crainiceanu et al. (2005) as

where K is the number of knots, χ_t,k is the year t and k th entry in the design matrix X(defined below), and β_s,k is the K-length vector of parameters that control the shapeof the trajectory in stratum s. Each β_s,k is estimated as a random effect, centered on a mean across all strata (a hyperparameter B_k)

and

where the variance σ²_B acts as the complexity penalty, shrinking the complexity and the overall change of the mean trajectory toward a flat line). It would be possible to add an additional slope parameter, as was done in the work of Crainiceanu et al. 2005, but we have found that the BBS data for most species are insufficient to allow for the separate estimation of the linear component to population change and the additive smooth. In addition, we see little benefit to including a linear component because the assumptions required to include a constant linear slope for a 53-yr time series are unlikely to be met for any continental-scale population. In combination, these variance parameters (σ²_β, σ²_B) control the complexity penalty of the species trajectories and the variation in pattern and complexity among strata and were given the following priors, following advice in the work of Crainiceanu et al (2005):

These prior parameters were chosen to ensure that the priors are sufficiently vague that they are overwhelmed by the data, particularly for σ²_B that controls the shape of the survey-wide trajectory (Crainiceanu et al 2005). We have so far had good results across a wide range of species using these priors, and in tests of alternative priors, there is no effect on posterior estimates ( Supplementary Material Figure S9). For example, estimates of B_K and σ_B for Chestnut-collared Longspur (a relatively data-poor species) are unchanged even if using a much more restrictive prior on σ_B that places 99% of the prior density for σ_B below 1.2 . However, these variance priors are an area of ongoing research, aimed at improving the efficiency of the Markov Chain Monte Carlo (MCMC) sampling.

The design matrix for the smoothing function (X) has a row for each year and a column for each of K knots. The GAM smooth represented a third-degree polynomial spline:χ_t,k = |t′ – t′_k|³, and was calculated in R, following the work of Crainiceanu et al (2005). We centered and re-scaled the year values to improve convergence, so that , where midyear is the middle year of the time series and T is the number of years in the time series. Here, we have used 13 knots (K = 13), across the 53-yr time series of the BBS (1966–2018), which results in approximately one knot for every 4 yr in the time series. With this number of knots, we have found that the 53-yr trajectories are sufficiently flexible to capture all but the shortest-term variation (i.e. variation on the scale of 3–53 yr, but not annual fluctuations). Models with more knots are possible but in the case of a penalized smooth, the overall patterns in the trajectory will be similar, as long as a sufficient number of knots is allowed (Wood 2017). The number of knots could be customized in a species-specific approach; however, because we are looking for a general model structure that can be applied similarly across the >500 species in the BBS, we have fixed the number of knots at 13. Our approach relies on the shrinkage of the smoothing parameters (B, β) to ensure that the trajectories are only as complex as the data support, and the limited number of knots constrains the complexity of the additive function (Fewster et al. 2000, Wood 2017).

GAMYE. The GAMYE was identical to the GAM, with the addition of random year effects (γ_t,s) estimated independently among strata, following the works of Sauer and Link (2011) and Smith et al. (2015), as

where σ²_γ,s are stratum-specific variances. Thus, the temporal component for the GAMYE is given by

The GAMYE trajectories are therefore an additive combination of the smooth and random annual fluctuations. The smooth components of the trajectory in the GAMYE are generally similar to those from the GAM, but tend to be slightly less variable (i.e. less complex) because the year effects components can account for single-year deviations from the longer-term patterns of population change. The full trajectories from the GAMYE (smooth plus the year effects) generally follow the same overall pattern as the GAM estimates and include abrupt single-year changes in abundance, which increases the capacity to model step-changes in abundance.

Alternative Models

For a selection of species, we compared the predictions and predictive accuracy of the 2 GAMs against 2 alternative models previously used for the BBS.

SLOPE. The SLOPE model includes a slope parameter and random year effects to model species trajectories. It is a linear year-effects model currently used by both the CWS (Smith et al. 2014) and the USGS (Sauer et al. 2017) as an omnibus model to supply status and trend estimates from the BBS (essentially the same as model SH, the SLOPE model with heavy-tailed error in the work of Link et al. 2017). The temporal component in the SLOPE model is

DIFFERENCE. The first difference model (DIFFERENCE) is based on a model described in the work of Link and Sauer (2016) and models the temporal component as

The DIFFERENCE model includes year effects that follow a random walk prior from the first year of the time series, by modeling the first-order differences between years as random effects with mean zero and an estimated variance.

All analyses in this article were conducted in R (R Core Team 2019), using JAGS to implement the Bayesian analyses (Plummer 2003), and an R-package bbsBayes (Edwards and Smith 2020) to access the BBS data and run all of the models used here. We used the same number of burn-in iterations (10,000), thinning rate (1/10), chains (3), and number of saved samples from the posterior (3,000) to estimate trends and trajectories for all models. We examined trace plots and the Rhat statistic to assess convergence. The graphs relied heavily on the package ggplot2 (Wickham 2016). BUGS-language descriptions of the GAM and GAMYE, as well as all the code and data used to produce the analyses in this study, are archived online (see Data availability in Acknowledgements).

Cross-Validation

We used a temporally and spatially stratified v-fold cross-validation (Burman 1989, often termed “k-fold,” but here we use Berman's original “v-fold” to distinguish it from “k” which is often used to describe the number of knots in a GAM) with v = 15, where we held-out random sets of counts, stratified across all years and strata so that each of the v-folds included some observations from almost every combination of strata and years. We did this by randomly allocating each count within a given stratum and year to one of the 15 folds. We chose this approach over a leave-one-out (loo) cross-validation approach using a random subset of counts (Link et al. 2017) because we wanted to assess the predictive success across all counts in the dataset, explore the temporal and spatial patterns in predictive success. Also, although full loo cross-validation minimizes bias and variance of the estimate of predictive accuracy (Zhang and Yang 2015), a full loo is not practical for computational reasons and cross-validation with k >10 is a reasonable approximation of loo (Kohavi 1995, Vehtari et al. 2017). We could also have chosen to conduct structured cross-validation (Roberts et al. 2017), but cross-validation in a Bayesian context has particularly large computational requirements; there are multiple levels of dependencies in the BBS data (dependences in time, space, and across observers); and models being compared vary in the way they treat some of those dependencies (models that share information differently in space and/or time). Therefore, we chose a relatively simple non-structured approach where the folds are balanced in time and space, and for a given species were identical across all models compared. We followed a similar procedure to that outlined in the work of Link et al. (2017) to implement the cross-validation in a parallel computing environment, using the R-package for each (Wallig and Weston 2020). We used the end-values from the model-run using the full dataset as initial values in each of the 15 cross-validation runs, ran a short burn-in of 1,000 samples, then used a draw of 3,000 samples of the posterior with a thinning rate of 1/10 spread across 3 chains. We did not calculate the widely applicable information criterion (WAIC), because previous work has shown that WAIC does not approximate loo well for the BBS data (Link et al. 2017, 2020).

We used the estimated log predictive density (elpd_i,M) to compare the observation-level, out-of-sample predictive success of all 4 models (Vehtari et al. 2017, Link et al. 2020). For a given model M, elpd is the estimated log posterior density for each observation i, for the model M fit to all data except those in the set v that includes i (Y_{—v, i—v}). That is,

Larger values of elpd indicate better predictive success, that is a higher probability of the observed data given the model M, the estimated parameters, the vector of covariates for observation i, such as the year, observer-route, etc. (X_i), and all of the data used to fit the model (Y_{—v, i—v}).

We have not summed elpd values to generate Bayesian predictive information criterion (BPIC) values; rather, we have compared model-based estimates of mean difference inelpd between pairs of models. We modeled the elpd values so that we could account for the imbalances in the BBS data in time and space (i.e. the variation in number of counts among strata and years). The raw sum of the elpd values would give greater weight to the regions with more data and to the recent years in the time series, which have more counts. Therefore, expanding on the approach in the work of Link et al. 2020 that used a z-score to estimate the significance of the difference in fit between 2 models, we used a hierarchical model to estimate the mean of the differences in predictive fit (). We first calculated the difference in the elpd of each observed count (Y_i) under models 1 and 2, as = log (f₁ (Y_i|Y_{—v, i—v}, X_i)) – log (f₂ (Y_i|Y_{–v, i∈v}, X_i)), so that positive values of indicate more support for model 1. We then analyzed these values using an additional Bayesian hierarchical model, with random effects for year and strata to account for the variation in sampling effort in time and space. These random effects account for the imbalances in the BBS data among years and regions, and the inherent uncertainty associated with any cross-validation statistic (Link et al. 2017, Vehtari et al. 2017). This model treated the elpd differences for a count from a given year t and stratum s () as having a t-distribution with an estimated variance () and degrees of freedom (ν). That is,

From the model, ϕ was our estimate of the overall comparison of the mean difference in predictive fit for Model 1 – Model 2 ( = ϕ), ϕ + Ψ_s was the estimate of the mean difference in stratum s, and ϕ + Ψ_t was the estimated difference in year t. The year and stratum effects (Ψ_s + Ψ_t) were estimated as random effects with a mean of zero and estimated variances given uninformative inverse gamma priors. We used this t-distribution as a robust estimation approach, instead of the z-score approach used by Link et al. (2020) because of the heavy tails in the distribution of the values ( Supplementary Material Figure S7). Given these heavy tails, a statistical analysis assuming a normal distribution in the differences would give an inappropriately large weight to a few counts where the prediction differences were large in magnitude (Gelman et al. 2014). In essence, our model is simply a “robust” version of the z-score approach (Lange et al. 1989) with the added hierarchical parameters to account for the spatial and temporal imbalance in the BBS data.

Trends and Population Trajectories

For all models, we used the same definition of trend following Sauer and Link (2011) and Smith et al. (2015); that is, an interval-specific geometric mean of proportional changes in population size, expressed as a percentage. Thus, the trend estimate for the interval from year a (t_a) through year b (t_b) is given by

where N_. represents the annual index of abundance in a given year (see below). Because this estimate of trend only considers the annual abundance estimates in the years at either end of the trend period, we refer to this estimate as an end-point trend. For the GAMYE model, we decomposed the trajectory (i.e. the series of annual indices of abundance) into long- and medium-term components represented by the GAM smooth and annual fluctuations represented by the random year effects. This decomposition allowed us to estimate 2 kinds of trend estimates: R_a:b that include all aspects of the trajectory, and R′_a:_b that removes the annual fluctuations, including only the GAM smooth components.

Population trajectories are the collection of annual indices of relative abundance across the time series. These indices approximate the mean count on an average BBS route, conducted by an average observer, in a given stratum and year. For all the models here, we calculated these annual indices for each year t and stratum s following the work of Smith et al. (2019) as

where each N_s,t is an exponentiated sum of the relevant components of the model (A_s,t), observer-route effects (ω_j), and count-level extra-Poisson variance (0.5 × σ²_ε), averaged over count-scale predictions across all of the n_Js observer-routes j in the set of observer-route combinations in stratum s (J_s), and then multiplied by a correction factor for the proportion of routes in the stratum on which the species has been observed (z_s, i.e. the proportion of routes on which the species has been observed; on all other routes species abundance is assumed to equal zero and they are excluded from the model, see Sauer and Link 2011). This is slightly different from the approach described in the works of Sauer and Link (2011) and Smith et al. (2015) and an area of ongoing research. We have found that this different annual index calculation ensures that the annual indices are scaled more similarly to the observed mean counts, which can affect the relative weight of different strata in trends estimated for broader regions (e.g., continental and national trends), but it has no effect on the trends estimated within each stratum and no effect on the cross-validation results presented here. For a discussion on the differences between these 2 ways of calculating annual indices, refer to  Supplementary Material Appendix A.

For the GAMYE model, we calculated 2 versions of the species trajectory (N_s): one that included the annual variation in the trajectory,

and a second that excluded the annual variations, including only the smoothing components of the GAM to estimate the time series,

We calculated population trajectories and trends from the GAMYE model using both sets of annual indices (N_s,t and Ng_s,t). When comparing predictions against the other models, we used the N_s,t values to plot and compare the population trajectories (i.e. including the year effects), and the Ng_s,t values to calculate the trends (i.e. removing the year-effect fluctuations).

Model Predictions

Population trajectories from the GAM, GAMYE, and DIFFERENCE are similar. All 3 of these models suggest that Barn Swallow populations increased from the start of the survey to approximately the early 1980s, compared to the SLOPE model predictions that show a relatively steady rate of decline (Figure 1). The trajectories for all species from both GAMs and the DIFFERENCE model were less linear overall than the SLOPE model trajectories and tended to better track nonlinear patterns, particularly in the early years of the survey and often in more recent years as well (Figure 1,  Supplementary Material Figures S1 and  S6). GAM and GAMYE trajectories vary a great deal among the strata, particularly in the magnitude and direction of the long-term change (Figure 2 for Barn Swallow). However, there are also many similarities among the strata, in the nonlinear patterns that are evident in the continental mean trajectory (e.g., the downward inflection in the early 1980s in Figure 2 and  Supplementary Material Figure S2). Figure 3 shows the estimated trajectories for Barn Swallow in the 6 strata that make up BCR 23 from the GAMYE, DIFFERENCE, and SLOPE models. The GAMYE estimates suggest that the species' populations increased in the early portion of the time series in all of the strata, and this is a pattern shared with the continental mean trajectory for the species (Figure 2). In contrast, the estimates from the SLOPE model only show an increase in the stratum with the most data (i.e. the most stacked gray dots along the x-axis indicating the number of BBS routes contributing data in each year, US-WI-23), the DIFFERENCE model shows more of the early increase in many strata, except those with the fewest data. In the other strata with fewer data, the SLOPE trajectories are strongly linear and the DIFFERENCE trajectories are particularly flat in the early years with particularly few data. The cross-validation results suggest that for Barn Swallow, the GAMYE is preferred over the SLOPE model, and generally preferred (some overlap with 0) to the DIFFERENCE model (Figure 4), particularly in the early years of the survey (pre-1975,  Supplementary Material Figure S6). Finally, the general benefits of sharing information among strata on the shape of the population trajectory are evident for the GAM, GAMYE, and the SLOPE models in Figure 5, where there is no relationship between the sample size and the absolute value of the long-term trend for Cooper's Hawk (more below).

FIGURE 1.

Survey-wide population trajectories for Barn Swallow (Hirundo rustica), Wood Thrush (Hylocichla mustelina), Cooper's Hawk (Accipiter cooperii), Carolina Wren (Thryothorus ludovicianus), Ruby-throated Hummingbird (Archilochus colubris), and Chimney Swift (Chaetura pelagica), estimated from the BBS using 2 models described here that include a GAM smoothing function to model change over time (GAM and GAMYE) the standard regression-based model used for BBS status and trend assessments since 2011 (SLOPE), and a first difference time series model (DIFFERENCE). The stacked dots along the x-axis indicate the approximate number of BBS counts used in the model in each year; each dot represents 50 counts.

FIGURE 2.

Variation among the spatial strata in the random-effect smooth components of the GAMYE model applied to Barn Swallow data from the BBS. Gray lines show the strata-level random-effect smooths, and the black lines show the survey-wide mean trajectory.

For most species here, the GAMs or the DIFFERENCE model generally were preferred over the SLOPE model (Figure 4,  Supplementary Material Figure S3). For the 2 species with population trajectories that are known to include strong year effects (Carolina Wren and Pine Siskin), the GAM model that does not include year effects performed poorly (Figure 4). For Carolina Wren, the DIFFERENCE model was preferred clearly over the GAMYE (Figure 4), and yet the predicted trajectories from the 2 models are similar (Figure 1). In contrast, for Pine Siskin, the DIFFERENCE and GAMYE were similar in their predictive accuracy (Figure 4) and yet the predicted trajectories are noticeably different in the first 10 yr of the survey ( Supplementary Material Figure S1). For Cooper's Hawk, the GAMYE model was generally preferred over the DIFFERENCE model, although there was some overlap with zero (Figure 4), but in this case, the predicted trajectories are different. The DIFFERENCE trajectory for Cooper's Hawk suggests much less change in the species' population over time than the GAM or GAMYE (Figure 1).

Cooper's Hawk provides an example of a species with sparse data, for which the sharing of information in space may be particularly relevant. In a single stratum, the model has relatively a few data with which to estimate changes in populations through time. For example, the mean counts for the species indicate that on average one bird was observed for every 40 BBS routes run in the 1970s, and since the species population has increased, it still requires more than 10 routes to observe a single bird. For this species, the models that share information among strata on population change (GAM, GAMYE, and SLOPE) suggest a greater change in populations than the DIFFERENCE. For these models, where the stratum-level population change parameters are able to share information across the species' range, the absolute change in the population does not depend on the sample size in the region. In addition, for each of these models, there is still large variability in the trends estimated for data-sparse regions, demonstrating that while the estimates benefit from the sharing of information among strata, the local trends are still influenced by the local data. In contrast, there is a strong relationship between the magnitude of the trend and the number of routes contributing data to the analysis for the DIFFERENCE model (Figure 5,  Supplementary Material Figure S4). In strata with fewer than 10 routes contributing data, the DIFFERENCE trends are almost all close to zero. In these relatively data-sparse strata, the DIFFERENCE model has little information available to estimate population change, and so the prior is more relevant and the population changes are shrunk toward zero. In contrast, the other models can integrate data from the local stratum with information on changes in the species' population across the rest of its range.

The decomposed trajectories from the GAMYE allow us to calculate trends from the smooth but also plot trajectories that show the annual fluctuations ( Supplementary Material Figure S5). For example, the smooth trajectory for the Carolina Wren captures the general patterns of increases and decreases well, while the full trajectory also shows the sharp population crash associated with the extreme winter in 1976 (Figure 6). Calculating trends from the smooth component generates short-term estimates that vary less from year to year for species with relatively strong annual fluctuations (Figure 7). For example, Figure 8 shows the series of short-term (10-yr) trend estimates for Wood Thrush in Canada, from the smooth component of the GAMYE, the GAMYE including the year effects, the DIFFERENCE model, and the SLOPE model used since 2011. In this example, the 10-yr trend estimate from the GAMYE with the year effects and the SLOPE model both cross the IUCN trend threshold criterion for Threatened (IUCN 2019) at least once in the last 12 yr, including 2011, when the species' status was assessed in Canada (COSEWIC 2012). In contrast, a trend calculated from the decomposed GAMYE model using only the smooth component (GAMYE—Smooth Only in Figure 8) fluctuates much less between years.

Cross-Validation Varies in Time and Space

The preferred model from the pairwise predictive fit comparisons varied in time and space (Figures 4, 9, and 10 and  Supplementary Material Figure S6). The contrast between GAMYE and DIFFERENCE for Barn Swallow provides a useful example: depending on the year or the region of the continent, either the GAMYE or the DIFFERENCE model was the preferred model, but overall, and in almost all regions and years, the 95% CI of the mean difference in fit between GAMYE and DIFFERENCE overlapped 0 (Figures 4, 9, and 10). For Barn Swallow, the GAMYE model has generally higher predictive fit during the first 5 yr of the time series, but then the DIFFERENCE model is preferred between approximately 1975 and 1983. The geographic variation in predictive fit is similarly complex. In the eastern parts of the Barn Swallow's range, the GAMYE model generally out-performed the DIFFERENCE model, whereas the reverse is generally true in the remainder of the species' range (Figure 10). Although the mapped colors only represent the point-estimates, they suggest an interesting spatial pattern in the predictive fit of these 2 models for this species. Many of the species considered here show similarly complex temporal and spatial patterns in the preferred model based on predictive fit ( Supplementary Material Figure S6).

FIGURE 3.

Stratum-level predictions for Barn Swallow population trajectories in BCR 23 from GAMYE against the predictions from the SLOPE and DIFFERENCE model. The similarity of the overall patterns in the GAMYE as compared to the SLOPE estimates demonstrates the inferential benefits of the sharing of information among regions on the nonlinear shape of the trajectory. In most strata, the similar patterns of observed mean counts and the GAMYE trajectories suggest a steep increase in Barn Swallows across all of BCR 23 during the first 10 yr of the survey. The GAMYE estimates show this steep increase in almost all of the strata, whereas the SLOPE predictions only show this pattern in the most data-rich stratum, US-WI-23. The DIFFERENCE trajectories model the nonlinear shapes well in all but the most data-sparse region (US-IL-23) and years (<1972). The facet strip labels indicate the country and state-level division of BCR 23 that makes up each stratum. The first 2 letters indicate all strata are within the United States, and the second 2 letters indicate the state. The stacked dots along the x-axis indicate the number of BBS counts in each year and stratum; each dot represents one count.

Predictive Accuracy

Overall, the cross-validation comparisons generally support the GAMYE, GAM, or DIFFERENCE model over the SLOPE model for the species considered here, in agreement with Link et al. (2020). For Barn Swallow, the overall difference in predictive fit, and particularly the increasing predictive error of the SLOPE model in the earliest years, strongly suggests that in the period between the start of the BBS (1966) and ∼1983 (Smith et al. 2015), Barn Swallow populations increased. All models agree, however, that since the mid-1980s populations have decreased.

Using all data in our cross-validations allowed us to explore the spatial and temporal variation in fit and to compare the fit across all data used in the model. We have not reported absolute values of predictive fit because estimates of fit from a random selection of BBS counts, or simple summaries of predictive fit from the full dataset, are biased by the strong spatial and temporal dependencies in the BBS data (Roberts et al. 2017). However, because our folds were identical across models, and we modeled the differences in fit with an additional hierarchical model that accounted for repeated measures among strata and years, we are reasonably confident that relative-fit assessments are unbiased within a species and among models. Alternative approaches, such as blocked cross-validation (Roberts et al. 2017) to assess the predictive success of models in time and space, and targeted cross-validation (Link et al. 2017) to explore the predictive success in relation to particular inferences (e.g., predictive accuracy in the end-point years used for short- and long-term trend assessments) are an area of ongoing research.

The overall predictive fit assessments provided some guidance on model selection for the species here, but not in all cases. The SLOPE model compared poorly against most other models in the overall assessment, similar to Link et al. 2020. However, among the other 3 models, many of the overall comparisons failed to clearly support one model, even in cases where the predicted population trajectories suggested different patterns of population change (e.g., Cooper's Hawk). For a given species, the best model varied among years and strata. These temporal and spatial patterns in predictive fit complicate the selection among models, given the varied uses of the BBS status and trend estimates (Rosenberg et al. 2017).

In general, estimates of predictive accuracy are one aspect of a thoughtful model building and assessment process, but are insufficient on their own (Gelman et al. 2013 p. 180, Burnham and Anderson 2002 p. 16). This is particularly true if there is little or no clear difference in overall predictive accuracy, but important differences in model predictions. For example, the overall cross-validation results do not clearly distinguish among the SLOPE, DIFFERENCE, and GAMYE for Cooper's Hawk, and yet predictions are different between the DIFFERENCE model and the others (Figures 1, 4, and 5). Interestingly, the cross-validation approach in the work of Link et al. 2020 suggested that the DIFFERENCE model was preferred over the SLOPE for Cooper's Hawk, but we did not find that here ( Supplemental Material Figure S3). The important differences in trend estimates and the equivocal cross-validation results suggest that further research is needed into the criteria for, and consequences of, model selection for BBS status and trend estimates. Model selection is also complicated when overall predictive accuracy appears to clearly support one model and yet the important parameters (trends and trajectories) are not noticeably different. For example, the overall cross-validation results for Carolina Wren suggest that the DIFFERENCE model is preferred over the GAMYE, and yet the trajectories are almost identical (Figures 1 and 4). Predictive accuracy is also complicated when robust predictions are required for years or regions with relatively few data against which predictions can be assessed (e.g., the earlier years of the BBS, or data-sparse strata that still have an important influence on the range-wide trend). Model selection is complicated, and predictive accuracy would never be the only criterion used to select a model for the BBS analyses. Limits to computational capacity and a desire to avoid a data-dredging all-possible-models approach ensure that some thoughtful process to select the candidate models is necessary.

We agree with Link et al. (2020) that we should not select models based on a particular pattern in the results. In fact, the necessary subjective process occurs before any quantitative analyses (Burnham and Anderson 2002) and relies on “careful thinking” to balance the objectives, the model, and the data (Chatfield 1995). The careful thinking required to select a BBS model or to interpret the BBS status and trend estimates is to consider the consequences of the potential conflicts between the model structures (“constraints on the model parameters” Chatfield 1995) and the objectives of the use of the modeled estimates. For example, one of the models that shares information on population change among strata is likely preferable to the DIFFERENCE model for species with relatively sparse data in any given stratum, because the prior of the DIFFERENCE model (stable population) will be more influential when the data are sparse. This prior dependency of the results may not be identified by the lower predictive accuracy of the estimates, as the results for Cooper's Hawk demonstrate (Figure 5). Similarly, a user of estimates from the DIFFERENCE model should carefully consider the conservation-relevant consequences of the prior and model structure when assessing potential changes in the population trends of declining and relatively rare species. These species' short-term rates of decline could appear to decrease, suggesting a stabilizing population, simply due to the increasing influence of the prior, if the species observations decline to a point where it is not observed in some years. In contrast, if a user wished to explicitly compare estimates of population change among political jurisdictions or ecological units, the sharing of information among those units in the GAM-based models here might be problematic. We suggest that the GAMYE's strong cross-validation performance, its sharing of information across a species range, its decomposition of the population trajectory, and its broad utility that suits the most common uses of the BBS status and trends estimates make it a particularly useful model for the sort of omnibus analyses conducted by the CWS and other agencies.