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:


# 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)

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")

MTCARS GAM Example: AIC vs Dev.Expl & model$aic vs AIC(model) MTCARS GAM Example

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).

MTCARS GAM SMOOTH mtcars gam smooth

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 

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)


> 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
[1] 6.494385

> # AIC = penalty - twice the log likelihood (2p-LL)

> 2*p - 2*LL
'log Lik.' 145.437 (df=6.494385)
> AIC(b)
[1] 145.437
> b$aic
[1] 145.437

> #DEVIANCE EXPLAINED = 1- (residual dev / total (null) deviance)

> 1-(b$deviance/b$null.deviance)  
[1] 0.895608
> summary(b)$dev.expl
[1] 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
[1] 117.5503

Any further thoughts would be greatly appreciated.


2 Answers 2


The respective formulas for these two quantities are: $$\text{deviance} = 2\log\mathcal{L}(\text{saturated model}\, |\, \text{data}) - 2\log\mathcal{L}(\text{model}\, |\, \text{data})$$ $$\text{AIC} = 2k- 2\log\mathcal{L}(\text{model}\, |\, \text{data})$$ where $\mathcal{L}$ is the likelihood and $k$ is the number of model parameters. For a fixed dataset and model family, the saturated model is fixed, and therefore for our purposes the equation for deviance is: $$\text{deviance} = \text{constant} - 2\log\mathcal{L}(\text{model}\, |\, \text{data})$$

Plotting AIC against deviance the way that you've done, we expect the data to fall along a straight line if there exist constants $c_1$ and $c_2$ such that: $$c_1 \cdot \text{AIC} + c_2 \approx \text{Deviance}$$

This can only be the case if $k \propto \log\mathcal{L}$. Although this is not a relationship that I have previously come across, it seems plausible.

However it could also be that a different formula for Deviance is being used altogether, as intimated here.

  • $\begingroup$ Thank you, Jon. This suggests both are penalized by a constant (e.g., 2*df model) correct? I'm a bit confused as I've also read that AIC is "based on the Deviance, but penalizes you for making the model more complicated." This inquiry began with a reviewer insisting that AIC was penalized and deviance is not, and this seems to agree with gam() help which indicates model aic is and deviance is not penalized. Howerver, the results above support your comment, and/or that the aic penalty is fairly negligible here. Still reading & learning; thanks again for any thoughts. $\endgroup$
    – L.Steele
    Jan 31, 2018 at 6:58
  • $\begingroup$ @L.Steele My initial answer was mistaken, and I've updated to clarify. I would recommend reading this post for information on how deviance is calculated for GAMs: stats.stackexchange.com/questions/190172/… $\endgroup$
    – jon_simon
    Feb 1, 2018 at 0:27
  • $\begingroup$ Another thought: to determine whether the observed effect is due to unexpected behavior on the AIC side side or the Deviance side, replace AIC and Deviance, each in turn, by logLik and see what happens to the plots $\endgroup$
    – jon_simon
    Feb 1, 2018 at 4:30

I realize this is a late response, but the documentation for logLik.gam at: https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/logLik.gam.html

explains that to obtain an AIC corrected for smoothness parameter estimation, it's necessary to use method="REML" or method="ML" when calling gam(), which doesn't appear to be the case in the code snippets posted so far.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.