Huge ΔAIC between GAM and GAMM models

Question

I'm working with spatial fisheries catch data and environmental variables, and I'm correlating the abundance in the catch to some oceanographic parameters. I'm using a Generalised Additive Model (GAM) and a Generalised Additive Mixed Model (GAMM) with one and two random effects (mgcv package in R), in particular:

m1 <- gam(kg ~ s(var1, k=10) + s(var2, k=10) + ... + offset(time), 
          family = Gamma(link=log), method=REML)

m2 <- gamm(kg ~ s(var1, k=10) + s(var2, k=10) + ... + offset(time), 
           family = Gamma(link=log), method=REML, random= list(vessel = ~ 1))

m3 <- gamm(kg ~ s(var1, k=10) + s(var2, k=10) + ... + offset(time), 
          family = Gamma(link=log), method=REML, 
          random= list(vessel = ~ 1, dayOfTheYear = ~ 1))

I get quite reasonable results in terms of residuals and covariates partial effects, but I get a huge difference in AIC between m1 and m2/m3, i.e. AIC(m1) = 58998, AIC(m2$lme) = 9410, and AIC(m3$lme) = 8751.

From what I understood, I should be able to compare GAM and GAMM from mgcv using AIC... Am I doing something wrong?

Thanks!

Answer 1

You shouldn't compare the AICs between objects fitted with different software. gam() is fitted via some fancy code fu in the mgcv package, whereas your gamm() fit is actually fitted via fancy code in the MASS (glmmPQL()) and then nlme (lme()) packages. It would be common for different constants to end up in the log likelihood.

When I read your question, I assumed you wanted to compare a GAM with no random effect to the same model with a random effect(s). To do that, fit the non-random effect model with gamm() too. For example, using @ACD's example data (from a now-deleted Answer):

set.seed(13)
x = rnorm(1000)
eff = as.factor(round(rnorm(1000)+5))
y = exp(x)*runif(1000)+as.numeric(eff)
plot(x,y)

gam_example = gamm(y ~ s(x), method="REML", family=Gamma(link="log"))
gamm_example = gamm(y ~ s(x), method="REML", random= list(eff = ~ 1), 
                    family = Gamma(link="log"))

Which gives two lines give:

> AIC(gam_example$lme)
    [1] -2.136317
    > AIC(gamm_example$lme)
[1] -1286.448

Hence there is strong support for the inclusion of the random effect and we can compare the AICs because they have been fitted via the same algorithm and code.

Technically, what @ACD shows (in a now-deleted Answer) is also incorrect. Whilst the two models are comparable in terms of both including a random effect for the eff variable, the AICs are not comparable because very different algorithms are used in the fitting:

## after running the code from @ACD's answer
> AIC(gam_example)
[1] 1721.197
> AIC(gamm_example$lme)
[1] -1286.448

The difference in AIC is meaningless; in this case, the gamm() is fitted using a penalised quasi likelihood criterion, where as the gam() is fitted using a standard REML (in this case) criterion.

Answer 2

Sorry for going necro on a relatively old question. I have just finished reading Simon Wood's book "Generalized Additive Models An Introduction with R" and I think I have the solution.

Generally, any penalized regression smoothers can be written as components of a mixed model (page 316). The smoothing parameters are treated as variance components and can be estimated by ML, REML, or PQL. In GAMM there is no need for explicit penalization as in GAMs, hence no GCV/UBRE needed to pick the best lambda (penalty coefficient). In GAMMs lambda is part of the variance components and is estimated in the model.

I suppose to perform a test on the effect of the random variable, you could compare two models

m2a <- gamm(kg ~ s(var1, k=10) + s(var2, k=10) + ... + offset(time), 
           family = Gamma(link=log), method=REML, random= list(vessel = ~ 1))

m2b <- gamm(kg ~ s(var1, k=10) + s(var2, k=10) + ... + offset(time), 
           family = Gamma(link=log), method=REML)

However, if you are to perform a generalized likelihood ratio test (GLRT), care must be taken because you are testing for sigma(vessel) which is on the boundary of the parameter space (variance cannot be negative). Only very small p-values should be taken as a sign of no vessel effect.

Another drawback of comparing your models m1 and m2 using REML is that m1 and m2 will probably have different fixed effects (although the formulas for fixed effects in R code are the same, my guess is that gam and gamm will handle these differently - you can look at the lists of the estimated effects and compare the fixed effects parameter names/numbers between the two models). If the fixed effects are not the same under REML, the models cannot be compared.