# Model selection and comparison in GAMM using R (mgcv)

I'm fitting a GAMM with correlation structure, using a non-Gaussian family. Here's an example of my global model: M0 <- gamm(response ~ var1*var2 + var3 + s(var4) + s(var5) + s(var6,var7), random=list(placeID= ~1), correlation= corAR1(form= ~ year | placeID), data=data, family=quasipoisson)

I'd like to do term selection on the global model to drop any nonrelevant variables, but I'm not sure how it should be done. Based on the help text of gam.selection - and hoping it applies to GAMM as well! - I've used backward stepwise selection so far. The text says "It is perfectly possible to perform backwards selection using p-values in the usual way". Thus, I've started from the smoothers (removing smoothers from any linear and/or nonsignificant variable) and then moved on to linear variables. I would be interested to hear, whether you think this method is useful in GAMM or not. Please, also tell if you have something else to suggest (e.g. GCV mentioned below).

I would've also liked to compare models that differ in their structure regarding coordinates (smoothed or not, interaction included or not). I suspect that s(latitude,longitude) might be overfitting, and I would like to check this somehow. Since the model uses PQL, I gather that AIC is not recommended (although included in lme part).

The gam.selection text mentions different scores: "In general the most logically consistent method to use for deciding which terms to include in the model is to compare GCV/UBRE/ML scores for models with and without the term." Is this method only for comparing models that differ by a single variable? If so, do you have any suggestions on how to compare models such as:

• M1 <- gamm(response ~ var1*var2 + s(lat,long)+ s(var5), random=list(placeID= ~1), correlation= corAR1(form= ~ year | placeID), data=data, family=quasipoisson)
• M2 <- gamm(response ~ var1*var2 + lat + s(var5), random=list(placeID= ~1), correlation= corAR1(form= ~ year | placeID), data=data, family=quasipoisson)

I had been under the impression that my global GAMM-model M0 uses GCV (the default), but when I call gam.check(M0$gam), I only get the four plots and no GCV score or other text output. Is there a way to acquire a GCV score of a GAMM object? ## 1 Answer You are not fitting a GAM but a GAMM. If you look at ?gamm you'll see that there are two options for model fitting an AMM and only one when fitting a GAMM:  method: Which of ‘"ML"’ or ‘"REML"’ to use in the Gaussian additive mixed model case when ‘lme’ is called directly. Ignored in the generalized case (or if the model has an offset), in which case ‘gammPQL’ is used.  As you are fitting with a quasipoisson family, the model fitting hierarchy is: gamm() ---> MASS::glmmPQL() ---> nlme::lme()  and the likelihood that is optimised is a penalised quasi-likelihood (PQL). As the splines are treated as part fixed and part random effects in this formulation of the GAMM, smoothness selection boils down to using PQL and I am not aware of results that can advise on how well PQL performs in that regard (unlike REML selection with regular GAMs). If you only need simple random effects, it might be possible to work in the GAM side of things: • In ?magic there is an example of fitting the GAM with correlated data (but in the Gaussian case. • In ?bam there is an option to specify an AR(1) parameter, but again, in the Gaussian case. You would need to use bs = "re" in an s() term for the random effect. • I had not properly acknowledged the differences in the underlying fitting methods of GAM vs GAMM. I infer from your reply, that gam.selection and gam.check are for GAM alone, not GAMM. From ?gamm, I now see that only predict, summary, print & vis.gam are available for gammobject$gam, ie. there's no information on GCV model selection. Hope I got it right now. So, from my original questions, it remains to be unclear whether the stepwise p-value-based selection is OK for GAMM. I do wish to include a correlation structure, which makes GAM with the bs-option not an ideal solution. – Mari L Apr 9 '15 at 18:45
• Does glmmPQL/gammPQL always optimise penalised quasi-likelihood? I am using a gamm with gamma family and I would like know the likelihood that is optimised in that case. – bee guy Feb 7 at 13:37
• If you fit using gamm() then yes, that model will be fitted via glmmPQL() and use PQL. If you fit using gamm4::gamm4() then ML will be used for this model. – Reinstate Monica - G. Simpson Feb 8 at 5:29