De acuerdo a Szretter (2017) “Una vez que se tienen todas las variables, es de interés contar con un criterio numérico para resumir la bondad del ajuste, que un modelo lineal con un cierto subconjunto de covariables da a la variable dependiente observada. A partir de este criterio se podrán ranquear los modelos y elegir un conjunto de unos pocos buenos candidatos para estudiar luego en detalle.”

En general en la regresión podemos establecer 3 valores para cada observación:

ŷ = el valor de la predicción de acuerdo al modelo

y̅ = el valor medio o esperado para las observaciones

y= el valor de la observación.

La figura 1 clarifica estos conceptos.

Figura 1

Relación entre SST, RSS, SSR

Donde:

SSTo =Suma de cuadrados totales,es la suma de los cuadrados de las desviaciones de todas las observaciones, y_i, de su media.

SSReg = Suma de cuadrados regresión, le indica cuánto de la variación en la variable dependiente explicó su modelo

SSRes = Suma de cuadrados de los residuos, es el valor que no es explicado por el modelo

Cuando tratamos de comparar modelos la primera aproximación sería mediante el R² obtenido para cada uno de los distintos modelos con distinta cantidad de covariables. La fórmula del coeficiente de determinación se muestra en la figura siguiente.

En la medida en la que el modelo sea más eficiente, la suma de cuadrado de los residuos va a ser más pequeña, llegando a 0 cuando el modelo es perfecto, en cuyo caso obtendríamos un R² sería igual a 1.

Resulta que comparar modelos usando el criterio de elegir aquél cuyo R² sea lo más grande posible equivale a elegir aquel que tenga la menor suma de cuadrados de residuos SSRes (ya que la suma de cuadrados total SSTo no depende de las covariables del modelo ajustado y por eso permanece constante).

Para mejorar el R² obtenido los analistas tratan de incluir una mayor cantidad de covariables, sean estas apropiadas para ajustar los datos o no. Es por esto, por lo que el criterio no es identificar el modelo con mayor R² (ese será siempre el modelo con todas las covariables disponibles) sino encontrar el punto a partir del cual no tiene sentido agregar más variables ya que estas no inciden en un aumento importante del R². Para encontrar este punto óptimo puede ser útil trazar un gráfico del R² en función de la cantidad de predictores. La siguiente figura explica el concepto.

Figura 3

R² obtenido en función de la cantidad de predictores

R² Ajustado

Para entender esta métrica, es importante indicar los grados de libertad de cada una de las sumas de cuadrados mencionadas anteriormente, esto nos ayudará a calcular la media de las sumas de cuadrados, es una forma de hacer un promedio que pueda ser comparable de una suma de cuadrados a otra. La siguiente tabla presenta los valores.

Tabla 1

Grados de libertad para las sumas de cuadrados

Como el R² no toma en cuenta el número de variables en el modelo de regresión, un criterio de decisión mucho más objetivo es calcular y comparar modelos por medio del R²a, para esto partimos de la fórmula de la figura 2 usando la media de las sumas de cuadrados. Entonces notamos que los grados de libertad del denominador no varía con respecto a la cantidad de coeficientes, mientras que el numerador si será afectado.

Por lo tanto, el R²a se calculará con la siguiente fórmula

Figura 4

Fórmula de R²a

Creación del autor Alfonso Prado

Esta sería la fórmula para un modelo sin intercepto donde p representa el número de predictores. Para modelos con intercept los grados de libertad del numerador se calcularía como n-p-1.

De acuerdo a Szretter (2017), entonces buscamos el subconjunto de p − 1 covariables que maximicen el R²a, o un subconjunto de muchas menos covariables para las cuales R²a el incremento sea tan pequeño que no justifique la inclusión de las covariables adicionales.

Veamos un ejemplo en la predicción de precios de vivienda

Del ejercicio anterior podemos notar que en la medida la cantidad de variables involucradas en cada modelo va aumentando la métrica de R²a aumenta, pero ya no es proporcional, en el modelo 4 tenemos 20 coeficientes y el incremento del R2a es insignificante.

