# Random wiggly smooths for 100s of factor levels

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.

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.

• This is super helpful thank you! Should I use a factor smooth or instead something like s(year, site, bs = c("cr", "re"))? I'm not sure of the difference (big picture), but found the latter in one of your posts from a few years back. Commented Aug 30, 2023 at 13:36
• Also I read a bit about the sum-to-zero smoother, and I think it said that makes the main smooth s(year) do more work and that the factor smooths would be random departures from that. Since my primary objective is estimating s(year) (while accounting for all this other variation), is the way I have the year vs year-x-site smooths ok? Commented Aug 30, 2023 at 14:00
• Re s(year, site, bs = c("cr", "re")), that simply won't work. If you did te() instead it would be OK, better with t2(..., full = TRUE), but then the latter is just the fs smooth, so use the fs smooth. Re the second question, the ti(year, k = 33) term will take care of the average effect. If you tried to use the sz basis for factor by terms, your model will take longer to fit as you will have lots more smoothing parameters to estimate as well as the same number of model coefficients as your current model formulation (mins the ranef smooth) Commented Aug 30, 2023 at 14:36
• Understood, thank you, this is a huge help. Commented Aug 30, 2023 at 14:46