This is my first post & I'm fairly new to GAMs; apologies.
I ran a series of 16 generalized additive negative binomial models (gam, family=nb, mgcv package) with increasing complexity that modeled the abundances of 12 fish species (independently) as functions of several environmental metrics (and their combinations):
It is a large dataset; a couple examples of models:
m4(s): ccount ~ s(sal, k=6, bs="tp") m12(ost): ccount ~ s(do, k=6, bs="tp") + s(sal, k=6, bs="tp") + s(temp, k=6, bs="tp")
I conducted a bunch of model checks, then plotted the estimated degrees of freedom (edf), deviance explained, and AIC for each model (in order from least to most complex, approximately) to assist with model selection for each species.
It appears that deviance explained is simply the inverse of AIC, and AIC does not seem to be penalized by variation in the model's EDF. I plotted deviance explained vs AIC and was surprised that they were almost perfectly correlated. This seemed counter intuitive.
AIC vs. Dev.Expl
Question/Commment 1 - I found three different ways to call AIC:
1) model$aic 2) AIC(model) 3) extractAIC(model)
(3) is documented only for parametric models (i.e., not GAM). It's unclear whether (1) or (2) is more appropriate for gam objects based on the mgcv pkg? Both gave practically the same result (see below), but appeared to use different df (edf vs reference df). I assumed (1) would be the best; however, it was unclear what df were used, and there is documentation of logLik used with (2) for gam objects. It was unclear if this is already built into (1)? https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/logLik.gam.html
Question2 - is it reasonable to expect deviance explained and AIC to be this strongly correlated, regardless of model complexity? I didn't think so. Are improvements in dev.expl large enough that changes in edf are inconsequential? Is deviance explained already penalized for the edf in the model?
Thanks for any help.
Here is a simple reproducible example from mtcars showing the same pattern:
library(mgcv) # 3 GAM models b1 <-gam(mpg~s(hp), data=mtcars) b2 <-gam(mpg~s(wt), data=mtcars) b3 <-gam(mpg~s(hp)+s(wt), data=mtcars) # AIC1 a1.1 <- b1$aic a1.2 <- b2$aic a1.3 <- b3$aic a1 <- c(a1.1,a1.2,a1.3) # AIC2 a2.1 <- AIC(b1) a2.2 <- AIC(b2) a2.3 <- AIC(b3) a2 <- c(a2.1, a2.2, a2.3) # dev.explained d1 <- summary(b1)$dev.expl d2 <- summary(b2)$dev.expl d3 <- summary(b3)$dev.expl d <- c(d1,d2,d3) par(mfrow=c(2,2)) plot(a1~d, type=c("b"), xlab="deviance explained", ylab="model$aic") plot(a1~a2, col="black", type="b", lwd=1, cex=1, xlab = "AIC(model)", ylab="model$aic")
Following the helpful comment by Jon (below), I explored further how gam terms are being calculated in mgcv using the 2-factor gam for mtcars (b3, above). I was able to calculate deviance explained and AIC, which matched model outputs (though calculation of df seemed to double-count the intercept to get the same df as the gam output; odd).
The only calculation I could not get to work was Deviance. The model provided suggests deviance=2logL(model|data)−constant). I, assuming the same constant 2*p as for AIC, could not produce the same result as the gam model output. This only proves that, though grasping AIC and dev.expl, I still don't understand how Deviance is calculated and its relationship to AIC.
> # GAM based on mtcars > library(mgcv) > b <-gam(mpg~s(hp)+s(wt), data=mtcars) > summary(b) Family: gaussian Link function: identity Formula: mpg ~ s(hp) + s(wt) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 20.0906 0.3723 53.97 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(hp) 2.409 3.016 6.305 0.00218 ** s(wt) 2.085 2.523 15.277 1.19e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.878 Deviance explained = 89.6% GCV = 5.3542 Scale est. = 4.4349 n = 32 > plot(b) > # LOG LIKELIHOOD OF GAM MODEL (AND DF) > LL <- logLik(b) > LL 'log Lik.' -66.22412 (df=6.494385) > # MODEL DF = no. of parameters (p) = sum of smooth (edf) + parametric (nsdf) terms > # note: intercept gets counted twice (in b$edf and b$nsdf) to attain same df as logLik(b) and aic = AIC(b) > # note: summary(b)$pTerms.pv(does not include intercept as a parameter term; gives zero result for unspecified param param) > p <- sum(b$edf)+sum(b$nsdf) > p  6.494385 > # AIC = penalty - twice the log likelihood (2p-LL) > 2*p - 2*LL 'log Lik.' 145.437 (df=6.494385) > AIC(b)  145.437 > b$aic  145.437 > #DEVIANCE EXPLAINED = 1- (residual dev / total (null) deviance) > 1-(b$deviance/b$null.deviance)  0.895608 > summary(b)$dev.expl  0.895608 > #DEVIANCE = the unexplained error in the model (residual deviance) > 2*LL - 2*p # DOES NOT MATCH b$deviance! 'log Lik.' -145.437 (df=6.494385) > b$deviance  117.5503 >
Any further thoughts would be greatly appreciated.