I'm working with a data set of annual bird migration counts, and I'd like to fit phenological models to determine how timing of migration has changed over the last few decades. Weather permitting, observers count the number of migrants passing each day in the fall (August-November; start and end dates vary by site). The total number of minutes spent counting birds each day (i.e., observer effort) is also recorded:
# A tibble: 366,158 x 7 site year cal_date doy minutes species count <chr> <fct> <date> <dbl> <dbl> <chr> <dbl> 1 siteA 1985 1985-09-06 249 510 spp1 0 2 siteA 1985 1985-09-07 250 480 spp1 0 3 siteA 1985 1985-09-08 251 500 spp1 0 4 siteA 1985 1985-09-09 252 570 spp1 0 5 siteA 1985 1985-09-10 253 480 spp1 0 6 siteA 1985 1985-09-11 254 180 spp1 0 7 siteA 1985 1985-09-12 255 540 spp1 0 8 siteA 1985 1985-09-13 256 495 spp1 0 9 siteA 1985 1985-09-14 257 480 spp1 0 10 siteA 1985 1985-09-15 258 465 spp1 0
From this analysis, I'm hoping to
- predict the daily count of each species by year:
- use the modeled daily count to estimate the cumulative number of birds over the course of the season by year:
- derive percentiles from the cumulative number of birds each year to represent early (0.1), peak (0.5), and late (0.9) stages of migration (this follows up on a recent question I posted here - calculating percentiles (quantiles) from GAM predictions in R). In a second analysis, I will model trends in the timing of each of these percentiles to understand whether birds are migrating earlier or later over time.
I've started by analyzing daily counts by year for a single species at a single site. This subset of the data (and the full data set) seems well suited to using a hierarchical GAM approach, and I've fit some of the model structures outlined in the Pedersen et al. 2019 HGAM paper. For example, I fit model 'GS' (global plus group-level smoothers) with smoothers for Julian date (or day-of-year, doy) and year:
modGS <- bam( # used bam() to speed things up count ~ s(doy, m = 2) + s(doy, year, bs = "fs", m = 2) + offset(log(minutes)), # account for daily observer effort data = df, method = "fREML", family = nb() )
I'd like to employ similar models to the full data set I'll actually be working with, but the structure is a bit more complex than the example above. Migrants are counted at a network of 8 sites across the western US (spanning a latitudinal gradient from northern Washington to southern Texas) and >25 species are recorded, although I'm only interested in ~12 species for this analysis (these dozen or so species are recorded at all sites). Monitoring duration also differs by site, with some locations beginning counts in the 1980s, 1990s, and early 2000s.
My question is how could I scale-up the HGAMs to do what I outlined above (provide predicted daily counts per year), but across the full multi-species, multi-site data set? Would it be a potential improvement to share information between species and sites? And how might that be specified in gam()? I expect that the effect of doy will be both species-specific and site-specific as some species migrate earlier than others, and the same species may pass through a northern site earlier than a southern site.