18
$\begingroup$

As an exercise to develop practical experience working with model selection criteria, I computed fits of the highway mpg vs. engine displacement data from the tidyverse mpg example data set using polynomial and B-spline models of increasing parametric complexity ranging from 3 to 19 degrees of freedom (note that the degree or df number in either model family counts the number of additional fitted coefficients besides the intercept).

The results of the fitting procedure are shown in the plot below. The blue line shows the predicted regression result across an evenly spaced set of engine displacement values along the range of the input data set, while the orange-red points show the result at the values within the original data set that were actually used to derive the fit:

Fitted polynomial and B-spline models with tidyverse mpg example data

Starting at around ~11 degrees of freedom, the blue lines (which are most visible where there were no observations within the input data set) begin to exhibit classic signs of overfitting: they wiggle and gyrate wildly, varying much more widely than the input data itself, indeed in some cases extending down into the non-phyiscal region (i.e., dipping down into negative mpg, which has no physical interpretation). Moreover, the two model classes (polynomial vs. B-spline) exhibit randomness in the locations of these dips and wiggles. An additional plot below shows the differences between the two model families vs. increasing number of model parameters. For simpler models with fewer parameters, the difference is uniformly small, usually < 1-2 mpg across the entire range of the data set, while for models with more parameters the difference is large, and generally becomes more divergent as the number of parameters increases:

Differences (polynomial - B-spline) between models with increasing numbers of parameters fitted to tidyverse mpg data set

Based upon the apparent overfitting that I can see with higher numbers of fitted model parameters, I would expect most model selection criteria to choose an optimal model as having < 10 fitted coefficients. However, I extracted the Akaike Information Criterion (AIC) values returned with each of the fitted models, and that's not actually what happens in this case. Instead, the most complex models with the largest number of fitted parameters exhibit the smallest AIC values, and are therefore favored by AIC:

Akaike Information Criterion (AIC) for models with increasing numbers of parameters

Edit: based upon another contributor's reply, I've modified the above plot to show both AICc as well as AIC. As expected, using AICc instead of AIC results in a correction which is indeed larger for models with a greater number of parameters, but not nearly large enough to make any difference in the final outcome:

Akaike Information Criterion (AIC and AICc) for models with increasing numbers of parameters

My question: what is happening here? Why does AIC give a counterintuitive result, apparently favoring overfitted models? And is there any well-established alternative criteria that might be expected to work better in this case, selecting a less complicated model that does not exhibit such obvious overfitting?

For reference, the code that produces these plots is here (edit: updated to produce version 2 of the AIC vs. input degrees of freedom plot):

library(tidyverse)
library(splines)
library(AICcmodavg)

# ---- Part 1: Fit various models to data and plot results 

# Select degrees of freedom (i.e., number of coefficients, 
# not counting
# the intercept coefficient) for polynomial or B-spline models
dflo <- 3
dfhi <- 19
inputdf <- seq(dflo,dfhi,2)
ndf <- length(inputdf)
dflist <- as.list(inputdf)
names(dflist) <- sprintf("%2d", inputdf)

# Fit all of the models, and save the fit objects 
# to a nested list
# (outer list level: model family (poly or spline),
#  inner list level: dof)
fitobj <- list()
fitobj$poly <- map(dflist, ~lm(formula=hwy~poly(displ, 
                   degree=.), data=mpg))
fitobj$bspl <- map(dflist, ~lm(formula=hwy~bs(displ, df=.), 
                   data=mpg))

# Make a list of data points at which to predict the 
#fitted models (grid:
# evenly spaced 1D series of points ranging from minimum engine 
# displacement
# in the original data set to maximum engine displacement, 
# orig: the input
# data set itself)
ngrid <- 200
dval <- list()
dval$grid <- data.frame(displ=seq(min(mpg$displ), max(mpg$displ),  
                              length.out=ngrid))
dval$orig <- data.frame(displ=mpg$displ)

# Key elements for a new list
mtype <- list("poly"="poly", "bspl"="spline")
dtype <- list("grid"="Evenly spaced grid", 
              "orig"="Original values")

