Random Effects Model – Why is There a Huge AIC Difference Between GAM and GAMM Models?

aicgeneralized-additive-modelrandom-effects-model

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!

Best Answer

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.

Related Question