¿Por qué el Criterio de Información de Akaike (AIC) favorece a veces un modelo sobreajustado?

Question

Como ejercicio para adquirir experiencia práctica trabajando con criterios de selección de modelos calculé los ajustes de los datos de mpg en carretera frente a la cilindrada del motor a partir del tidyverse mpg conjunto de datos de ejemplo utilizando polinomio y B-spline modelos de complejidad paramétrica creciente que van de 3 a 19 grados de libertad (obsérvese que el degree o df en cada familia de modelos cuenta el número de coeficientes adicionales ajustados además de la intercepción).

Los resultados del procedimiento de ajuste se muestran en el siguiente gráfico. La línea azul muestra el resultado de la regresión prevista a través de un conjunto de valores de desplazamiento del motor espaciados uniformemente a lo largo de la gama del conjunto de datos de entrada, mientras que los puntos naranja-rojo muestran el resultado en los valores dentro del conjunto de datos original que se utilizaron realmente para derivar el ajuste:

A partir de alrededor de ~11 grados de libertad, las líneas azules (que son más visibles cuando no hay observaciones en el conjunto de datos de entrada) empiezan a mostrar signos clásicos de sobreajuste: se contonean y giran salvajemente, variando mucho más que los propios datos de entrada, de hecho en algunos casos se extienden hasta la región no-fiscal (es decir, cayendo hasta mpg negativo, que no tiene interpretación física). Además, las dos clases de modelos (polinómico frente a B-spline) muestran aleatoriedad en las ubicaciones de estas caídas y ondulaciones. Un gráfico adicional a continuación muestra las diferencias entre las dos familias de modelos frente al aumento del número de parámetros del modelo. Para modelos más simples con menos parámetros, la diferencia es uniformemente pequeña, normalmente < 1-2 mpg en todo el rango del conjunto de datos, mientras que para modelos con más parámetros la diferencia es grande, y generalmente se vuelve más divergente a medida que aumenta el número de parámetros:

Basándome en el sobreajuste aparente que puedo ver con un mayor número de parámetros del modelo ajustado, esperaría que la mayoría de los criterios de selección de modelos eligieran un modelo óptimo con < 10 coeficientes ajustados. Sin embargo, extraje los valores del criterio de información de Akaike (AIC) devueltos con cada uno de los modelos ajustados, y eso no es realmente lo que ocurre en este caso. En su lugar, los modelos más complejos con el mayor número de parámetros ajustados presentan los valores AIC más pequeños y, por lo tanto, se ven favorecidos por el AIC:

Edita: Basándome en la respuesta de otro colaborador, he modificado el gráfico anterior para mostrar tanto el AICc como el AIC. Como era de esperar, el uso de AICc en lugar de AIC da lugar a una corrección que es de hecho mayor para los modelos con un mayor número de parámetros, pero no lo suficientemente grande como para hacer ninguna diferencia en el resultado final:

Mi pregunta: ¿qué está pasando aquí? ¿Por qué el AIC da un resultado contrario a la intuición, favoreciendo aparentemente a los modelos sobreajustados? ¿Existe algún criterio alternativo bien establecido que pueda funcionar mejor en este caso, seleccionando un modelo menos complicado que no presente un sobreajuste tan obvio?

Como referencia, el código que produce estos gráficos está aquí ( editar: actualizado para producir la versión 2 del gráfico AIC frente a los grados de libertad de entrada):

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

Answer 1

El AIC es sensible al tamaño de la muestra utilizada para entrenar los modelos. Con tamaños de muestra pequeños, "existe una probabilidad sustancial de que el AIC seleccione modelos que tengan demasiados parámetros, es decir, que el AIC sobreajuste". [1] En este caso, la referencia sugiere el AICc, que introduce una penalización adicional por el número de parámetros.

Esta respuesta de Artem Kaznatcheev sugiere un umbral de $n/K < 40$ como punto de corte para utilizar o no el AICc, basándose en Burnham y Anderson. Aquí $n$ significa el número de muestras y $K$ el número de parámetros del modelo. Sus datos tienen 234 filas disponibles (que figuran en la página web que ha enlazado). Esto indicaría que el punto de corte se sitúa aproximadamente en 6 parámetros, a partir de los cuales debería considerar el AICc.

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

Answer 2

Descargo de responsabilidad: No he revisado tu código línea por línea. A primera vista, parece legítimo, así que asumiré que lo es.

AIC no es más que la log-verosimilitud penalizada por el número de parámetros $k$

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

