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:

enter image description here

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

enter image description here


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

  • $\begingroup$ HI Jason, that's a pretty open-ended question as you could allow all terms to covary with all the others implying several hierarchical effects. I'm a little swamped with work today but I'll write something tomorrow, which addresses some of your points and which should (hopefully) indicate how to proceed building up to a full model with all the terms in it. $\endgroup$ Commented Feb 26, 2020 at 16:23
  • $\begingroup$ Hi Gavin, sorry for the ambiguity. I would imagine that running separate analyses by species and site (as I had done above) is pretty inefficient, and that there's some information that is being thrown away that could be shared across sites or species. Maybe one approach would be analyzing each species separately, and having a hierarchical structure in the model to borrow information across sites about the relationship between day-of-year and year? Especially if one or more sites doesn't have many counts in a season? Not sure if that make sense, curious to hear what you recommend. $\endgroup$
    – Jason
    Commented Feb 28, 2020 at 17:43
  • $\begingroup$ It's not ambiguous; it's more that these models get complicated fast and there are many permutations for how to represent all the effects. Also, when I wrote my initial comment I was thinking you had a more complex spatial arrangement than you do. Even with the simpler "space-as-site-ranef" the model complexity has blown up a bit as per my answer below. Sorry it took me a few days longer to get to this; dealing with a lot of work-related stuff recently. $\endgroup$ Commented Mar 5, 2020 at 22:53

1 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) +

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) +

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.


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.