15 Todas las regresiones posibles
En este capítulo se muestra como elegir el mejor modelo de regresión cuando es posible ajustar o entrenar todos los posibles modelos relacionados a un problema.
Si un modelo tiene \(k\) variables cuantitativas, se pueden construir \(2^k\) modelos diferentes con subconjuntos de las variables originales. A continuación un ejemplo de cómo crece el número de todas las regresiones posibles en función del número \(k\) de covariables.
## k num_de_regresiones
## [1,] 2 4
## [2,] 4 16
## [3,] 8 256
## [4,] 16 65536
## [5,] 32 4294967296Criterios para elegir modelos
A continuación una lista de los posibles indicadores y el ideal de cada uno de ellos.
- \(R^2\): coeficiente de determinación, entre más grande mejor.
- \(R^2_A\): coeficiente de determinación ajustado, entre más grande mejor.
- \(C_p\) de Mallows: el mejor modelo es aquél para el cual \(C_p\) es el más pequeño posible.
- \(AIC\): Akaike Information Criterium, entre más pequeño mejor.
Paquete olsrr
El paquete olsrr creado por Hebbali (2024) contiene varias funciones útiles para modelación, en particular se destacan dos funciones.
ols_step_all_possible(model)
ols_step_best_subset(model,
max_order = NULL,
include = NULL,
exclude = NULL,
metric = c("rsquare", "adjr", "predrsq",
"cp", "aic", "sbic", "sbc",
"msep", "fpe", "apc", "hsp"))Ejemplo
En este ejemplo vamos a usar la base de datos mtcars para crear todas las regresiones posibles para explicar mpg en función de las covariables disp, hp, wt, qsec, a continuación una parte de la base de datos.
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1Solución
Lo primero que se debe hacer es ajustar el modelo saturado (o modelo full) para luego crear todas las regresiones posibles que incluyan las covariables disp, hp, wt, qsec.
model <- lm(mpg ~ disp + hp + wt + qsec, data=mtcars)
library(olsrr)
res <- ols_step_all_possible(model)
res # Este objeto contiene los resultados## Index N Predictors R-Square Adj. R-Square Mallow's Cp
## 3 1 1 wt 0.7528328 0.7445939 0.70869536
## 1 2 1 disp 0.7183433 0.7089548 0.67512054
## 2 3 1 hp 0.6024373 0.5891853 0.50969578
## 4 4 1 qsec 0.1752963 0.1478062 0.07541973
## 8 5 2 hp wt 0.8267855 0.8148396 0.78108710
## 10 6 2 wt qsec 0.8264161 0.8144448 0.77856272
## 6 7 2 disp wt 0.7809306 0.7658223 0.72532105
## 5 8 2 disp hp 0.7482402 0.7308774 0.69454380
## 7 9 2 disp qsec 0.7215598 0.7023571 0.66395284
## 9 10 2 hp qsec 0.6368769 0.6118339 0.52014395
## 14 11 3 hp wt qsec 0.8347678 0.8170643 0.78199548
## 11 12 3 disp hp wt 0.8268361 0.8082829 0.76789526
## 13 13 3 disp wt qsec 0.8264170 0.8078189 0.76988533
## 12 14 3 disp hp qsec 0.7541953 0.7278591 0.68301440
## 15 15 4 disp hp wt qsec 0.8351443 0.8107212 0.77102968Los resultados anteriores se pueden ver de forma gráfica así:
Ejemplo
En este ejemplo vamos a usar la base de datos mtcars para identificar el mejor modelo que explique mpg en función de las covariables disp, hp, wt, qsec. La restricción que se impone es que el modelo TIENE que incluir la variable hp y la métrica a usar será el \(R^2_{Adj}\).
Solución
model <- lm(mpg ~ disp + hp + wt + qsec, data=mtcars)
library(olsrr)
res <- ols_step_best_subset(model, include="hp", metric="adjr")
res # Este objeto contiene los resultados## Best Subsets Regression
## ------------------------------
## Model Index Predictors
## ------------------------------
## 1 hp
## 2 hp wt
## 3 hp wt qsec
## 4 disp hp wt qsec
## ------------------------------
##
## Subsets Regression Summary
## ----------------------------------------------------------------------------------------------------------------------------------
## Adj. Pred
## Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
## ----------------------------------------------------------------------------------------------------------------------------------
## 1 0.6024 0.5892 0.5097 37.1126 181.2386 87.8752 185.6358 477.5836 15.8551 0.5146 0.4506
## 2 0.8268 0.8148 0.7811 2.3690 156.6523 66.5755 162.5153 215.5104 7.3563 0.2402 0.2091
## 3 0.8348 0.8171 0.782 3.0617 157.1426 67.7238 164.4713 213.1929 7.4756 0.2461 0.2124
## 4 0.8351 0.8107 0.771 5.0000 159.0696 70.0408 167.8640 220.8882 7.9497 0.2644 0.2259
## ----------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria
## SBIC: Sawa's Bayesian Information Criteria
## SBC: Schwarz Bayesian Criteria
## MSEP: Estimated error of prediction, assuming multivariate normality
## FPE: Final Prediction Error
## HSP: Hocking's Sp
## APC: Amemiya Prediction CriteriaLos resultados anteriores se pueden ver de forma gráfica así:
Se le recomienda al lector explorar las viñetas del paquete olsrr para conocer todas las bondades del paquete.
Ejemplo
Volver a realizar el ejemplo anterior pero usando la función myAllRegTable del curso de Estadística II.
Solución
Para accerder a la función myAllRegTable se corre el siguiente código y se carga el paquete leaps creado por Lumley (2020).
library(leaps)
source("https://raw.githubusercontent.com/fhernanb/Trabajos_Est_2/main/Trabajo%201/funciones.R")Luego se ajusta el modelo y luego si se aplica la función myAllRegTable así:
## k R_sq adj_R_sq SSE Cp Variables_in_model
## 1 1 0.753 0.745 278.322 12.481 wt
## 2 1 0.718 0.709 317.159 18.130 disp
## 3 1 0.602 0.589 447.674 37.113 hp
## 4 1 0.175 0.148 928.655 107.070 qsec
## 5 2 0.827 0.815 195.048 2.369 hp wt
## 6 2 0.826 0.814 195.464 2.429 wt qsec
## 7 2 0.781 0.766 246.683 9.879 disp wt
## 8 2 0.748 0.731 283.493 15.233 disp hp
## 9 2 0.722 0.702 313.537 19.603 disp qsec
## 10 2 0.637 0.612 408.894 33.472 hp qsec
## 11 3 0.835 0.817 186.059 3.062 hp wt qsec
## 12 3 0.827 0.808 194.991 4.361 disp hp wt
## 13 3 0.826 0.808 195.463 4.429 disp wt qsec
## 14 3 0.754 0.728 276.788 16.258 disp hp qsec
## 15 4 0.835 0.811 185.635 5.000 disp hp wt qsec