$2$ es una constante . $\ln(\hat L)$ es una suma de las evaluaciones no normalizadas de la función de verosimilitud sobre todos los puntos de datos. $2\ln(\hat L)$ puede ser realmente lo que sea y no hay garantías de que $2$ es la ponderación "adecuada" para el número de parámetros, de forma que le permita elegir el mejor modelo. Lo mismo se aplica al BIC y a todos los demás criterios de este tipo, funcionan bajo una serie de supuestos y tienden a funcionar bien en muchos casos, pero no hay ninguna garantía de que cosas como $2k$ son la penalización que funcionaría para todos los modelos posibles.

Answer 3

El problema del AIC es que no tiene en cuenta la estocástica del vector de parámetros ${\boldsymbol { \beta}}$ . Recordemos que en la regresión múltiple, cada estimación de los parámetros de regresión $\beta_0,\ldots,\beta_p$ siguen la distribución $\;\hat{\beta_j} \; \sim T(n-p-1)$ . Aquí $n$ es el número de puntos de datos y $p$ el tamaño del vector de parámetros (más uno, la constante $\beta_0$ ).

En concreto, define la inversa al cuadrado de la matriz de datos: $\bf{W} = (\bf{X}^T\bf{X})^{-1}$ y elemento $w_{j,j}$ siendo el $j$ elemento diagonal en $\bf{W}$ . Dejemos que $s=SRS/(n-p-1)$ donde SRS es la suma de los residuos al cuadrado: $SRS=(\bf{y}-\bf{X}^T{\boldsymbol { \beta}})^T(\bf{y}-\bf{X}^T{\boldsymbol { \beta}})$ . El vector $\bf{y}$ contiene los valores que se intentan predecir. Por último, $t_j=\hat{\beta}_j/(\sqrt{s}\,\sqrt{w_{j,j}})$ que es el estadístico de prueba para el parámetro $\hat{\beta}_j$ .

Claramente, la varianza normalizada de los grados de libertad $s$ no se tiene en cuenta en el AIC. Por eso este criterio de complejidad del modelo no es adecuado para la selección de modelos. Le sugiero que busque en la literatura un criterio de complejidad de modelos más avanzado. Al final, lo que quiere es hacer la elección óptima entre el sesgo del modelo y la varianza del modelo.

En lugar del AIC, se recomienda utilizar un criterio de selección de modelos que tenga en cuenta la incertidumbre del ajuste del modelo. Tomado de [1], el $SIC_f$ está bien definida y ni siquiera es compleja de calcular para un modelo de regresión lineal. Defina $SIC_f$ como \begin{equation} SIC_f = (n-p-2) \ln(s) + \ln \mid \bf{X}^T\bf{X} \mid \end{equation} donde $s$ es la varianza residual estimada, tal como se ha definido anteriormente y $\mid \; \cdot \; \mid$ es el determinante de la matriz.

Así que la propuesta es comparar $SIC_f(Model_a)$ con $SIC_f(Model_b)$ y elija el modelo preferido $a$ o $b$ . En $Model_{z}$ que tiene el más bajo $SIC_f$ es preferible. Recomiendo la lectura del referido artículo al respecto.

1. A. A. Neathy, J. E. Cavanaugh, Regression and time series model selection using variants of the Schwarz information criterion. Comunicaciones en Estadística , 26(3), 1997, p. 559-580.

Answer 4

No tengo una respuesta para ti, pero sí algunas cosas a tener en cuenta.

En primer lugar, hay algunas filas duplicadas en los datos que está utilizando, por lo que obviamente el sobreajuste no es un problema tan grande en ese caso.

Sin embargo, eliminar las filas duplicadas no resolverá realmente el extraño comportamiento del 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')

Ahora gam seleccionará automáticamente df alrededor de 5 cuando le damos el mismo problema, con el ajuste que parece bastante suave (es splines cúbicos en lugar de bsplines, pero eso no debería importar)

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

y cuando comparamos diferentes complejidades manualmente usando AIC, seleccionamos 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))

¿Por qué ocurre esto? No lo sé. Se dice que el AIC prefiere modelos complicados, pero yo no esperaría que fallara tan espectacularmente en un modelo tan sencillo. Se supone que se aproxima al rendimiento fuera de muestra, o LOOCV, por lo que los modelos con demasiados parámetros podrían no funcionar tan mal si se dispone de suficientes datos. Su modelo spline sobreajustado parece que funcionaría bien fuera de la muestra, aunque no perfectamente, lo que definitivamente no es el caso de su modelo polinómico. Al fin y al cabo, es sólo una aproximación, y existen muchas otras aproximaciones y criterios de selección de modelos. Por ejemplo, para los modelos lineales, se puede calcular el LOOCV sin necesidad de reajustar los modelos. Lo sorprendente es que en algunos campos, el AIC es el número más importante que se utiliza para juzgar sus modelos.