# For both models (poly and spline), and both sets of prediction points (grid
# and original), predict all the models, and dump all of the results to a list 
# of data frames, keyed by the cross product of the 2 pairs of input keys
pred <- list()
for(mkey in c("poly", "bspl")) {
    for(dkey in c("grid", "orig")) {
        crosskey <- paste(mkey, dkey, sep="_")
        pred[[crosskey]] <- map_dfr(map(fitobj[[mkey]], 
                                        predict.lm,
                                        newdata=dval[[dkey]],
                                        interval="confidence"),
                                    as_tibble, .id="inputdf") %>%
                            bind_cols(displ=rep(dval[[dkey]]$displ, ndf),
                                      modeltype=mtype[[mkey]],
                                      prediction_points=dtype[[dkey]])
    }
}
# Reorganize the list of data frames into a single giant unified data frame
dfpred <- map_dfr(pred, bind_rows) %>%
          mutate(modelspec=paste(inputdf, modeltype, sep=" "))

# Plot all of the fitted models. Evenly spaced 1D grid 
# is shown as a blue
# line, the original data points from the parent data set are 
# orange-red
# dots superimposed on top of it.  Gray ribbons are confidence 
# intervals 
# produced by predict.lm in previous step above.
plt_overview <- ggplot() +
                geom_ribbon(aes(x=displ, ymin=lwr, ymax=upr), 
                            data=dfpred,
                            fill="grey70") +
                geom_point(aes(x=displ, y=hwy), data=mpg, 
                size=1.5) +
                geom_line(aes(x=displ, y=fit, 
                color=prediction_points),
                          data=filter(dfpred, 
                          prediction_points==dtype$grid),
                      size=2) +
            geom_point(aes(x=displ, y=fit, 
                       color=prediction_points),
                       data=filter(dfpred, 
                       prediction_points==dtype$orig),
                           size=3) +
                scale_color_manual(breaks=dtype, 
                values=c("blue", "tomato")) +
                xlab("Engine Displacment (liters)") +
                ylab("Highway Miles Per Gallon (mpg)") +
                coord_cartesian(ylim=c(0,50)) +
                facet_wrap(~modelspec, ncol=4) +
                theme(text=element_text(size=32),
                      legend.position="bottom")
png(filename="fit_overview.png", width=1200, height=1600)
print(plt_overview)
dev.off()

# ---- Part 2: Plot differences between poly / spline model families ----

# For each input degree of freedom, calculate the difference between the
# poly and B-spline model fits, at both the grid and original set of 
# prediction points
dfdiff <- bind_cols(filter(dfpred, modeltype=="poly") %>%
                      select(inputdf, displ, fit, prediction_points) %>%
                      rename(fit_poly=fit),
                    filter(dfpred, modeltype=="spline") %>%
                      select(fit) %>%
                      rename(fit_bspl=fit)) %>%
          mutate(fit_diff=fit_poly-fit_bspl)

# Plot differences between two model families
plt_diff <- ggplot(mapping=aes(x=displ, y=fit_diff, color=prediction_points)) +
            geom_line(data=filter(dfdiff, prediction_points==dtype$grid),
                  size=2) +
        geom_point(data=filter(dfdiff, prediction_points==dtype$orig),
                       size=3) +
            scale_color_manual(breaks=dtype, values=c("blue", "tomato")) +
            xlab("Engine Displacment (liters)") +
            ylab("Difference (Poly - B-Spline) of Fit Results (mpg)") +
            coord_cartesian(ylim=c(-10,10)) +
            facet_wrap(~inputdf, ncol=4) +
            theme(text=element_text(size=32),
                  legend.position="bottom")
png(filename="fit_diff.png", width=1200, height=800)
print(plt_diff)
dev.off()

# ---- Part 3: Plot Akaike Information Criterion (AIC) for all models ----

# Compute both AIC and AICc for both model families (polynomial and B-spline)
# and organize the result into a single unified data frame
sord <- list("AIC"=FALSE, "AICc"=TRUE)
aictbl <- map(sord,
              function(so) map(fitobj,
                               function(fo) map(fo, AICc, second.ord=so) %>%
                                 unlist())) %>%
          map_dfr(function(md) map_dfr(md, enframe,
                                       .id="mt"), .id="aictype") %>%
          mutate(modeltype=unlist(mtype[mt]),
                 inputdf=as.numeric(name)) %>%
          rename(aic=value) %>%
          select(inputdf, aic, modeltype, aictype)

# Plot AIC
plt_aic <- ggplot(aictbl, aes(x=inputdf, y=aic, color=aictype)) +
           geom_line() +
           geom_point(size=3) +
           xlab("Input Degrees of Freedom") +
           ylab("Akaike Information Criterion (AIC)") +
           labs(color="AIC Correction Type") +
           facet_wrap(~modeltype, ncol=1) +
           theme(text=element_text(size=32),
                 legend.position="bottom")