El término likelihood a veces se lo traduce como probabilidad o verosimilitud. En general se refiere a la probabilidad de que una muestra provenga de una distribución dada. Recordando el concepto de la función PDF: Dada una distribución θ que tiene una media=μ y var= σ², la probabilidad de un valor específico es dada por la función de densidad de dicha distribución.

El likelihood se calcula como el producto de densidades de las observaciones. Cuando se hace una predicción se esperaría que el valor predicho caiga dentro de la distribución de la variable Y, en la medida en la que la probabilidad difiera de la distribución θ entendemos que el modelo no es bueno. Por ejemplo, suponga que usted tiene una variable de respuesta Y que sigue una distribución gaussiana θ como en la figura 5, suponga que usted está interesado en predecir cierto valor de la variable basado en determinadas características y al hacerlo encuentra que la salida Y tendría una probabilidad muy baja, obviamente deducimos que el modelo ha fallado en predecir adecuadamente.

Figura 5

Distribución de la variable Y y su predicción

Pero, como la probabilidad PDF para un valor dado es generalmente baja, el producto de estas va a ser muy muy bajo, por lo que se estila hacer una transformación logarítmica conocida como el “log-likelihood”. Para esto recordemos una de las propiedades de los logaritmos.

Figura 6

Propiedad de la multiplicación de logaritmos

Si se usa logaritmos naturales (base e) obtendríamos las fórmulas del likelihood como mostramos en la siguiente figura, donde π representa la función de producto de las probabilidades.

Figura 7

Fórmulas likelihood como producto (a) y log-likelihood como suma (b). Prado ( 2025)

Adaptación de la fórmula para la resolución de likelihood Faraway (2014)

En otras palabras, responde a la pregunta: "¿Cuál es la probabilidad de observar estos datos, dados estos parámetros?

Estimadores

Existen otros estimadores que nos ayudan a visualizar el mejor modelo basado en el total de las K covariables (el más grande posible) que incluya las covariables importantes de modo que en un modelo completo la estimación de la varianza del error sea insesgada.

Los principales estimadores son:

AIC : Criterio de Akaike

BIC : Criterio Bayesiano de Schwartz

CP: Cp de Mallows

Todos intentan resolver este problema introduciendo un término de penalización para el número de parámetros en el modelo.

Desafortunadamente diferentes estimadores pueden recomendar diferentes modelos. Sin embargo, tomados en conjunto estos criterios permiten identificar un conjunto pequeño de modelos de regresión que pueden ser construidos a partir de las variables independientes relevadas. (Szretter, 2017)

AIC y BIC

Cuando se usa un modelo para representar el proceso que generó los datos, la representación casi nunca será exacta, por lo que se perderá cierta cantidad de información al usar el modelo para representar el proceso.

De acuerdo a Faraway (2020) , AIC y BIC estiman la cantidad relativa de información perdida por un modelo dado, siendo el mejor aquel cuya pérdida de información sea la menor.

Al ajustar modelos, es posible aumentar la probabilidad agregando variables, pero hacerlo puede provocar un sobre ajuste. Los estimadores intentan resolver este problema introduciendo un término de penalización proporcional al número de variables en el modelo.

Entonces dado una familia de modelos M_α y una cantidad de observaciones, queremos saber cuál es el mejor modelo, donde α representa el número de predictor(es), ejemplo modelo 1 usa el primer predictor, modelo 2 usa los dos primeros predictores, etc. y donde los posibles predictores se ordenan por correlación con la variable de respuesta.

AIC y BIC parten del likelihood (L(θ)), pero como vemos en la siguiente tabla sus fórmulas son diferentes.

Tabla 2

Fórmulas matemáticas y funciones programáticas de AIC y BIC

Donde:

L=likelihood p=# variables n=# muestras

K= la penalización por parámetro a utilizar; el valor predeterminado k = 2

Dado una cantidad de modelos con distinta cantidad predictores podemos obtener la siguiente gráfica.

