Why does the Akaike Information Criterion (AIC) sometimes favor an overfitted model? 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:

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:

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:

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:

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

 A: 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.
A: 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.

*

*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.
A: 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.

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


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)


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


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.
A: 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