png(filename="aic_vs_inputdf.png", width=800, height=800)
print(plt_aic)
dev.off()
$\endgroup$
0
12
$\begingroup$

The AIC is sensitive to the sample size used to train the models. At small sample sizes, "there is a substantial probability that AIC will select models that have too many parameters, i.e. that AIC will overfit". [1] The reference goes on to suggest AICc in this scenario, which introduces an extra penalty term for the number of parameters.

This answer by Artem Kaznatcheev suggests a threshold of $n/K < 40$ as cutoff point for whether to use AICc or not, based on Burnham and Anderson. Here $n$ signifies the number of samples and $K$ the number of model parameters. Your data has 234 rows available (listed on the webpage you linked). This would indicate that the cutoff exists at roughly 6 parameters, beyond which you should consider AICc.

[1] https://en.m.wikipedia.org/wiki/Akaike_information_criterion#modification_for_small_sample_size

$\endgroup$
2
  • 1
    $\begingroup$ i think that the recommendation is to not bother and always use AICc $\endgroup$
    – rep_ho
    May 14 at 21:34
  • 1
    $\begingroup$ @CursedHotdog: Welcome to CrossValidated, and thanks for taking the time to answer my question. Responding to your suggestion, I've augmented the 3rd plot with an updated version so that it now shows both AICc as well as AIC. As you described, the size of the resulting correction is indeed larger and more prominent for models with a greater number of adjustable parameters, however it is nowhere near large enough to make any difference in the final outcome: even AICc still favors an overfitted model. $\endgroup$
    – stachyra
    May 15 at 3:07
8
$\begingroup$

Disclaimer: I didn't go through your code line-by-line. At first sight, it seems legit, so I'll assume it is.

AIC is just the log-likelihood penalized by the number of parameters $k$

$$ 2k - 2\ln(\hat L) $$

$2$ is a constant. $\ln(\hat L)$ is a sum of the unnormalized likelihood function evaluations over all the data points. $2\ln(\hat L)$ can really be whatever and there are no guarantees that $2$ is the "appropriate" weight for the number of parameters so that it would enable you to pick the best model. The same applies to BIC and all the other criteria like this, they work under a set of assumptions and tend to work well in many cases, but there is no guarantee that things like $2k$ are the penalty that would work for every possible model.

$\endgroup$
6
  • 2
    $\begingroup$ I disagree that AIC is as arbitrary as you make it sound. The likelihood is well-defined up to a multiplicative constant, so the log-likelihood is well-defined up to an additive constant. That is, $\ln(\hat L)$ is in fact normalized in that it's a log to a specific base ($e$), and the additive change in $\ln(\hat L)$ from one model to another is absolute. $\endgroup$
    – nanoman
    May 16 at 7:41
  • 2
    $\begingroup$ @nanoman I never said it’s arbitrary, if it was it wouldn’t make any sense as a metric. I said that you have no guarantees that $2k$ is the one-size-fits-all penalty that is guaranteed to be appropriate for any model. “Up to a constant” is an important point if you put a constant in the equation. $\endgroup$
    – Tim
    May 16 at 8:17
  • 1
    $\begingroup$ “Up to a constant” is an important point if you put a constant in the equation. An additive constant in AIC drops out in comparing AIC for any set of models. The $2k$ term is not an additive constant. It says that an additional parameter is "worth" a factor of $e$ in likelihood. Indeed, this is not universally correct, but not for a trivial reason like being normalization-dependent as you imply. $\endgroup$
    – nanoman
    May 16 at 19:46
  • 2
    $\begingroup$ @nanoman I know that and not claiming that it is only about normalization. Just saying that this is a constant, while likelihood depends on data and other things, so the "weight" of the penalty may not be appropriate in every scenario. Nothing more than this. $\endgroup$
    – Tim
    May 16 at 19:51
  • 2
    $\begingroup$ Just to add to @tim compare to regularised regression with crossvalidation. Instead of 2k you would have $\lambda ||\beta||^2$, where the lambda is estimated from your data by crossvalidation $\endgroup$
    – seanv507
    May 20 at 6:38
6
+100
$\begingroup$

The problem with AIC is that it does not take into account the stochastics of the parameter vector ${\boldsymbol { \beta}}$. Recall that in multiple regression, each estimate of the regression parameters $\beta_0,\ldots,\beta_p$, follow the distribution $\;\hat{\beta_j} \; \sim T(n-p-1)$. Here $n$ is the number of data points and $p$ the size of the parameter vector (plus one, the constant $\beta_0$).