Figura 8

Relación de los estimadores en función de predictores

Continuando con nuestro ejemplo de predicción de precios de vivienda

Figura 9

Comparación de 5 modelos con estimadores AIC y BIC

Note que en este caso los 2 estimadores recomiendan el mismo modelo, aunque sus valores son distintos.

De la tabla de los estimadores vale la pena aclarar el termino df. En primer lugar, df significa grados de libertad, que en este caso se traduce a la cantidad de variables predictoras. Si se fija en la creación del modelo 1, este toma 3 variables Units +SqFt +Boro, las 2 primeras son variables continuas y la tercera es una variable categórica que tiene 5 niveles, los cuales han sido automáticamente codificados mediante one-hot, por lo tanto, glm ha creado 5 variables codificadas para la variable Boro. Además, tenemos el intercept del modelo con el cual suma o considera que el modelo tiene 8 variables.

En el modelo 2 tenemos, además de las anteriores, una variable que incluye el efecto combinado de Units:SqFt. Y así sucesivamente con el resto de modelos. Estos grados de libertad vienen a ser la cantidad de variables sobre las cuales se hará la penalización en ambos estimadores.

Selección automática

De los ejercicios presentados anteriormente entendemos que existen diferentes formas de evaluar la bondad de un modelo, sin embargo, es interesante notar que esto ha requerido una gran cantidad de tiempo de parte del analista, por lo que es importante conocer de mecanismos automáticos que nos permitan ubicar los mejores modelos en una forma oportuna.

Para esto existen 4 enfoques:

• Todos los subconjuntos posibles
• Eliminación hacia atrás
• Selección hacia adelante incorporando variables
• Regresión de a pasos

Todos los conjuntos posibles

Si el dataset tiene p predictores teóricamente se podría crean 2^p modelos. Mediante este algoritmo crearemos efectivamente esta cantidad de modelos y evaluará mediante cualquier métrica descrita anteriormente cual es el subconjunto óptimo. Por supuesto, esto significa que se realizarán 2^p regresiones, lo cual podría ser problemático para valores de p > 10.

Eliminación hacia atrás

En este algoritmo, iniciamos con un modelo que contiene todas las covariables, y analizamos los valores t y pvalue para cada coeficiente, si todos los pvalues < α, entonces el modelo completo es el mejor, caso contrario elimina la variable que tenga el mayor pvalue, manteniendo las restantes en el modelo. Se evalúa las métricas del modelo. Si el retiro de dicha variable ha provocado un descenso considerable de la métrica entonces el último modelo sería el mejor.

Selección hacia adelante

Es el inverso del anterior, iniciamos con un modelo con una variable y vamos incorporando variables adicionales hasta el punto en el cual las ganancias en el valor de la métrica ya no son significativas. En este caso la métrica que se utiliza es el estadístico F que nos muestra la relación del modelo con respecto a otro que solo tiene el intercept. Recuerde que este estadístico debe compararse con un Fcritico y cuando esta relación supera el valor α el proceso se detiene

Paso a paso

Es una modificación del procedimiento hacia adelante, que elimina una variable en el modelo si ésta pierde significancia cuando se agregan otras variables. El algoritmo es similar al de selección hacia adelante excepto que, a cada paso, después de incorporar una variable, el procedimiento elimina del modelo las variables que ya no tienen significancia.

Es importante mencionar que no siempre el algoritmo hacia adelante y el algoritmo hacia atrás terminarán seleccionando el mismo modelo. Note que la selección se hace en base a AIC.

Veamos un ejemplo de movilidad con mtcars

Para entender la presente clase debemos profundizar en los siguientes conceptos vistos en clases anteriores. Específicamente, habíamos tratado el concepto del estadístico R², como una medida de la bondad del modelo.

