Access AIC for Generalized Additive Mixed Models (GAMM) generated by r package `gamm4`

I am using the r pacakge gamm4 to run a gamm on a binary response varible with eight predictor variabels, two of which are smoothed using a thin plate spline. This model runs fine and I can varify it has generated results from summery()

Here is a subset of my data: structure(list(DATE = c("1/1/1990", "1/2/1990", "1/3/1990", "1/4/1990", "1/5/1990", "1/6/1990"), RESPONSE = c(0L, 0L, 0L, 0L, 0L, 0L), PREDICTOR1 = c(100L, 80L, 60L, 100L, 100L, 100L), PREDICTOR2 = c(5.3, 3.8, 2.5, 2.2, 2.8, 4), PREDICTOR3 = c(0.001016002, 0, 0, 0.001778004, 0.021590043, 0), PREDICTOR4 = c(315, 90, 90, 315, 360, 225), PREDICTOR5 = c(1019.326087, 1028.770833, 1027.920833, 1024.4625, 1022.25, 1018.891667), PREDICTOR6 = c(17.36086957, 12.1, 13.02083333, 16.5375, 16.48333333, 19.6375), PREDICTOR7 = c(13.47916667, 12.775, 12.525, 12.3875, 12.725, 13.09166667), PREDICTOR8 = c(0.28335247, 0.212789468, 0.239630838, 0.386550722, 0.221811468, 0.379079404), random_effect1 = c(1L, 1L, 1L, 2L, 2L, 2L),

Here is my model formula:

model_1 <- gamm(RESPONSE ~ PREDICTOR1 + PREDICTOR2 + PREDICTOR3 + s(PREDICTOR4) + 
           PREDICTOR5 + s(PREDICTOR6) + PREDICTOR7 + PREDICTOR8, data = data, 
           random = random_effect1, family = binomial())

However, when I attempt to run AIC on the gamm outcome object it fails

aic_model1 <- AIC(model_1)
Error in UseMethod("logLik") :
no applicable method for 'logLik' applied to an object of class
 "c('gamm', 'list')"

I have attemted to break down the AIC parts and calculate the fitted values using predict.gam and

deviance_model <- 
  -2 * sum(log(dbinom(model_1$RESPONSE, size = 1, prob = fitted_values)))

Both of these generate a NULL result. I also attempted AIC.gam and AIC.gamm those these failed. The issue (though I could be mistaken) appears to be the class of the object gamm, list. Does anyone know how to get around this issue? Thank you for your help in advance!

Solution

Firstly, you are not using gamm4() in your code. The package gamm4 provides a function with name gamm4(), not gamm().

If you're willing to use the AIC from the mixed model formulation (read ?mgcv::logLik.gam, what I show below yields the marginal likelihood-based AIC, which integrates out the random effects and hence the [wiggly bits of the] smooths) then you can apply the AIC() function to the $mer component of the object returned by gamm4()

# use example from ?gamm4::gamm4
set.seed(0) 
dat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth

## Now add 20 level random effect `fac'...
dat$fac <- fac <- as.factor(sample(1:20,400,replace=TRUE))
dat$y <- dat$y + model.matrix(~fac-1)%*%rnorm(20)*.5
br <- gamm4(y~s(x0)+x1+s(x2),data=dat,random=~(1|fac))

Now apply AIC() to br$mer such that the logLik() method for the mixed model class ("merMod") is used:

> AIC(br$mer)                                                                
[1] 1811.787

You can't do this with gamm() models with non_Gaussian families as those models are actually fitted via MASS::glmmPQL() and that functions does't maximise a likelihood (it uses a penalized quasilikelihood) and as such there is no likelihood and no AIC.

The point about the returned object containing multiple model fits remains; these functions return the mixed model form as one component and a GAM form as the second component. Use str(my_model, max = 1) to see the names of these two components as they differ between model fitting functions.