Specifically, define the squared inverse of the data matrix: $\bf{W} = (\bf{X}^T\bf{X})^{-1}$, and element $w_{j,j}$ being the $j$'th diagonal element in $\bf{W}$. Let further $s=SRS/(n-p-1)$, where SRS is the sum of the squared residuals: $SRS=(\bf{y}-\bf{X}^T{\boldsymbol { \beta}})^T(\bf{y}-\bf{X}^T{\boldsymbol { \beta}})$. The vector $\bf{y}$ contains the values you are trying to predict. Finally, $t_j=\hat{\beta}_j/(\sqrt{s}\,\sqrt{w_{j,j}})$, which is the test statistic for parameter $\hat{\beta}_j$.

Clearly, the degrees-of-freedom normalized variance $s$ is disregarded in the AIC. That is why this model complexity criterion is unsuited for model selection. I suggest that you look further in the literature for a more advanced model complexity criterion. In the end, you want to make the optimal choice between model-bias and model-variance.

Instead of the AIC, it is recommended to use a model selection criterion which accounts for the uncertainty of the model fit. Taken from [1], the $SIC_f$ is well-defined and not even complex to compute for a linear regression model. Define $SIC_f$ as \begin{equation} SIC_f = (n-p-2) \ln(s) + \ln \mid \bf{X}^T\bf{X} \mid \end{equation} where $s$ is the estimated residual variance, as defined above and $\mid \; \cdot \; \mid$ is the matrix determinant.

So the proposal is to compare $SIC_f(Model_a)$ with $SIC_f(Model_b)$ and choose the preferred model $a$ or $b$. The $Model_{z}$ that has the lowest $SIC_f$ is preferred. I recommend reading the referred article on this score.

  1. A. A. Neathy, J. E. Cavanaugh, Regression and time series model selection using variants of the Schwarz information criterion. Communications in Statistics, 26(3), 1997, p. 559-580.
$\endgroup$
6
  • 2
    $\begingroup$ "That is why this model complexity criterion is unsuited for model selection. I suggest that you look further in the literature for a more advanced model complexity criterion." O.K., but AIC (or its corrected version, AICc) is in fact widely used for exactly this purpose. If you have in mind an alternative criterion that you suspect might work better, please say what it is, and what's different or special about it that makes it less likely to overfit the data. $\endgroup$
    – stachyra
    May 17 at 18:24
  • $\begingroup$ In short, Fisher's Information Criterion (FIC) is a better suited alternative. See for example faculty.washington.edu/yenchic/17Sp_403/Lec3_MLE_MSE.pdf. FIC explicitly formulates a model variance and a squared bias term - these two terms add up, and a trade-off between them has to be made. So, I recommend FIC for model selection in multiple linear regression. See the extensive paper on this: A.A. Neathy, J.E. Cavanaughz, REGRESSION AND TIME SERIES MODEL SELECTION USING VARIANTS OF THE SCHWARZ INFORMATION CRITERION, Communications in Statistics, Vol. 26, 1997, p. 559-580. $\endgroup$ May 17 at 20:59
  • 2
    $\begingroup$ These articles suggest Fisher Information is not a model selection metric itself, but rather a correction term that may be applied to the Schwarz Information Criterion (SIC). It'd be helpful if you could point to a function or package (preferably R, but other frameworks like python/scipy are good too) which implements the selection criteria described in Neath & Cavanaugh. (Or if it's just a trivial one-liner, then edit your response to show how to implement within the context of the code sample provided in the question.) Otherwise your advice is probably not actionable by most future readers. $\endgroup$
    – stachyra
    May 18 at 18:27
  • $\begingroup$ As currently written, this answer appears to contain an error: the expression for $SIC_{f}$ does not seem to match the expression from eq. 2.9 of the cited reference (Neath & Cavanaugh), which includes a term $\ln{|I_{n}(_{k}\underline{\hat{\theta}}_{n},Y_{n})|}$ in which the item inside the log-determinant is identified elsewhere in the paper as representing the Fisher Information Matrix. By contrast, the closest-looking corresponding term in the answer, $\ln{|\mathbf{X^{T}X}|}$, actually represents the log-determinant of the data matrix--not the same thing. $\endgroup$
    – stachyra
    May 24 at 15:25
  • 2
    $\begingroup$ This answer states, "In the end, you want to make the optimal choice between model-bias and model-variance", which seems to suggest the contributor believes $SIC_{f}$ is a superior metric because it explicitly implements a bias-variance tradeoff, whereas AIC and AICc do not. This answer could be further improved by clarifying which term in $SIC_{f}$ represents the bias, and which the variance, and how we know how to identify each. Also, AIC/AICc is obviously trading off something; if it's not bias-variance, then what alternative trade is being made, and why is it conceptually inferior? $\endgroup$
    – stachyra
    May 24 at 15:43
