Specifying hierarchical GAM for ecological count data---annual bird migration counts

Question

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

1. predict the daily count of each species by year:

2. use the modeled daily count to estimate the cumulative number of birds over the course of the season by year:

and

3. 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.

Answer 1

The way to extend these models to higher order terms is to use tensor product smooths. You can get exactly the same smooth as a bs = 'fs' term by using t2(x, f, bs = c('cr','re'), full = TRUE) so you could write your model as:

count ~ s(doy, m = 2) +  t2(doy, year, bs = c('cr', 're'), full = TRUE) +
  offset(log(minutes))

This allows us to extend these effects to higher-order effects because so long as you have the data, memory, and compute power to do it, tensor products can take as many marginal smooths as you want. So say we wanted to extend the model to have the DoY functional effect vary by Year and by Species, we could use:

count ~ s(doy, m = 2) +
  t2(doy, year, species, bs = c('cr', 're', 're'), full = TRUE) +
  offset(log(minutes))

where we're just bolting on another level of random effects, and are assuming that year and species are coded as factors.

You can also use the te() or, importantly in the case of your question here, the ti() smooths, but the parameteristation won't be exactly the same as the 'fs' smooth.

Why ti()? Well, once you start having multiple smooth effects of DoY occurring in the model, you can run into problems because as clever as mgcv is, it can't always remove all the redundant terms from the bases/model matrix that arise from having a covariate pop up in many smooth terms.

It can also help to use ti() smooths because you can partition up the problem and use the summary() output for the model to determine which levels of the hierarchy you see effects.

## first order effects
s(doy) +
s(species, bs = 're') +
s(year, bs = 're') +
s(site, bs = 're') +
## second & third order random effects
s(year, site, bs = 're') +
s(year, species, bs = 're') +
s(site, species, bs = 're') +
s(year, site, species, bs = 're') +
## second order functional effects
s(doy, year, bs = 'fs') +
s(doy, site, bs = 'fs') +
s(doy, species, bs = 'fs') +
## higher order functional effects
t2(doy, species, year, bs = c('cr', 're', 're'), full = TRUE) +
t2(doy, species, site, bs = c('cr', 're', 're'), full = TRUE) +
t2(doy, site, year, bs = c('cr', 're', 're'), full = TRUE) +
t2(doy, species, site, year, bs = c('cr', 're', 're', 're'), full = TRUE)

If year was treated as a continuous variable then you're likely to want to change this to something like:

s(doy) +
s(year) +
s(species, bs = 're') +
s(site, bs = 're') +
## second & third order random effects
s(site, species, bs = 're') +
ti(year, site, species, bs = c('cr', 're', 're')) +
## second order functional effects
s(doy, site, bs = 'fs') +
s(doy, species, bs = 'fs') +
s(year, site, bs = 'fs') +
s(year, species, bs = 'fs') +
ti(doy, year) +
## higher order functional effects
ti(doy, year, species, bs = c('cr', 'cr', 're')) +
ti(doy, species, site, bs = c('cr', 're', 're')) +
ti(doy, year, site, bs = c('cr', 'cr', 're')) +
ti(doy, year, species, site, bs = c('cr', 'cr', 're', 're'))

which is where ti() can come in useful as we've got a tonne of terms that now include smooth functions of doy and or year.

In both these models we see some duplication; the fs basis in

s(doy, year, bs = 'fs')

contains random intercepts for year but you'll notice that I also include the first order ranef for year via s(year, bs = 're'). Again, I think this ins fine and it's how we did it in Pedersen et al (2019). The random intercepts in the factor-smooths are fully penalised so the first order year ranef should account for the overall between year variation, and the random intercept aspects of the factor-smooth or higher-order functional effects should not be very large at all.

I've assumed that you want site as a random effect and not a smooth effect; the latter may be useful when you more sites along a gradient where you expect smooth effects along the spatial gradient/transect.

If you had more sites over a range of lat and long, you could also model the spatial effect via a spatial smooth s(long, lat, bs = 'ds') for example.

References

Pedersen, E.J., Miller, D.L., Simpson, G.L., Ross, N., 2019. Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7, e6876. https://doi.org/10.7717/peerj.6876