https://doi.org/10.5061/dryad.280gb5msz

Hierarchical Breeding Bird Survey Abundance Model
To evaluate the effect of winter rainfall on population abundances of American redstarts on the breeding grounds we relied on 30 years (1990-2019) years of the long-term Breeding Bird Survey (BBS) data to coincide with our 30-year banding and isotope dataset. Data from 554 BBS routes included in this analysis were selected based on the likely breeding origins of this population (SI Appendix, Figure S4). Stable-hydrogen isotope analysis indicated that this wintering population originates from a broad latitudinal band (35 to 50 °N). To determine the longitudinal band of our population we calculated the 95% CI of breeding longitude estimated from a sample of geolocator tagged individuals (-93 to -83 °W).

The model includes estimates of log-linear population trends ().  represents the effect of week-of-year on abundance, β₃ accounts for breeding season (May – July) Total Rainfall, β₄ accounts for breeding season (May – July) average max temperature, and β₅ represents the effect of total winter rainfall (Jan – Mar) in Jamaica.  accounts for the effect of Latitude on abundance and  allows for a varying effect of winter rainfall by breeding latitude. We used uninformative normal priors with mean 0 and variance 10² for all beta parameters and diffuse uniform distributions for all variance parameters. To account for the non-independence of observations on a given route, we have included a random effect on the intercept for each route (Route_i) and observer (Observer_i,t), both with a Gaussian prior centered on 0. To compare the relative contributions of each factor on breeding bird abundance and to aid in computations, we centered and scaled all predictors so that values were in units of SD about the mean.  
 
We ran three chains of 25000 iterations, discarded the first 5000, thinned by a factor of 10 to give 6,000 samples from the posterior distribution for inference. We assessed model convergence through the parameter histogram plots and Gelman-Rubin (R-hat) convergence diagnostics. We interpreted all effects by examining posterior mean and associated 95% credible intervals (CI) for all β_i parameters.
 
Hierarchical Multi-Season Occupancy Model
We modeled changes in occupancy probability with multi-season Bayesian hierarchical occupancy model implemented in the spOccupancy package in R (60). We utilized the BBS data from the abundance model above as occurrence data (presence/absence) by treating each of the 5 segments along a BBS route (pooled 10-stop summaries, 50 stops in total per route) as a replicate and converting all counts > 1 as a detections and counts < 1 as an absence. Conventional practice in occupancy modeling is to utilize temporal replicates to estimate detection probability. BBS routes are only sampled once per year and so occupancy models have been adapted to utilize spatial replicates (stops) as a substitute (61). However, it’s important to note that this space-for-time substitution introduces potential spatial autocorrelation that may bias occupancy if not properly accounted for. We attempt to account for the correlation by including a route-level random effect on detection probability which would allow for stops within a route to be correlated. Ideally, modeling the spatial correlation directly would provide a more exact approach towards accounting for the survey structure (61). We modeled detection as a function the week of the year the route was run (centered and standardized) while accounting for random observer affects and route effects. Occupancy was modeled as a function of proportion of forest, a year trend (centered and standardized), and random route effects. The proportion of forest represents the area within 40km of a BBS route with forest cover > 10%. These data were extracted from the Global Forest Change dataset (v1.6; 62) using Google Earth Engine. We incorporated a AR(1) temporal covariance structure to formally account for its effect on occupancy dynamics. We fit the model with default non-informative on all  and weakly informative priors on the AR(1) variance and correlation parameters following Doser et al. (2022). We ran 3 chains for 10000 iterations with a burn in of 2000 iterations and thinning rate of 10 resulting in 2400 posterior samples. We assessed model convergence through the parameter histogram plots and Gelman-Rubin (R-hat) convergence diagnostics. The significance of parameter estimates was evaluated by whether or not the 95% credible interval overlapped 0.

Description of the data and file structure

Survival Data are available upon request; however, the code to run survival analysis is included here as Survival Analysis_Final.r.

BBS_Abundance_Data.rds contains all necessary data to replicate the abundance analysis, including winter and breeding season environmental covariates scaled and centered. Below is a quick description of the variables:

Count = # of Redstarts observed
Route = Route Identification Number
Observer = Observer Identification Number
Year = Year of Route scaled to number of years since 1989. Year 1989 would be coded as 0.
Week = Week of the Year with 0 indicating the first week of the year
Latitude = scaled and centered latitude (degrees north converted to units of standard deviation)
BreedingRain = Total Rainfall in mm (May - July) scaled and centered (units of standard deviation)
BreedingTemp = Average Daily Max Temperature Celsius (May - July) scaled and centered (units of standard deviation)
WinterRain = Total Winter Rain in Jamaica mm (Jan - Mar) scaled and centered (units of standard deviation)

NAs in count, observer, and week represent instances where the route survey was not conducted.

The model formulation is presented in the BBS_Abundance_Model.r script but can easily be implemented in the brms package with model description in manuscript methods.

BBS_Occupancy_Data.rds contains all the necessary data in the proper format to replicate the occupancy analysis in BBS_Occupancy_Analysis.r.

BBS Occupancy data are presented in list format. The first list item are the observed counts in an array. Rows represent the different routes, columns represent the different years, and matrix slices represent 10-stop groups (e.g., Count10, Count20, Count30, Count40, Count50). NA represent instances when the route was not surveyed.

The next dataset within the list is the spatial coordinates of each route. This is given in latitude (degree North) and longitude format (degree East).

The next dataset within the list includes all of the detection covariates. This includes Observer, site, and week of the year. Observer is stored as a matrix organized by route and year. The site includes a unit identifier for the route, and the week is organized like an observer but with a week of the year denoted as above (weeks since Jan 1) scaled and centered (units of standard deviation).

The final dataset within the occupancy data list includes all of the occupancy covariates. These include year, proportion of forest, and site. The year is organized by route and year, the site is organized as it was in the detection covariates, and the proportion of forest is organized at the route level and represents the total area covered with forest (fraction from 0-1) within 40km of the route center point.

Sharing/Access information

Isotope data are available upon request. Geolocator data are available in Dossman et al. (2022) https://doi.org/10.1002/ecy.3938.