I am studying patterns of bird abundances with certain habitat variables and how they vary over time. I am interested in using GAMs with smooth terms for some of the variables. I am, however, confused when it comes to model selection here. So far, I went with the default in mgcv::gam (using gam(..., select=TRUE), which uses a GCV fitting procedure, but some reading pointed out to me that this might not be the best choice, especially for models that include both smooth and non-smooth terms (e.g. this answer or this one). This is because the default GCV fitting procedure does not penalise non-smooth terms enough and hence could bias the estimation towards linearity and undersmoothing. In such cases, REML is offered as a better solution.
Now, here is where I am confused. I am aware of REML and ML in the context of mixed models, with REML being useful for comparing models with the same fixed effects but different random effects (see this or this; this in the context of GLS). But I am unsure what exactly they mean in the context of GAMs? I wish to understand what the differences (and similarities) are, and which cases call for one over the other. I believe that the GCV computes prediction error from multiple versions of the data (i.e. the usual cross-validation procedure) while REML compares the likelihood with a whole random distribution (in my case, Poisson) and this is why the REML does better? Am I correct, or am I mixing up different things? And in this case, the GCV predictions should approximate the ones from REML with higher sample sizes?
Finally, adding to on to the last point, I understand that in the case of mixed models, the different methods converge at high sample sizes. Is this true in the context of GAMs as well?
(I have tried reading the ?gam.selection documentation, but I'm having trouble wrapping my head around all of it. But the documentation suggests additionally that ML methods for smoothness selection are better for inference using the approximate p-values provided in the summary.gam() function. And lastly, I think I am indeed mixing up different things, namely term selection and smoothness selection, so would it be better to build multiple models with just single-term differences and compare using AIC or AICc as suggested by the author of the package elsewhere?)
