Regresión Lineal Múltiple
Ángel Rmz.
Ejemplo 3
Se tiene información sobre el Rendimiento de la gasolina para 32 automóviles.
library("xlsx")
datos=read.xlsx("rendimiento_autos.xlsx",1)
datos=datos[,c(-1)]
head(datos)## y x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11
## 1 18.90 350 165 260 8.00 2.56 4 3 200.3 69.9 3910 A
## 2 17.00 350 170 275 8.50 2.56 4 3 199.6 72.9 2860 A
## 3 20.00 250 105 185 8.25 2.73 1 3 196.7 72.2 3510 A
## 4 18.25 351 143 255 8.00 3.00 2 3 199.9 74.0 3890 A
## 5 20.07 225 95 170 8.40 2.76 1 3 194.1 71.8 3365 M
## 6 11.20 440 215 330 8.20 2.88 4 3 184.5 69.0 4215 A\(y:\) Millas/galón
\(x_1:\) Cilindrada (pulgadas cúbicas)
\(x_2:\) Potencia (Hp)
\(x_3:\) Par de torsión (pies-lb)
\(x_4:\) Relación de compresión
\(x_5:\) Relación de eje trasero
\(x_6:\) Carburador (gargantas)
\(x_7:\) Número de velocidades en la transmisión
\(x_8:\) Longitud total (pulgadas)
\(x_9:\) Ancho (pulgadas)
\(x_{10}:\) Peso (lb)
\(x_{11}:\) Tipo de transmisión (A = automática, M = manual)
Ajuste un modelo de RLM
Solución
Ajustaremos un modelo de RLM considerando lo siguiente:
Emplearemos el método forward.
El critero de selección del mejor modelo estará basado en la Suma de Cuadrados Residual (RSS), el \(BIC\) y \(C_p\) de Mallow.
El método forward (selección hacia adelante) comienza con la hipótesis que no hay regresores en el modelo además de la ordenada al origen. Se trata de determinar un subconjunto óptimo insertando regresores, uno por uno, en el modelo. La idea es, aquella variable que mejore en mayor medida el modelo se selecciona. A continuación se intenta incrementar el modelo probando a introducir una a una las variables restantes. Si introduciendo alguna de ellas mejora, también se selecciona.
La Suma de Cuadrados Residual (Suma de Cuadrados del Error) se puede emplear para comparar modelos que tienen el mismo número de variables regresoras.
El BIC (Criterio de Información Bayesiana) es una medida de la calidad de un modelo estadístico, considera la bondad de ajuste del modelo y su complejidad.
\[ BIC=-2log(L)+log(n)*n_{parametros} \]
El mejor modelo es aquel que tenga un \(BIC\) pequeño.
EL \(C_p\) de Mallow es una estadística que nos permite comparar modelos en términos de sesgo.
\[ C_p=\frac{RSS(p)}{\hat{\sigma}^2}-n+2p \]
Los modelos con poco sesgo tendrán valores de \(C_p\) que caen cerca de la recta \(C_p = p\). En general, diremos que el mejor modelo es aquel que tenga un \(C_p\) pequeño.
Paso 0: Sobre la Base de Datos y Estadística Descriptiva
Siempre es importante ver que la base de datos esté ad hoc para poder empezar a hacer un análisis. Es importante ver si hay datos faltantes o nulos,etcétera.
Después de establecer que la base está ad hoc resulta de gran ayuda hacer estadística descriptiva. Es importante observar si hay variables numéricas, categóricas, etc.
Primer Paso: Relación entre Variables (Indicios)
Podemos hacer gráficas de dispersión para ver el comportamiento de los datos.
library("tidyverse")
library("GGally")
attach(datos)
ggpairs(datos, upper = list(continuous = "smooth"),
diag = list(continuous = "bar"),lower=list(continuous="cor"))
Segundo Paso: Generación de Modelos y Elección del Mejor Modelo
Antes de generar los modelos y elegir el mejor, vamos a hacer el ajuste considerando todas las variables regresoras, esto con la finalidad de ver lo que está pasando y mencionar algunas cosas importantes.
#generamos el modelo con todas las variables regresoras
modelo_completo=lm(y~.,data=datos)
#hagamos un análisis del modelo completo (nadamás para ver que está pasando)
summary(modelo_completo)##
## Call:
## lm(formula = y ~ ., data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8330 -2.1977 -0.3141 1.3560 5.8153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.674837 24.353685 1.342 0.195
## x1 -0.041375 0.032731 -1.264 0.221
## x2 -0.007591 0.058346 -0.130 0.898
## x3 0.020348 0.026026 0.782 0.443
## x4 0.620928 2.437856 0.255 0.802
## x5 4.397581 3.214251 1.368 0.186
## x6 0.030013 1.279192 0.023 0.982
## x7 -2.544539 3.016119 -0.844 0.409
## x8 0.138900 0.097132 1.430 0.168
## x9 -0.447410 0.304626 -1.469 0.157
## x10 -0.002515 0.002788 -0.902 0.378
## x11M -0.087429 3.092854 -0.028 0.978
##
## Residual standard error: 3.342 on 20 degrees of freedom
## Multiple R-squared: 0.8195, Adjusted R-squared: 0.7203
## F-statistic: 8.258 on 11 and 20 DF, p-value: 2.889e-05Nota importante: La función lm() acepta variables categóricas y las codifica en automático,esto es una gran ventaja de R (en particular de la libreria stats) respecto a otros programas como Python en donde es necesario hacer a mano las variables dummys (one-hot encoding) para cada nivel en la variable categórica.
Ahora bien, en este caso la variable categórica \(x_{11}\) tiene dos niveles (dos categorías), entonces lo que sucede es que \(x_{11}\) se codifica a 0 y 1, donde el valor de 1 representa cuando \(x_{11}\)=“Manual” y el valor de 0 cuando \(x_{11}\)=Automático. En general, cuando una variable tiene \(c\) categorías (niveles) la función lm() crea \(c-1\) regresores nuevos para poder ajustar el modelo; en este caso es importante identificar cómo se está codificando las variables.
Ahora si, procedemos a generar los modelos y elegir el mejor. Emplearemos el método forward y aplicaremos diversos criterios de selección.
library("leaps")
#la función regsubsets() nos permite generar los mejores modelos de acuerdo al número de regresores y tomando como criterio RSS, además nos permite elegir el método que queramos:forward o backward.
mejores_modelos <- regsubsets(y~., data = datos,method = "forward",nvmax = 11)
mejores=summary(mejores_modelos)
mejores## Subset selection object
## Call: regsubsets.formula(y ~ ., data = datos, method = "forward", nvmax = 11)
## 11 Variables (and intercept)
## Forced in Forced out
## x1 FALSE FALSE
## x2 FALSE FALSE
## x3 FALSE FALSE
## x4 FALSE FALSE
## x5 FALSE FALSE
## x6 FALSE FALSE
## x7 FALSE FALSE
## x8 FALSE FALSE
## x9 FALSE FALSE
## x10 FALSE FALSE
## x11M FALSE FALSE
## 1 subsets of each size up to 11
## Selection Algorithm: forward
## x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11M
## 1 ( 1 ) "*" " " " " " " " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " "*" " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " " " "*" " " " " "*" " " " "
## 4 ( 1 ) "*" " " " " " " " " "*" " " "*" "*" " " " "
## 5 ( 1 ) "*" " " " " " " "*" "*" " " "*" "*" " " " "
## 6 ( 1 ) "*" " " " " " " "*" "*" " " "*" "*" "*" " "
## 7 ( 1 ) "*" " " " " " " "*" "*" "*" "*" "*" "*" " "
## 8 ( 1 ) "*" " " "*" " " "*" "*" "*" "*" "*" "*" " "
## 9 ( 1 ) "*" " " "*" "*" "*" "*" "*" "*" "*" "*" " "
## 10 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" " "
## 11 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"La función regsubsets() genera los mejores modelos de acuerdo al número de regresores y tomando como criterio de selección RSS, como un segundo filtro nosotros podemos elegir lo mejor de lo mejor de acuerdo al \(BIC\), \(C_p\), \(R^2\) ajustado.
El mejor modelo de acuerdo al BIC
Si tomamos como criterio el \(BIC\) entonces tenemos lo siguiente:
library("tidyverse")
ggplot(data = data.frame(n_predictores = 1:11, BIC= mejores$bic),
aes(x = n_predictores, y = BIC)) +
geom_line() +
geom_point(size=3) +
geom_point(aes(x = n_predictores[which.min(BIC)],
y = BIC[which.min(BIC)]),
colour = "red", size = 4) +
scale_x_continuous(breaks = c(0:19))+
theme_bw() +
labs(title = "Método Forward bajo el criterio BIC",
x = "Número Predictores")
El mejor modelo de acuerdo al Cp
Si tomamos como criterio el \(C_p\) de Mallow entonces tenemos lo siguiente:
library("tidyverse")
ggplot(data = data.frame(n_predictores = 1:11, Cp= mejores$cp),
aes(x = n_predictores, y = Cp)) +
geom_line() +
geom_point(size=3) +
geom_point(aes(x = n_predictores[which.min(Cp)],
y = Cp[which.min(Cp)]),
colour = "blue", size = 4) +
scale_x_continuous(breaks = c(0:19))+
theme_bw() +
labs(title = "Método Forward bajo el criterio Cp",
x = "Número Predictores")
El mejor modelo de acuerdo al Coeficiente de Determinación Ajustado
Si tomamos como criterio el \(R^2\) ajustado entonces tenemos lo siguiente:
library("tidyverse")
ggplot(data = data.frame(n_predictores = 1:11, R2_ajustado= mejores$adjr2),
aes(x = n_predictores, y = R2_ajustado)) +
geom_line() +
geom_point(size=3) +
geom_point(aes(x = n_predictores[which.max(R2_ajustado)],
y = R2_ajustado[which.max(R2_ajustado)]),
colour = "pink", size = 4) +
scale_x_continuous(breaks = c(0:19))+
theme_bw() +
labs(title = "Método Forward bajo el criterio R2 ajustado",
x = "Número Predictores")
De acuerdo a lo visto en las gráficas anteriores hay un top 2: el modelo de un sólo regresor (\(x_1\)) y el modelo con dos regresores (\(x_1\) y \(x_6\)).
Veamos primero el modelo que considera dos regresores (\(x_1\) y \(x_6\)):
modelo2=lm(y~x1+x6)
summary(modelo2)##
## Call:
## lm(formula = y ~ x1 + x6)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.0456 -1.6368 -0.3348 1.6503 6.2540
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.910041 1.540929 21.357 < 2e-16 ***
## x1 -0.053025 0.006145 -8.628 1.68e-09 ***
## x6 0.929500 0.670108 1.387 0.176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.021 on 29 degrees of freedom
## Multiple R-squared: 0.7862, Adjusted R-squared: 0.7714
## F-statistic: 53.31 on 2 and 29 DF, p-value: 1.934e-10Lo anterior sugiere que \(x_6\) no es significativa. Tomaremos el modelo de un sólo regresor (\(x_1\)).
modelo1=lm(y~x1)
summary(modelo1)##
## Call:
## lm(formula = y ~ x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7875 -1.9616 0.0206 1.7878 6.8182
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.727439 1.445559 23.33 < 2e-16 ***
## x1 -0.047428 0.004706 -10.08 3.82e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.067 on 30 degrees of freedom
## Multiple R-squared: 0.772, Adjusted R-squared: 0.7644
## F-statistic: 101.6 on 1 and 30 DF, p-value: 3.82e-11A través de la tabla ANOVA se tiene:
anova(modelo1)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 955.34 955.34 101.56 3.82e-11 ***
## Residuals 30 282.20 9.41
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Tercer Paso: Validación de Supuestos
Ya que hemos obtenido el mejor modelo debemos hacer validación de supuestos.
Relación lineal entre regresores y la variable respuesta
library(ggplot2)
library(gridExtra)
plot1 <- ggplot(data = datos, aes(x1,y)) +
geom_point() + geom_smooth(color = "firebrick") +
theme_bw()
plot1
Todos los supuestos sobre los Residuales
library("tidyverse")
residuales=data.frame(ajustados=modelo1$fitted.values,residuales_estand=rstandard(modelo1))
ggplot(data=residuales,aes(x=ajustados,y=residuales_estand))+
geom_point(color="darkblue",size=8)+geom_hline(yintercept = 0)+
geom_hline(yintercept=2,color="red")+
geom_hline(yintercept = -2,color="red")
También podemos hacer gráficas como QQPLOT e histogramas para ver gráficamente la normalidad de los residuales. Eso les toca a ustedes.
Pruebas Formales
library("nortest")
ad.test(rstandard(modelo1))##
## Anderson-Darling normality test
##
## data: rstandard(modelo1)
## A = 0.2048, p-value = 0.861library("lmtest")
bptest(modelo1)##
## studentized Breusch-Pagan test
##
## data: modelo1
## BP = 5.2642, df = 1, p-value = 0.02177library("lmtest")
dwtest(modelo1)##
## Durbin-Watson test
##
## data: modelo1
## DW = 1.6676, p-value = 0.1634
## alternative hypothesis: true autocorrelation is greater than 0
Multicolinealidad
Dado que llegamos a un modelo con un sólo regresor, no tiene caso checar si hay problemas de Multicolinealidad.
Falta ver si hay datos influyentes y atípicos; eso ya les toca a ustedes
Cuarto Paso: Corrección del Modelo
Hemos Visto que los supuestos del modelo se cumplen,excepto el de la relación lineal entre \(x_1\) y \(y\) y la homocedasticidad.
A continuación se muestran transformaciones que se pueden aplicar para estabilizar la varianza.
Caption for the picture.
Aplicaremos el \(logaritmo\) a la variable respuesta \(y\) y veamos qué pasa.
library(ggplot2)
library(gridExtra)
plot1 <- ggplot(data = datos, aes(x1,log(y))) +
geom_point() + geom_smooth(color = "firebrick") +
theme_bw()
plot1
modelo1_transformado=lm(log(y)~x1)
summary(modelo1_transformado)##
## Call:
## lm(formula = log(y) ~ x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27360 -0.09393 0.01329 0.08760 0.26873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.5936771 0.0631215 56.93 < 2e-16 ***
## x1 -0.0022104 0.0002055 -10.76 8.16e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1339 on 30 degrees of freedom
## Multiple R-squared: 0.7941, Adjusted R-squared: 0.7872
## F-statistic: 115.7 on 1 and 30 DF, p-value: 8.159e-12
Pruebas Formales
library("nortest")
ad.test(rstandard(modelo1_transformado))##
## Anderson-Darling normality test
##
## data: rstandard(modelo1_transformado)
## A = 0.31839, p-value = 0.5202library("lmtest")
bptest(modelo1_transformado)##
## studentized Breusch-Pagan test
##
## data: modelo1_transformado
## BP = 2.0621e-06, df = 1, p-value = 0.9989library("lmtest")
dwtest(modelo1_transformado)##
## Durbin-Watson test
##
## data: modelo1_transformado
## DW = 1.4602, p-value = 0.05474
## alternative hypothesis: true autocorrelation is greater than 0
Conclusiones Respecto al Problema
¿Cuáles son las conclusiones respecto al problema?
El mejor modelo que elegimos fue:
\[ log(y)=3.5936-0.0022104*x1 \]
No obstante, el modelo que considera las variables \(x_1\) y \(x_6\) podría ser igualmente bueno.
Como Tarea Moral verifiquen que efectivamente el modelo que considera las variables \(x_1\) y \(x_6\) cumple con todos los supuestos de RLM.
Conclusiones Respecto a la Metodología Empleada
El método forward no garantiza que se seleccione el mejor modelo de entre todos los posibles ya que no se evalúan todas las posibles combinaciones. Sin embargo, suele llegar a modelos muy optimizados consiguiendo un buen rendimiento computacional.
El \(BIC\) es una buena medida para comparar modelos en términos de ajuste y complejidad. Agrega una mayor penalización que el \(AIC\).
El \(C_p\) de Mallow es una buena estadística para comparar modelos en términos de sesgo.
Conclusiones Generales
Para ajustar un modelo de RLM se recomienda emplear diversos métodos (backward, forward y todos los posibles modelos) en conjunto bajo varios criterios: \(AIC\), \(BIC\), \(R^2\) ajustada y \(C_p\).
Al emplear diversos métodos y diversos criterios de selección, no necesariamente se llega al mismo mejor modelo. Es decir, puede haber varios modelos igualmente buenos.
El Pilón
Podemos emplear el método forward bajo el criterio del \(AIC\) para el problema anterior:
library("MASS")
#tenemos que crear un modelo vacio(sin regresores) y el horizonte en donde pare el algoritmo
modelo_vacio=lm(y~1)
horizonte=formula(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11)
#la función stepAIC nos permite generar modelos paso a paso con selección hacia adelante y toma el criterio de AIC.
modelo_forward_AIC= stepAIC(modelo_vacio,direction = "forward",scope=horizonte)## Start: AIC=118.96
## y ~ 1
##
## Df Sum of Sq RSS AIC
## + x1 1 955.34 282.20 73.661
## + x10 1 871.62 365.93 81.974
## + x2 1 805.71 431.83 87.274
## + x3 1 795.96 441.59 87.988
## + x9 1 739.68 497.86 91.827
## + x8 1 704.72 532.82 93.998
## + x11 1 686.97 550.58 95.048
## + x7 1 645.12 592.42 97.391
## + x5 1 437.81 799.74 106.994
## + x6 1 293.50 944.04 112.302
## + x4 1 157.32 1080.22 116.614
## <none> 1237.54 118.965
##
## Step: AIC=73.66
## y ~ x1
##
## Df Sum of Sq RSS AIC
## + x6 1 17.5580 264.65 73.605
## <none> 282.20 73.661
## + x9 1 15.6268 266.58 73.838
## + x10 1 7.3255 274.88 74.819
## + x2 1 6.3666 275.84 74.930
## + x4 1 5.7410 276.46 75.003
## + x7 1 4.0104 278.19 75.203
## + x5 1 3.3683 278.83 75.276
## + x3 1 2.3190 279.88 75.397
## + x11 1 0.5238 281.68 75.601
## + x8 1 0.0067 282.20 75.660
##
## Step: AIC=73.6
## y ~ x1 + x6
##
## Df Sum of Sq RSS AIC
## <none> 264.65 73.605
## + x9 1 5.0100 259.64 74.993
## + x3 1 2.2897 262.36 75.327
## + x8 1 2.2628 262.38 75.330
## + x4 1 2.0773 262.57 75.353
## + x10 1 0.8435 263.80 75.503
## + x11 1 0.3979 264.25 75.557
## + x5 1 0.1324 264.51 75.589
## + x7 1 0.0158 264.63 75.603
## + x2 1 0.0040 264.64 75.605summary(modelo_forward_AIC)##
## Call:
## lm(formula = y ~ x1 + x6)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.0456 -1.6368 -0.3348 1.6503 6.2540
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.910041 1.540929 21.357 < 2e-16 ***
## x1 -0.053025 0.006145 -8.628 1.68e-09 ***
## x6 0.929500 0.670108 1.387 0.176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.021 on 29 degrees of freedom
## Multiple R-squared: 0.7862, Adjusted R-squared: 0.7714
## F-statistic: 53.31 on 2 and 29 DF, p-value: 1.934e-10¿Qué se puede concluir?
Existe otra metodología de generación de modelos, es una combinación del método forward y backward: método por segmentos. Ejemplificaremos dicho método con la función stepAIC:
library("MASS")
#tenemos que crear un modelo vacio y el horizonte en donde pare el algoritmo
modelo_vacio=lm(y~1)
horizonte=formula(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11)
#la función stepAIC nos permite generar modelos por segmentos y toma el criterio de AIC.
modelo_segmentos_AIC= stepAIC(modelo_vacio,direction = "both",scope=horizonte)## Start: AIC=118.96
## y ~ 1
##
## Df Sum of Sq RSS AIC
## + x1 1 955.34 282.20 73.661
## + x10 1 871.62 365.93 81.974
## + x2 1 805.71 431.83 87.274
## + x3 1 795.96 441.59 87.988
## + x9 1 739.68 497.86 91.827
## + x8 1 704.72 532.82 93.998
## + x11 1 686.97 550.58 95.048
## + x7 1 645.12 592.42 97.391
## + x5 1 437.81 799.74 106.994
## + x6 1 293.50 944.04 112.302
## + x4 1 157.32 1080.22 116.614
## <none> 1237.54 118.965
##
## Step: AIC=73.66
## y ~ x1
##
## Df Sum of Sq RSS AIC
## + x6 1 17.56 264.65 73.605
## <none> 282.20 73.661
## + x9 1 15.63 266.58 73.838
## + x10 1 7.33 274.88 74.819
## + x2 1 6.37 275.84 74.930
## + x4 1 5.74 276.46 75.003
## + x7 1 4.01 278.19 75.203
## + x5 1 3.37 278.84 75.276
## + x3 1 2.32 279.88 75.397
## + x11 1 0.52 281.68 75.601
## + x8 1 0.01 282.20 75.660
## - x1 1 955.34 1237.54 118.965
##
## Step: AIC=73.6
## y ~ x1 + x6
##
## Df Sum of Sq RSS AIC
## <none> 264.65 73.605
## - x6 1 17.56 282.20 73.661
## + x9 1 5.01 259.64 74.993
## + x3 1 2.29 262.36 75.327
## + x8 1 2.26 262.38 75.330
## + x4 1 2.08 262.57 75.353
## + x10 1 0.84 263.80 75.503
## + x11 1 0.40 264.25 75.557
## + x5 1 0.13 264.51 75.589
## + x7 1 0.02 264.63 75.603
## + x2 1 0.00 264.64 75.605
## - x1 1 679.39 944.04 112.302summary(modelo_segmentos_AIC)##
## Call:
## lm(formula = y ~ x1 + x6)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.0456 -1.6368 -0.3348 1.6503 6.2540
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.910041 1.540929 21.357 < 2e-16 ***
## x1 -0.053025 0.006145 -8.628 1.68e-09 ***
## x6 0.929500 0.670108 1.387 0.176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.021 on 29 degrees of freedom
## Multiple R-squared: 0.7862, Adjusted R-squared: 0.7714
## F-statistic: 53.31 on 2 and 29 DF, p-value: 1.934e-10¿Qué se puede concluir?
La paquetería olsrr contiene varias funciones útiles para la modelación, en particular se destaca la función ols_step_all_possible que se puede utilizar para obtener TODAS LAS REGRESIONES POSIBLES. Dicha función es computacionalmente costosa cuando el número de regresores es grande.
library("olsrr")
#creamos un modelo con ciertas variables regresoras para ejemplificar el uso de ols_step_all_possible()
modelo=lm(y~x1+x6+x8+x9,data=datos)
#corremos la función ols_step_all_possible()
resultados=ols_step_all_possible(modelo)
resultados## Index N Predictors R-Square Adj. R-Square Mallow's Cp
## 1 1 1 x1 0.7719647 0.7643635 2.711947
## 4 2 1 x9 0.5976999 0.5842899 26.182049
## 3 3 1 x8 0.5694537 0.5551022 29.986264
## 2 4 1 x6 0.2371662 0.2117384 74.738962
## 5 5 2 x1 x6 0.7861525 0.7714044 2.801122
## 7 6 2 x1 x9 0.7845920 0.7697362 3.011300
## 6 7 2 x1 x8 0.7719701 0.7562439 4.711219
## 9 8 2 x6 x9 0.6632559 0.6400321 19.352925
## 10 9 2 x8 x9 0.6205824 0.5944157 25.100213
## 8 10 2 x6 x8 0.6031480 0.5757789 27.448298
## 13 11 3 x1 x8 x9 0.7927686 0.7705653 3.910058
## 12 12 3 x1 x6 x9 0.7902009 0.7677224 4.255884
## 11 13 3 x1 x6 x8 0.7879810 0.7652646 4.554865
## 14 14 3 x6 x8 x9 0.6682582 0.6327144 20.679207
## 15 15 4 x1 x6 x8 x9 0.7995258 0.7698259 5.000000Podemos graficar lo obtenido en la función ols_step_all_possible() de la siguiente manera:
plot(resultados)
¿Qué se puede concluir?