Random wiggly smooths for 100s of factor levels

Question

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.

Answer 1

A couple of things you can do:

1. Turn on a form of multithreaded fitting: use the nthreads argument or a multithreaded BLAS - see the Details section of ?bam for details, or cluster the cluster argument (again see the Details) to estimate the chunked model in parallel.
2. Drop the s(site, bs = "re") term as random intercepts are already included in the s(year, site, bs = "fs", k = 5) term. This will save you 400+ columns in the model matrix, and these kinds of terms often crank the computational time of a GAM considerably.
3. Use the samfrac argument to select a subset of the 20,000 rows for initial estimation of the coefficients, before finishing off fitting with the full data.

I would start with 2 but spend time figuring out 1 especially if you are running your models on a modern machine with plenty of cores (the gain is less if your main machine is a laptop as you'll want to fit on all the cores but then won't be able to use your machine while it is fitting) — but do remember to turn off hyperthreading or only use the number of physical cores on your system. 3 can help if you really need it, but the other options should help you more.

When Simon says the fs terms are efficient for gamm(), I think he means compared with their use in gam(); for modest data set sizes like yours and above, I've never had any success using gamm(); I just make sure that I have a lot of available RAM and use multithreaded computation to the extent allowed by gam() or bam(), and leave the models running over lunch or in the background while I'm doing something else.