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.
Best Answer
They are the same thing; we just prefer the terminology "hierarchical" over "mixed", because the salient practical feature of these models is that they can model variation in the response that occurs at multiple levels, rather than the fact that they have fixed and random effects.