3
$\begingroup$

I don't have an answer for you, but a few things to consider.

First, there are some duplicates rows in the data you are using, so obviously overfitting is not such a big issue in that case.

enter image description here

However, removing duplicate rows won't really solve the weird behavior of AIC

    library(tidyverse)
    library(splines)
    
    data(mpg)
    mpg_unique <- distinct(mpg[,c('displ', 'hwy')])
    
    df <- 1:19
    models_mpg_poly <- lapply(df, function(x) {lm(hwy ~ poly(displ, 
                                 x), data=mpg)})
    models_mpg_poly_unique <- lapply(df, function(x) {lm(hwy ~ 
                poly(displ, x), data=mpg_unique)})
    
    models_mpg_splines <- lapply(df, function(x) {lm(hwy ~ bs(displ, 
                             x), data=mpg)})
    models_mpg_splines_unique <- lapply(df, function(x) {lm(hwy ~ 
                  bs(displ, x), data=mpg_unique)})
    
    par(mfrow=c(2,2))
    plot(df, sapply(X = models_mpg_poly, FUN = AICc), 
             main='AICc poly', ylab = 'AICc')
    plot(df, sapply(X = models_mpg_poly_unique, FUN = AICc), 
          main='AICc poly unique rows', ylab = 'AICc')
    plot(df, sapply(X = models_mpg_splines, FUN = AICc), 
         main='AICc splines', ylab = 'AICc')
    plot(df, sapply(X = models_mpg_splines_unique, FUN = AICc), 
           main='AICc splines unique rows', ylab = 'AICc')

enter image description here

Now gam will automatically select df around 5 when we give it the same problem, with the fit that looking pretty smooth (it's cubic splines instead of bsplines, but that shouldn't matter)

    library(mgcv)
    gam_model <- gam(hwy ~ s(displ, k=20, bs='cr'), data=mpg_unique)
    summary(gam_model)
    plot(gam_model)

enter image description here

and when we compare different complexities manually using AIC, we select df = 6.

    models_mpg_gam <- lapply(1:19, function(x) gam(hwy ~ 
                s(displ, k=x, bs='cr', fx=T), data=mpg_unique))
    plot(sapply(X = models_mpg_gam, FUN = AICc))

enter image description here

So why does this happen? I don't know. It is said that AIC prefers complicated models, but I wouldn't expect it to fail as spectacularly in such a simple model. It is supposed to approximate out of sample, or LOOCV performance, so models with too many parameters might not actually perform that bad if you have enough data. Your overfitted spline model looks like it would work fine out of sample, although not perfectly, this is definitely not the case for your polynomial model. In the end, it's just an approximation, and there are many other approximations and model selection criteria out there. For example, for linear models, you can calculate LOOCV without the need to refit your models. What is striking is that in some fields, the AIC is the most important number that is used to judge your models.

$\endgroup$
4
  • 2
    $\begingroup$ Note that the result on the AIC approximating the LOOCV log score is an asymptotic result. It says nothing about the finite sample case where we blow model complexity up beyond all reasonableness. $\endgroup$ May 15 at 10:29
  • $\begingroup$ @StephanKolassa good point $\endgroup$
    – rep_ho
    May 15 at 10:59
  • $\begingroup$ @rep_ho: The de-duplicated gam result is interesting, and looks more like what I was expecting to see among the other model families. However, I'm not convinced the de-duplication procedure that got you there is entirely defensible. Many of the rows you're calling duplicates are in fact separate manufacturers and models; for example, the 1999 Chevrolet Malibu automatic transmission and 1999 Hyundai Sonata manual trans both had 2.4 liters displacement and got 27 mpg. In some sense, these are not truly duplicates, but rather, independent observations that both got the same measurement result. $\endgroup$
    – stachyra
    May 17 at 17:56
  • $\begingroup$ @stachyra you are right, in that case it's not justifiable to remove duplicates $\endgroup$
    – rep_ho
    May 21 at 14:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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