Supongamos por un momento que cada punto y_i estuviera muy cerca de la media y̅, esto significaría que cada y_i no dependería de la variación de x_i, y que tampoco habría mucho error aleatorio (ε) en el valor. Dado que esperamos que esto no sea así, entendemos cuánto contribuyen la predicción de x_i y el error aleatorio a y_i. En particular, observemos la distancia entre y_i y la media y̅. Escribiremos esta diferencia como SSTo.

En particular, el residuo se define como y_i - ŷ_i, que es la distancia desde el punto de datos original hasta el valor predicho en la línea de regresión. Se puede considerar como el error restante después de que el modelo haya realizado su trabajo. Esta diferencia se muestra gráficamente en la figura siguiente.

Figura 1

Relación entre SSTo, SSres, SSReg

De la figura podemos deducir que para cada punto i (y - y̅) =( ŷ - y̅) + (y - ŷ). El primer término del lado derecho puede entenderse como el valor explicado por el modelo , mientras que el segundo término puede entenderse como un valor que no es explicado por el modelo o el error cometido por el modelo. Si sumamos todas las diferencias y las elevamos al cuadrado (para evitar el neteo de valores positivos y negativos), obtenemos la fórmula de la figura 2a y poniendo la misma en relación con SSto la figura 2b.

Pero SSreg/SSto es precisamente la definición de R2 o coeficiente de determinación mencionado en clases anteriores. Aquí vemos porque R² puede interpretarse como la fracción de variabilidad en los datos que explica el modelo, por lo tanto, ahora tenemos una definición distinta de que es R², en realidad representa cuanto de la variable de respuesta está explicada por las predictoras y cuanto no. Cuando utilizamos el algoritmo de mínimos cuadrados ordinarios (MCO) lo que éste trata de obtener son los coeficientes (Betas) que minimizan el segundo componente (SSres/SSto).

Por otro lado, en los modelos lineales generalizables (GLM) su métrica principal para explicar la relación entre variables predictoras y de respuesta es conocida como “deviance”, (desviación o desviancia en español) es una medida de la bondad del ajuste, cuanto menor sea la desviación, mejor será el ajuste. Más específicamente podemos decir que si, para un modelo glm en particular, denotamos L_m como el máximo likelihood (o verosimilitud) alcanzable bajo este modelo, y a L_s como el likelihood del "modelo saturado" (que se define como el modelo perfecto cuya predicción se ajusta exactamente a las observaciones), entonces el deviance residual se calcularía como el doble negativo del logaritmo de la relación entre L_m y L_s . Y generalizamos diciendo que este modelo también tiene la máxima verosimilitud alcanzable entre todos los modelos posibles dentro del marco glm.

Figura 3

Fórmula de desviación en glm

Creación del autor: Alfonso Prado

Entonces la desviación puede verse como una generalización de la suma de cuadrados residuales de los modelos lineales.

Por otro lado, respecto a anova, habíamos indicado que existen 2 tipos de pruebas en las que usamos este concepto. La prueba conocida como anova de 1 vía cuya función es aov() la cual nos permite contestar la siguiente pregunta: ¿La varianza entre las medias de dos poblaciones es significativamente diferente?", y la segunda que, es de nuestro interés, es la función anova() que nos permite obtener una tabla de análisis de varianza de uno o más objetos de la clases lm o glm incluyendo una prueba F. A continuación, presentamos la sintaxis de anova.

Figura 4

Sintaxis de la función anova

En la función anova(), al especificar un solo objeto, se obtiene una tabla de análisis de varianza secuencial para ese ajuste.

Por ejemplo, en clase anteriores habíamos generado un modelo para explicar los valores de viviendas en NY en función la cantidad de unidades habitacionales, los pies cuadrados de la misma y el vecindario (valxSqFt ~TUnits+ GSqFt +Boro ). Aplicando anova() obtenemos la siguiente información:

Cada fila de la tabla corresponderá a un término del modelo, más una fila adicional en la parte superior para el modelo nulo (solo considera intersección).

Df (grados de libertad): Los grados de libertad asociados a cada término, normalmente 1 en caso de variable numéricas y numero de niveles -1 en caso de variable nominales.

