I have a GAM of presence of a bird species from a large dataset. There are over 400 sites in my data, where each site corresponds to a unique observer. Site-level variation is the biggest source of variation in the data. There are 30+ years of data, and 20,000 rows.
My goal is to capture the overall trend in presence through time, and the yearly smooths by site should be like random departures from that overall trend.
Here's the model I have:
m_occ <- bam(present ~ s(site, bs = "re") + s(year, site, bs = "fs", k = 5) + te(expertise, llcent) + fye + offset(log.n.periods) + s(tmin) + ti(visit) + ti(year, k = 33) + ti(visit, year), data = df_train, family = "binomial", gamma = 1.4, weights = neg_pos_weights, select = T, discrete = T, control = gam.control(trace = T))
Clearly there's a lot going on, but I want to focus on this part:
s(site, bs = "re") + s(year, site, bs = "fs", k = 5)
Here I have a random intercept for each site and a random wiggly smooth by year for each site. This is great, it does what I want. But, it takes a long time to fit, and I need to do 50 bootstraps across 60 species across 13 study areas.
Is there another approach that would be similar that might speed things up? A linear trend fits efficiently (several seconds vs. 10+ minutes), but I just want to add a little extra wiggle.
s(site, bs = "re") + s(site, by = year, bs = "re")
GAMs are so powerful, I'm optimistic there must be a solution!! In Simon Wood's GAM book this is covered in 7.7.4, but he doesn't talk about situations where you have lots of data. Also, in the factor.smooth vignette/documentation, Simon Wood also mentions that one formulation is especially fast in
gamm(), so I tried that, but it sat for 10+ minutes on the first iteration, so I think gamm() can't handle lots of data like bam() can.