Comparing GAM models with/without fixed effect interactions using REML versus ML

Question

I have several GAM models fit with package mgcv that share the same smooths and random effects groups. I would like to compare support for whether interactions between the fixed effects are supported using AIC.

My understanding is that we should use maximum likelihood rather than REML to compare between models differing only in their fixed effects. However, fitting the GAMs using ML instead of REML to my dataset takes much, much longer (more than a day versus tens of minutes).

I'd like to know if it is a sign of something wrong that changing the fitting method from REML to ML takes orders of magnitude longer. For reference my dataset has 19,395 entries. Here is the code for the most complex model. I only change method "REML" to method = "ML".

level_bam7 <- bam(pa ~ scale(prev.rain) * level * roya_present + 
                        scale(days.elapse) +
                        s(month, k=5, bs="cc") +
                        s(month, level, k=5, bs="fs") +
                        s(plantID, k=582, bs="re") +
                        s(quadrat, k=127, bs="re") + s(year, k=4, bs="re"),
                family = binomial, 
                method = "REML", cluster = cl,
                data = roya_1_long)

Answer 1

bam() uses the default "fREML" for the method argument not "REML", which indicates to the software to use a higly optimised fast REML calculation. Using either "REML" or "ML" I would expect to increase model fit time over "fREML".

Why do you want to use AIC? If the aim is to fit a model and find the subset of terms that has the best predictive performance? If so, you can probably figure that out without having to use all 582 plants and all your data. Choosing a subset that is manageable with plain ol' gam() and running on a few random subsets would likely work best to suggest which model form worked best.

If the aim is to learn whether an interaction is important, I wouldn't both with AIC or even any formal p value computation. I'd just fit the model using "fREML" to the full data set and just see what the estimated value for the interaction term is and report the point estimate plus an uncertainty interval.