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?