Deviance: La desviación, una medida del ajuste del modelo a los datos para cada término.

DF residual (grados de libertad residuales): Los grados de libertad residuales tras añadir el término.

Deviance Residual (desviación residual): La desviación residual tras añadir el término.

Figura 5

Anova de 1 modelo logístico, la salida varia en modelos de regresión lineal

Es decir, las reducciones en la suma de cuadrados residuales a medida que se añade cada término de la fórmula (variables predictoras) y se presentan como las filas de la tabla, por último, presenta la suma de cuadrados residuales.

Esa tabla nos permite analizar cuál es el efecto de ir aumentando cada predictora, específicamente, sobre el ejemplo, vemos que al incluir la variable Boro el deviance residual cae de 136.05 a 64 indicado que esta variable ha mejorado la calidad del modelo.

Si se especifica más de un objeto como argumento en la función anova() ésta permite comparar distintos modelos. La tabla incluye los estadísticos usados en la prueba. Normalmente, la prueba F es la más adecuada, ya que compara la media cuadrática de una fila con la suma de cuadrados residuales del modelo más grande considerado, pero opcionalmente podrían especificarse otros como chi-cuadrado.

La prueba F se construye a partir de dos tipos de cantidades: sumas de cuadrados (SS) y grados de libertad (gl). Estos dos elementos definen un valor cuadrático medio (MS = SS/gl), y obtenemos nuestro estadístico F contrastando el valor de MS asociado con "todo lo demás" que serían los residuos. Esto nos permitirá establecer un marco para comparar modelos diferentes. En otras palabras, lo que queremos hacer es averiguar cómo describir el valor de SS asociado con la diferencia entre dos modelos.

Es importante entender cómo anova trabaja para comprender su salida. La función está hecha para comparar distintos modelos, pero estos deben ser incrementales (algunos autores lo llaman modelos anidados). Por ejemplo, el modelo 1 tendría 2 predictoras, el modelo 2 debe tener las mismas variables del modelo 1 más un incremento de variable(s). Entonces podemos decir que el modelo 1 está anidado dentro del segundo. Para cada modelo va a calcular los grados de libertad, que en este caso representa la cantidad de variables -1.

En el caso de comparar 2 modelos podemos despejar los SSreg para cada modelo de la siguiente forma SSreg_M0=SSTo−SSres_M0 y SSreg_M1=SSTo−SSres_M1 y su diferencia como SSΔ= SSreg_M0 - SSreg_M1 que sería igual a SSres_M0 - SSres_M1 . Ahora estos valores tenemos que ponerlos en función de los grados de libertad de cada modelo (también llamado la media de las sumas de cuadrados), entonces para el primer modelo MS_m1= SSRes_m1/df_m1 y para el modelo completo MSΔ= SSΔ/dfΔ. Finalmente obtenemos el estadístico F de la siguiente forma:

Figura 6

Estadístico F (Fratio) para el modelo 1

Creación de autor Alfonso Prado

Veamos un ejemplo para la comparación de los modelos de viviendas

La tabla tiene una fila para los grados de libertad residuales y la suma de cuadrados de cada modelo. Para todos los modelos, excepto el primero, también se proporciona la diferencia en los grados de libertad con respecto al modelo 1 y la suma de cuadrados obtenida. Esto solo tiene sentido estadístico si los modelos están anidados. Es habitual enumerar los modelos de menor a mayor, pero esto queda a criterio del analista.

Anova es en realidad una prueba de hipótesis cuyo H0 indica que los nuevos coeficientes (con respecto al modelo anterior) son 0 y H1 que al menos un coeficiente es distinto de 0.

Si los nuevos coeficientes son 0, entonces entendemos que las correspondientes variables no tienen valor de predicción y anova recomendará el modelo más sencillo. En principio, si los nuevos coeficientes mejoran el RSS, anova preferirá el modelo más complejo (mayor df), pero si no hay mejora anova seleccionará el modelo más simple.

Prueba de Significancia:

