I'm fitting a GAM using mgcv::gam() with spatial smoothes, similar to the model in this example (data here). I've found that using AIC(model) gives a different result to model$aic. Why is this? Which is correct?
library('mgcv')
galveston <- read.csv('gbtemp.csv')
galveston <- transform(galveston,
datetime = as.POSIXct(paste(DATE, TIME),
format = '%m/%d/%y %H:%M', tz = "CDT"))
galveston <- transform(galveston,
STATION_ID = factor(STATION_ID),
DoY = as.numeric(format(datetime, format = '%j')),
ToD = as.numeric(format(datetime, format = '%H')) +
(as.numeric(format(datetime, format = '%M')) / 60))
knots <- list(DoY = c(0.5, 366.5))
M <- list(c(1, 0.5), NA)
m <- bam(MEASUREMENT ~
s(ToD, k = 10) +
s(DoY, k = 30, bs = 'cc') +
s(YEAR, k = 30) +
s(LONGITUDE, LATITUDE, k = 100, bs = 'ds', m = c(1, 0.5)) +
ti(DoY, YEAR, bs = c('cc', 'tp'), k = c(15, 15)) +
ti(LONGITUDE, LATITUDE, ToD, d = c(2,1), bs = c('ds','tp'),
m = M, k = c(20, 10)) +
ti(LONGITUDE, LATITUDE, DoY, d = c(2,1), bs = c('ds','cc'),
m = M, k = c(25, 15)) +
ti(LONGITUDE, LATITUDE, YEAR, d = c(2,1), bs = c('ds','tp'),
m = M, k = c(25, 15)),
data = galveston, method = 'fREML', knots = knots,
nthreads = 4, discrete = TRUE)
AIC(m)
[1] 57073.08
m$aic
[1] 57053.21
Note that the example I've given uses bam() instead of gam() but the result is the same.
I couldn't replicate this with a simpler model (example from here):
set.seed(2)
dat <- gamSim(1,n=400,dist="normal",scale=2)
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
AIC(b)
[1] 1696.143
b$aic
[1] 1696.143
The discrepancy is because what is stored in $aic takes as the degrees of freedom for the complexity correction in AIC is the effective degrees of freedom (EDF) of the model. This has been shown to be too liberal or conservative and can result in AIC always selecting the more complex model or the simpler model depending on whether a marginal or conditional AIC is used.
There are approaches to correct for this behaviour and mgcv implements the one of Wood et al (2016), which applies a correction to the degrees of freedom. This is done via the logLik.gam() function, which is called by AIC.gam(). This also explains the difference as $aic is the standard AIC without the correction applied and IIRC is a component of the GAM object that significantly predates the work of Wood et al (2016).
As to why you couldn't replicate this with the simple example, that's because the correction requires the use of components of the fit that are only available when the method used to fit is "REML" or "ML" (including "fREML" for bam() and also not when the Extended Fellner Schall or BFGS optimizers are used:
> library('mgcv')
Loading required package: nlme
This is mgcv 1.8-34. For overview type 'help("mgcv-package")'.
> set.seed(2)
> dat <- gamSim(1,n=400,dist="normal",scale=2)
Gu & Wahba 4 term additive model
> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat, method = 'REML')
> AIC(b)
[1] 1698.504
> b$aic
[1] 1696.894
where the default for method is to use GCV.
Wood, S.N., Pya, N., Säfken, B., 2016. Smoothing Parameter and Model Selection for General Smooth Models. J. Am. Stat. Assoc. 111, 1548–1563. https://doi.org/10.1080/01621459.2016.1180986