Los criterios usados para evaluar el ajuste de los modelos como el \(R^2\) y el error cuadrático medio \((\text{RMSE})\), tienen una tendencia a mejorar el modelo al agregar predictores, es decir, a través de estos criterios generalmente se escoge el modelo con más variables, lo que genera un ajuste con variables que pueden estar generando ruido.
Es por esto, que se requiere un criterio de calidad que tenga en cuenta el tamaño del modelo (cantidad de covariables), dado que es preferible un modelo pequeño que ajuste bien, en algunos casos se podría sacrificar una cantidad pequeña de “bondad de ajuste” para obtener un modelo más pequeño.
El logaritmo de la verosimilitud en un modelo de regresión se puede escribir como:
\[ \log L(\boldsymbol{\hat{\beta}}, \hat{\sigma}^2) = -\frac{n}{2}\log(2\pi) - \frac{n}{2}\log\left(\frac{\text{RSS}}{n}\right) - \frac{n}{2}, \]
donde \(\text{RSS} = \sum_{i=1}^n (y_i - \hat{y}_i) ^ 2\) y \(\boldsymbol{\hat{\beta}}\) y \(\hat{\sigma}^2\) son seleccionados por el método de máxima verosimilitud:
Luego, se puede definir \(\text{AIC}\) como
\[ \text{AIC} = -2 \log L(\boldsymbol{\hat{\beta}}, \hat{\sigma}^2) + 2p = n + n \log(2\pi) + n \log\left(\frac{\text{RSS}}{n}\right) + 2p, \]
El componente de penalidad del \(\text{AIC}\) es:
\[ 2p, \]
donde \(p\) es el número de parámetros \(\beta\) en el modelo, porque entre más grande sea \(p\), el \(\text{AIC}\) va aumentar pero en realidad se busca un valor pequeño \(\text{AIC}\). De esta forma, un modelo con bajo \(\text{AIC}\), tendrá un buen balance entre el ajuste y el número pequeño de parámetros.
El BIC o \(\text{BIC}\), es similar a \(\text{AIC}\), pero tiene una penalidad mayor. \(\text{BIC}\) también compensa entre un modelo que ajuste bien y el número de parámetros en el mismo, sin embargo, generalmente tiende a seleccionar modelos con menor número de parámetros que el \(\text{AIC}\). Para la selección, se escoge el modelo con menor \(\text{BIC}\).
\[ \text{BIC} = -2 \log L(\boldsymbol{\hat{\beta}}, \hat{\sigma}^2) + \log(n) p = n + n\log(2\pi) + n\log\left(\frac{\text{RSS}}{n}\right) + \log(n)p. \]
Note que en \(\text{AIC}\) la penalidad era
\[ 2p, \]
Mientras que \(\text{BIC}\), la penalidad es
\[ \log(n) p. \]
Entonces, para cualquier set de datos donde \(log(n) > 2\) la penalidad del \(\text{BIC}\) será mayor que \(\text{AIC}\), por tanto \(\text{BIC}\) seleccionará un modelo más pequeño.
Recordando que,
\[ R^2 = 1 - \frac{\text{SSE}}{\text{SST}} = 1 - \frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2}. \]
Se define ahora:
\[ R_a^2 = 1 - \frac{\text{SSE}/(n-p)}{\text{SST}/(n-1)} = 1 - \left( \frac{n-1}{n-p} \right)(1-R^2) \]
Denominado \(R^2\) ajustado.
Debido a que \(R^2\) no se puede volver más pequeño a medida que aumentan el número de predictores, \(R^2\) ajustado penaliza efectivamente por predictores adicionales y puede disminuir si se aumentan, como \(R^2\), entre mayor es mejor.
El proceso de selección de variable involucra los criterios de calidad del modelo además del proceso de búsqueda y selección. Aqui es importante tener en cuenta que, un modelo de regresión permite explicar y predecir, en la explicación el objetivo es poder establecer relaciones y la magnitud de estas, mientras en la predicción se buscan modelos que puedan entender el comportamiento del set de datos y replicarlo ante nuevas observaciones. Los métodos de selección que aqui se presentan principalmente aplican para objetivos de predicción.
Teniendo en cuenta los datos seatpos de la librería
faraway, los cuales muestran información sobre la posición
del cinturón de seguridad (variable dependiente hipcenter: distancia
horizontal del punto medio de las caderas desde un punto fijo en el
coche en mm) y variables relacionadas con esto como la edad, peso, talla
con zapatos, talla sin zapatos, altura de la persona sentada, longitud
del brazo, longitud de la pierna, etc. en total son 8 predictores
(variables independientes).
library(faraway)
hipcenter_mod = lm(hipcenter ~ ., data = seatpos)
coef(hipcenter_mod)## (Intercept) Age Weight HtShoes Ht Seated
## 436.43212823 0.77571620 0.02631308 -2.69240774 0.60134458 0.53375170
## Arm Thigh Leg
## -1.32806864 -1.14311888 -6.43904627Este método inicia con el modelo ajustado por todos los predictores
posibles, luego considera la eliminación de cada uno de acuerdo con una
métrica. Usando la función step() en R la cual
por defecto utiliza \(\text{AIC}\).
hipcenter_mod_back_aic = step(hipcenter_mod, direction = "backward")## Start: AIC=283.62
## hipcenter ~ Age + Weight + HtShoes + Ht + Seated + Arm + Thigh +
## Leg
##
## Df Sum of Sq RSS AIC
## - Ht 1 5.01 41267 281.63
## - Weight 1 8.99 41271 281.63
## - Seated 1 28.64 41290 281.65
## - HtShoes 1 108.43 41370 281.72
## - Arm 1 164.97 41427 281.78
## - Thigh 1 262.76 41525 281.87
## <none> 41262 283.62
## - Age 1 2632.12 43894 283.97
## - Leg 1 2654.85 43917 283.99
##
## Step: AIC=281.63
## hipcenter ~ Age + Weight + HtShoes + Seated + Arm + Thigh + Leg
##
## Df Sum of Sq RSS AIC
## - Weight 1 11.10 41278 279.64
## - Seated 1 30.52 41297 279.66
## - Arm 1 160.50 41427 279.78
## - Thigh 1 269.08 41536 279.88
## - HtShoes 1 971.84 42239 280.51
## <none> 41267 281.63
## - Leg 1 2664.65 43931 282.01
## - Age 1 2808.52 44075 282.13
##
## Step: AIC=279.64
## hipcenter ~ Age + HtShoes + Seated + Arm + Thigh + Leg
##
## Df Sum of Sq RSS AIC
## - Seated 1 35.10 41313 277.67
## - Arm 1 156.47 41434 277.78
## - Thigh 1 285.16 41563 277.90
## - HtShoes 1 975.48 42253 278.53
## <none> 41278 279.64
## - Leg 1 2661.39 43939 280.01
## - Age 1 3011.86 44290 280.31
##
## Step: AIC=277.67
## hipcenter ~ Age + HtShoes + Arm + Thigh + Leg
##
## Df Sum of Sq RSS AIC
## - Arm 1 172.02 41485 275.83
## - Thigh 1 344.61 41658 275.99
## - HtShoes 1 1853.43 43166 277.34
## <none> 41313 277.67
## - Leg 1 2871.07 44184 278.22
## - Age 1 2976.77 44290 278.31
##
## Step: AIC=275.83
## hipcenter ~ Age + HtShoes + Thigh + Leg
##
## Df Sum of Sq RSS AIC
## - Thigh 1 472.8 41958 274.26
## <none> 41485 275.83
## - HtShoes 1 2340.7 43826 275.92
## - Age 1 3501.0 44986 276.91
## - Leg 1 3591.7 45077 276.98
##
## Step: AIC=274.26
## hipcenter ~ Age + HtShoes + Leg
##
## Df Sum of Sq RSS AIC
## <none> 41958 274.26
## - Age 1 3108.8 45067 274.98
## - Leg 1 3476.3 45434 275.28
## - HtShoes 1 4218.6 46176 275.90De acuerdo con la salida anterior, se observa que se inicia con el
modelo
hipcenter ~ Age + Weight + HtShoes + Ht + Seated + Arm + Thigh + Leg
con una AIC 283.62. Luego, se van eliminando cada predictor
hasta que para, o hasta que llega al modelo nulo (sin covariables).
El signo- al inicio de cada fila indica que se considera
eliminar esa variable en el modelo. La fila con
<none> indica el modelo actual. Se debe notar que
esta fila tiene el \(\text{RSS}\) más
pequeño cuando en realidad es el modelo con más predictores. Cada fila
por encima de <none> tiene un \(\text{AIC}\) más pequeño, se inicia con el
primero, la variable Ht que da el \(\text{AIC}\) más bajo. Se remueve la
variable y se continua con el proceso.
Se continúa con el proceso hasta llegar al modelo
hipcenter ~ Age + HtShoes + Leg. En este paso, la fila
<none> está de primera lo que muestra que remover
cualquier variable adicional no mejorará \(\text{AIC}\), y este es finalmente guardado
en hipcenter_mod_back_aic.
coef(hipcenter_mod_back_aic)## (Intercept) Age HtShoes Leg
## 456.2136538 0.5998327 -2.3022555 -6.8297461Se puede también usar el criterio de \(\text{BIC}\) con la función
step(), especificando k = log(n), donde
n guarda el número de observaciones en la base de
datos.
n = length(resid(hipcenter_mod))
hipcenter_mod_back_bic = step(hipcenter_mod, direction = "backward", k = log(n))## Start: AIC=298.36
## hipcenter ~ Age + Weight + HtShoes + Ht + Seated + Arm + Thigh +
## Leg
##
## Df Sum of Sq RSS AIC
## - Ht 1 5.01 41267 294.73
## - Weight 1 8.99 41271 294.73
## - Seated 1 28.64 41290 294.75
## - HtShoes 1 108.43 41370 294.82
## - Arm 1 164.97 41427 294.88
## - Thigh 1 262.76 41525 294.97
## - Age 1 2632.12 43894 297.07
## - Leg 1 2654.85 43917 297.09
## <none> 41262 298.36
##
## Step: AIC=294.73
## hipcenter ~ Age + Weight + HtShoes + Seated + Arm + Thigh + Leg
##
## Df Sum of Sq RSS AIC
## - Weight 1 11.10 41278 291.10
## - Seated 1 30.52 41297 291.12
## - Arm 1 160.50 41427 291.24
## - Thigh 1 269.08 41536 291.34
## - HtShoes 1 971.84 42239 291.98
## - Leg 1 2664.65 43931 293.47
## - Age 1 2808.52 44075 293.59
## <none> 41267 294.73
##
## Step: AIC=291.1
## hipcenter ~ Age + HtShoes + Seated + Arm + Thigh + Leg
##
## Df Sum of Sq RSS AIC
## - Seated 1 35.10 41313 287.50
## - Arm 1 156.47 41434 287.61
## - Thigh 1 285.16 41563 287.73
## - HtShoes 1 975.48 42253 288.35
## - Leg 1 2661.39 43939 289.84
## - Age 1 3011.86 44290 290.14
## <none> 41278 291.10
##
## Step: AIC=287.5
## hipcenter ~ Age + HtShoes + Arm + Thigh + Leg
##
## Df Sum of Sq RSS AIC
## - Arm 1 172.02 41485 284.02
## - Thigh 1 344.61 41658 284.18
## - HtShoes 1 1853.43 43166 285.53
## - Leg 1 2871.07 44184 286.41
## - Age 1 2976.77 44290 286.50
## <none> 41313 287.50
##
## Step: AIC=284.02
## hipcenter ~ Age + HtShoes + Thigh + Leg
##
## Df Sum of Sq RSS AIC
## - Thigh 1 472.8 41958 280.81
## - HtShoes 1 2340.7 43826 282.46
## - Age 1 3501.0 44986 283.46
## - Leg 1 3591.7 45077 283.54
## <none> 41485 284.02
##
## Step: AIC=280.81
## hipcenter ~ Age + HtShoes + Leg
##
## Df Sum of Sq RSS AIC
## - Age 1 3108.8 45067 279.89
## - Leg 1 3476.3 45434 280.20
## <none> 41958 280.81
## - HtShoes 1 4218.6 46176 280.81
##
## Step: AIC=279.89
## hipcenter ~ HtShoes + Leg
##
## Df Sum of Sq RSS AIC
## - Leg 1 3038.8 48105 278.73
## <none> 45067 279.89
## - HtShoes 1 5004.4 50071 280.25
##
## Step: AIC=278.73
## hipcenter ~ HtShoes
##
## Df Sum of Sq RSS AIC
## <none> 48105 278.73
## - HtShoes 1 83534 131639 313.35El procedimiento es el mismo, excepto que se busca mejorar \(\text{BIC}\), el cual R
etiqueta como \(\text{AIC}\) en la
salida.
coef(hipcenter_mod_back_bic)## (Intercept) HtShoes
## 565.592659 -4.262091Se debe notar que este modelo tiene menos predictores, lo cual es lo
esperado. Se puede usar summary() para comparar los valores
ajustado de \(R^2\).
summary(hipcenter_mod)$adj.r.squared## [1] 0.6000855summary(hipcenter_mod_back_aic)$adj.r.squared## [1] 0.6531427summary(hipcenter_mod_back_bic)$adj.r.squared## [1] 0.6244149En este proceso se inicia con el modelo sin predictores
hipcenter ~ 1, luego en cada paso se adiciona predictores
hasta que se encuentra un buen modelo o alcanzan
hipcenter ~ Age + Weight + HtShoes + Ht + Seated + Arm + Thigh + Leg.
hipcenter_mod_start = lm(hipcenter ~ 1, data = seatpos)
hipcenter_mod_forw_aic = step(
hipcenter_mod_start,
scope = hipcenter ~ Age + Weight + HtShoes + Ht + Seated + Arm + Thigh + Leg,
direction = "forward")## Start: AIC=311.71
## hipcenter ~ 1
##
## Df Sum of Sq RSS AIC
## + Ht 1 84023 47616 275.07
## + HtShoes 1 83534 48105 275.45
## + Leg 1 81568 50071 276.98
## + Seated 1 70392 61247 284.63
## + Weight 1 53975 77664 293.66
## + Thigh 1 46010 85629 297.37
## + Arm 1 45065 86574 297.78
## <none> 131639 311.71
## + Age 1 5541 126098 312.07
##
## Step: AIC=275.07
## hipcenter ~ Ht
##
## Df Sum of Sq RSS AIC
## + Leg 1 2781.10 44835 274.78
## <none> 47616 275.07
## + Age 1 2353.51 45262 275.14
## + Weight 1 195.86 47420 276.91
## + Seated 1 101.56 47514 276.99
## + Arm 1 75.78 47540 277.01
## + HtShoes 1 25.76 47590 277.05
## + Thigh 1 4.63 47611 277.06
##
## Step: AIC=274.78
## hipcenter ~ Ht + Leg
##
## Df Sum of Sq RSS AIC
## + Age 1 2896.60 41938 274.24
## <none> 44835 274.78
## + Arm 1 522.72 44312 276.33
## + Weight 1 445.10 44390 276.40
## + HtShoes 1 34.11 44801 276.75
## + Thigh 1 32.96 44802 276.75
## + Seated 1 1.12 44834 276.78
##
## Step: AIC=274.24
## hipcenter ~ Ht + Leg + Age
##
## Df Sum of Sq RSS AIC
## <none> 41938 274.24
## + Thigh 1 372.71 41565 275.90
## + Arm 1 257.09 41681 276.01
## + Seated 1 121.26 41817 276.13
## + Weight 1 46.83 41891 276.20
## + HtShoes 1 13.38 41925 276.23En este caso, las filas con + indican adición de
predictores al modelo.
hipcenter_mod_forw_bic = step(
hipcenter_mod_start,
scope = hipcenter ~ Age + Weight + HtShoes + Ht + Seated + Arm + Thigh + Leg,
direction = "forward", k = log(n))## Start: AIC=313.35
## hipcenter ~ 1
##
## Df Sum of Sq RSS AIC
## + Ht 1 84023 47616 278.34
## + HtShoes 1 83534 48105 278.73
## + Leg 1 81568 50071 280.25
## + Seated 1 70392 61247 287.91
## + Weight 1 53975 77664 296.93
## + Thigh 1 46010 85629 300.64
## + Arm 1 45065 86574 301.06
## <none> 131639 313.35
## + Age 1 5541 126098 315.35
##
## Step: AIC=278.34
## hipcenter ~ Ht
##
## Df Sum of Sq RSS AIC
## <none> 47616 278.34
## + Leg 1 2781.10 44835 279.69
## + Age 1 2353.51 45262 280.05
## + Weight 1 195.86 47420 281.82
## + Seated 1 101.56 47514 281.90
## + Arm 1 75.78 47540 281.92
## + HtShoes 1 25.76 47590 281.96
## + Thigh 1 4.63 47611 281.98Al comparar el \(R^2\) ajustado:
summary(hipcenter_mod)$adj.r.squared## [1] 0.6000855summary(hipcenter_mod_forw_aic)$adj.r.squared## [1] 0.6533055summary(hipcenter_mod_forw_bic)$adj.r.squared## [1] 0.6282374Los métodos descritos anteriormente son útiles sin embargo dejan por
fuera combinaciones de modelos y por tanto, se podría perder el mejor
modelo posible. Con un número alto de predictores revisar cada posible
modelo puede consumir mucho tiempo, sin embargo, con un número razonalbe
de predictores no es tan complejo usando la función
regsubsets() de la librería leaps.
library(leaps)
all_hipcenter_mod = summary(regsubsets(hipcenter ~ ., data = seatpos))Dentro de la función regsubsets() se específica el
modelo hipcenter ~ .. Este corresponde al modelo mas grande
usando todos los predictores y R revisará todos los
posibles subsets. Tiene regsubsets() tiene un argumento de
nvmax para ajustar el máximo número de subsets que por
defecto es \(8\).
all_hipcenter_mod$which## (Intercept) Age Weight HtShoes Ht Seated Arm Thigh Leg
## 1 TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## 2 TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE
## 3 TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE TRUE
## 4 TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE TRUE
## 5 TRUE TRUE FALSE TRUE FALSE FALSE TRUE TRUE TRUE
## 6 TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE
## 7 TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE
## 8 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEUsando $which da el mejor modelo, teniendo en cuenta
\(\text{RSS}\), para un modelo con cada
tamaño posible, en este caso teniendo en cuenta de 1 a 8 predictores.
Por ejemplo, el mejor modelo para (\(p =
5\)) usaría Age, HtShoes,
Thigh, and Leg.
Se puede obtener el \(\text{RSS}\)
para cada modelo usando $rss. Se debe notar que este
disminuye dado que inician desde el más pequeño al más grande.
all_hipcenter_mod$rss## [1] 47615.79 44834.69 41938.09 41485.01 41313.00 41277.90 41266.80 41261.78Se puede obtener el \(R^2\)
all_hipcenter_mod$adjr2## [1] 0.6282374 0.6399496 0.6533055 0.6466586 0.6371276 0.6257403 0.6133690
## [8] 0.6000855Se puede extraer el modelo con el \(R^2\) más alto:
(best_r2_ind = which.max(all_hipcenter_mod$adjr2))## [1] 3Se pueden extraer los predictores de ese modelo:
all_hipcenter_mod$which[best_r2_ind, ]## (Intercept) Age Weight HtShoes Ht Seated
## TRUE TRUE FALSE FALSE TRUE FALSE
## Arm Thigh Leg
## FALSE FALSE TRUEPara calcular el \(\text{AIC}\) y \(\text{BIC}\) para cada uno de los modelos con el mejor \(\text{RSS}\), se requieren a \(n\) y \(p\) para el modelo más grande.
p = length(coef(hipcenter_mod))
n = length(resid(hipcenter_mod))Se utilizará la siguiente expresión del \(\text{AIC}\):
\[ \text{AIC} = n\log\left(\frac{\text{RSS}}{n}\right) + 2p. \]
Como se tiene \(\text{RSS}\) para cada modelo, se puede calcular.
hipcenter_mod_aic = n * log(all_hipcenter_mod$rss / n) + 2 * (2:p)Se pueden extraer los predictores del modelo con el mejor \(\text{AIC}\).
best_aic_ind = which.min(hipcenter_mod_aic)
all_hipcenter_mod$which[best_aic_ind,]## (Intercept) Age Weight HtShoes Ht Seated
## TRUE TRUE FALSE FALSE TRUE FALSE
## Arm Thigh Leg
## FALSE FALSE TRUEAhora se ajusta este modelo y se compara con los anteriores usando \(\text{AIC}\):
hipcenter_mod_best_aic = lm(hipcenter ~ Age + Ht + Leg, data = seatpos)extractAIC(hipcenter_mod_best_aic)## [1] 4.0000 274.2418extractAIC(hipcenter_mod_back_aic)## [1] 4.0000 274.2597extractAIC(hipcenter_mod_forw_aic)## [1] 4.0000 274.2418Ejemplo tomado de: https://book.stat420.org/model-diagnostics.html#checking-assumptions