El ANOVA también puede comprobar la significancia de términos individuales o grupos de términos dentro de un modelo examinando el cambio en la desviación estándar al eliminarlos. Esto se realiza comparando el estadístico F calculado, con un valor crítico o valor p obtenido de la distribución F. Si el valor p es inferior al nivel de significancia (normalmente 0,05), se rechaza la hipótesis nula y las diferencias entre las medias de los grupos se consideran estadísticamente significativas.

prueba de razón de verosimilitud:

El análisis de la desviación estándar está estrechamente relacionado con la prueba de razón de likelihood, una prueba estadística que se utiliza para comparar el ajuste de dos modelos anidados.

Veamos un ejemplo

En clases anteriores habíamos revisado la técnica de validación cruzada como un mecanismo para prevenir el sobreajuste. Pero la validación cruzada sirve principalmente como una medida de la bondad del modelo mediante el entrenamiento del modelo con distintos subconjuntos de datos.

Para empezar, diremos que la validación cruzada es una técnica de re muestreo, esto quiere decir que partiendo de un dataset grande puedo generar múltiples dataset pequeños.

Para propósito de evitar el sobreajuste habíamos indicado la necesidad de tener 2 subconjuntos conocidos como data de entrenamiento y prueba. Pero para propósito de la evaluación de la bondad del modelo entrenaremos el modelo usando varios subconjuntos de datos y los modelos resultantes se evalúan posteriormente en otros subconjuntos, que no se utilizaron durante el entrenamiento. El promedio del rendimiento que alcanzan los modelos con estos subconjuntos es una estimación más precisa del rendimiento del modelo.

Existen muchas técnicas de re-muestreo, a continuación, algunas de las más comunes.

• Re muestreo simple: Es el método más sencillo de validación consiste en repartir aleatoriamente las observaciones disponibles en dos grupos, uno se emplea para entrenar al modelo y otro para evaluarlo. La estimación del error es altamente variable dependiendo de qué observaciones se incluyan como conjunto de entrenamiento y cuáles como conjunto de validación.
• Dejar uno afuera (LOOCV): Es un método iterativo que se inicia empleando como conjunto de entrenamiento todas las observaciones disponibles excepto una, que se excluye para emplearla como validación. Si se emplea una única observación para calcular el error, este varía mucho dependiendo de qué observación se haya seleccionado. El proceso requiere que el modelo sea reajustado y validado tantas veces como observaciones disponibles (n) lo que en algunos casos puede ser muy complicado.
• Re muestreo aleatorio k-fold Consiste en dividir los datos de forma aleatoria en k grupos de aproximadamente el mismo tamaño, k-1 grupos se emplean para entrenar el modelo y uno de los grupos se emplea como validación. Este proceso se repite k veces utilizando un grupo distinto como validación en cada iteración. El proceso genera k estimaciones del error cuyo promedio se emplea como estimación final. En la práctica valores de k entre 5 y 10 son recomendados, lo cual no implica un costo computacional alto. La principal ventaja de K-fold CV es que consigue una estimación precisa del error de test gracias a un mejor balance entre bias y varianza

Realizar el re muestreo, entrenar el modelo y probar con otro conjunto parecería una ardua labor, pero R lo hace muy fácil a través del uso de la función cv.glm() del paquete boot. Esta función retorna una lista con algunos valores entre los cuales se halla el estadístico delta. Delta es un vector de longitud dos. El primer componente es el error de predicción normal y el segundo componente es un poco más difícil de entender. Recuerde que el modelo se creó con la data de entrenamiento menos 1 fold. Se requiere un ajuste por este motivo y eso es el segundo componente (error de cross validación ajustado por efecto del k-1).

cv.glm toma como argumento el valor k (cantidad de grupos en los que se deben dividir los datos kfold) si el valor es igual al número de observaciones asume que debe usar LOOCV.

En base a la métrica delta preferimos el modelo con el menor rango de validación cruzada.

Figura 7

Comparación de modelos usado validación cruzada