En este ejercicio, generaremos datos simulados y luego utilizaremos los métodos de selección de predictores para seleccionar predictores. El objetivo es comprobar las resultados proporcionados por estos métodos en un ejemplo donde se conocen los predictores a incluir en el modelo lineal.

set.seed(1)

#200 observaciones de una normal con media 0 y desviación estandar 1
X<-rnorm(200,0,1)
#vector de errores (suponiendo normalidad)
e<-rnorm(200,0,0.1)

#definimos los predictores
x1<-X
x2<-X^2
x3<-X^3

#creamos el vector de respuestas
Y<-3+2*x1-7*x2+4*x3+e
#definimos las variables faltantes
x4 <- X^4
x5 <- X^5
x6 <- X^6
x7 <- X^7
x8 <- X^8

#creamos nuestro data frame
df <- data.frame(Y, X1 = x1, X2 = x2, X3 = x3, X4 = x4, X5 = x5, X6 = x6, X7 = x7, X8 = x8)

head(df)
##           Y         X1         X2           X3          X4            X5
## 1 -1.942471 -0.6264538 0.39244438 -0.245848275 0.154012589 -0.0964817733
## 2  3.324873  0.1836433 0.03372487  0.006193347 0.001137367  0.0002088698
## 3 -5.734519 -0.8356286 0.69827518 -0.583498718 0.487588224 -0.4074426711
## 4  4.582478  1.5952808 2.54492084  4.059863355 6.476622070 10.3320308510
## 5  2.813570  0.3295078 0.10857537  0.035776429 0.011788611  0.0038844391
## 6 -5.312603 -0.8204684 0.67316837 -0.552313364 0.453155653 -0.3717998868
##             X6            X7           X8
## 1 6.044137e-02 -3.786373e-02 2.371988e-02
## 2 3.835755e-05  7.044108e-06 1.293603e-06
## 3 3.404708e-01 -2.845071e-01 2.377423e-01
## 4 1.648249e+01  2.629420e+01 4.194663e+01
## 5 1.279953e-03  4.217544e-04 1.389714e-04
## 6 3.050501e-01 -2.502839e-01 2.053500e-01
#definimos nuestra muestra de entrenamiento y nuestra muestra test
entrenamiento <- Y[1:150]
test <- Y[151:200]

df_entrenamiento <- df[1:150, ]
df_test <- df[151:200, ]
library(ggplot2)
#Diagrama de dispersión muestra de entrenamiento: X1 frente a Y
ggplot(df_entrenamiento, aes(x = Y, y = X1)) + geom_point(color = "blue")+
  labs(title = "Diagrama de Dispersión X1 frente a Y", x = "Y", y = "X1")

Como era de esperarse, existe una clara correlación entre las variables X1 e Y. Notemos que sabiamos que Y = 3+2X1-7X12+4X13+e. Los errores no son lo suficientemente grandes para que esta relación no sea visible.

Selección de predictores

Mejor selección de predictores

#La mejor selección de predictores

install.packages("leaps")
## 
## The downloaded binary packages are in
##  /var/folders/jd/zn7_v9ls5n35srxc2v42_ymh0000gn/T//RtmpHjYvwi/downloaded_packages
library(leaps)

# Ajuste de modelo 
mejor_modelo <- regsubsets(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, data = df_entrenamiento, nbest = 1)

summary(mejor_modelo)
## Subset selection object
## Call: regsubsets.formula(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, 
##     data = df_entrenamiento, nbest = 1)
## 8 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
## 1 subsets of each size up to 8
## Selection Algorithm: exhaustive
##          X1  X2  X3  X4  X5  X6  X7  X8 
## 1  ( 1 ) " " " " "*" " " " " " " " " " "
## 2  ( 1 ) " " "*" "*" " " " " " " " " " "
## 3  ( 1 ) "*" "*" "*" " " " " " " " " " "
## 4  ( 1 ) "*" "*" "*" "*" " " " " " " " "
## 5  ( 1 ) "*" "*" "*" "*" "*" " " " " " "
## 6  ( 1 ) "*" "*" "*" " " "*" "*" "*" " "
## 7  ( 1 ) "*" "*" "*" "*" "*" " " "*" "*"
## 8  ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
#Con los resultados comparamos los r2 ajustado, AIC, y BIC de los modelos
summary(mejor_modelo)$adjr2 #R2 ajustado
## [1] 0.6300958 0.9942557 0.9999405 0.9999414 0.9999410 0.9999408 0.9999404
## [8] 0.9999400
summary(mejor_modelo)$cp #AIC
## [1] 9.118069e+05 1.392216e+04 2.640549e+00 1.609006e+00 3.592199e+00
## [6] 5.122250e+00 7.081923e+00 9.000000e+00
summary(mejor_modelo)$bic #BIC
## [1]  -140.1655  -760.9274 -1442.5080 -1440.6747 -1435.6819 -1431.1699 -1426.2021
## [8] -1421.2786
#Obtenemos el mejor modelo de acuerdo a cada uno de los criterios de información

nPredictoresR2<-which.max(summary(mejor_modelo)$adjr2)
nPredictoresAIC<-which.min(summary(mejor_modelo)$cp)
nPredictoresBIC<-which.min(summary(mejor_modelo)$bic)

nPredictoresR2
## [1] 4
nPredictoresAIC
## [1] 4
nPredictoresBIC
## [1] 3
#obtenemos los mejores modelos
M4 <- lm(Y ~ X1 + X2 + X3 + X4, data = df_entrenamiento)
M3 <- lm(Y ~ X1 + X2 + X3, data = df_entrenamiento)

coef(M4)
##  (Intercept)           X1           X2           X3           X4 
##  2.988660451  1.962508630 -6.959992964  4.015385247 -0.008652045
coef(M3)
## (Intercept)          X1          X2          X3 
##    2.999295    1.968817   -6.993233    4.012021
mat_testM4 <- model.matrix(Y ~ X1 + X2 + X3 + X4, data = df_test)
prediccionesM4 <- mat_testM4[, names(coef(M4))] %*% coef(M4)


mat_test_M3 <- model.matrix(Y ~ X1 + X2 + X3, data = df_test)
prediccionesM3 <- mat_test_M3[, names(coef(M3))] %*% coef(M3)


#Comparamos errores cuadraticos medios
ECM_M4 <- mean((df_test$Y - prediccionesM4)^2)
ECM_M3 <- mean((df_test$Y - prediccionesM3)^2)

# Mostrar los resultados
ECM_M4 
## [1] 0.01060874
ECM_M3 
## [1] 0.01064443
cat("El mejor modelo es M4 con un ECM de", ECM_M4)
## El mejor modelo es M4 con un ECM de 0.01060874

Selección progresiva hacia adelante

modelosDelante <- regsubsets(Y~.,df_entrenamiento,nvmax=8,method="forward")
summary(modelosDelante)
## Subset selection object
## Call: regsubsets.formula(Y ~ ., df_entrenamiento, nvmax = 8, method = "forward")
## 8 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
## 1 subsets of each size up to 8
## Selection Algorithm: forward
##          X1  X2  X3  X4  X5  X6  X7  X8 
## 1  ( 1 ) " " " " "*" " " " " " " " " " "
## 2  ( 1 ) " " "*" "*" " " " " " " " " " "
## 3  ( 1 ) "*" "*" "*" " " " " " " " " " "
## 4  ( 1 ) "*" "*" "*" "*" " " " " " " " "
## 5  ( 1 ) "*" "*" "*" "*" "*" " " " " " "
## 6  ( 1 ) "*" "*" "*" "*" "*" " " "*" " "
## 7  ( 1 ) "*" "*" "*" "*" "*" " " "*" "*"
## 8  ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
#hacemos las mismas comparaciones que antes 
nPredictoresR2<-which.max(summary(modelosDelante)$adjr2)
nPredictoresAIC<-which.min(summary(modelosDelante)$cp)
nPredictoresBIC<-which.min(summary(modelosDelante)$bic)

nPredictoresR2
## [1] 4
nPredictoresAIC
## [1] 4
nPredictoresBIC
## [1] 3

Desde aqui podemos notar que el mejor modelo va a ser el mismo que en la mejor seleccion de predictores, puesto que el numero y predictores coinciden con los del inciso anterior

#En otras palabras el mejor modelo dado el metodo de seleccion hacia adelante es M4
summary(M4)
## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_entrenamiento)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28181 -0.05164  0.00891  0.06682  0.22545 
## 
## Coefficients:
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  2.988660   0.011749  254.374   <2e-16 ***
## X1           1.962509   0.016874  116.302   <2e-16 ***
## X2          -6.959993   0.020261 -343.521   <2e-16 ***
## X3           4.015385   0.006178  649.943   <2e-16 ***
## X4          -0.008652   0.004911   -1.762   0.0802 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09984 on 145 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 6.353e+05 on 4 and 145 DF,  p-value: < 2.2e-16

Selección progresiva hacia atras

modelosDetras <- regsubsets(Y~.,df_entrenamiento,nvmax=8,method="backward")
summary(modelosDetras)
## Subset selection object
## Call: regsubsets.formula(Y ~ ., df_entrenamiento, nvmax = 8, method = "backward")
## 8 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
## 1 subsets of each size up to 8
## Selection Algorithm: backward
##          X1  X2  X3  X4  X5  X6  X7  X8 
## 1  ( 1 ) " " " " "*" " " " " " " " " " "
## 2  ( 1 ) " " "*" "*" " " " " " " " " " "
## 3  ( 1 ) "*" "*" "*" " " " " " " " " " "
## 4  ( 1 ) "*" "*" "*" " " " " " " " " "*"
## 5  ( 1 ) "*" "*" "*" " " " " " " "*" "*"
## 6  ( 1 ) "*" "*" "*" " " "*" " " "*" "*"
## 7  ( 1 ) "*" "*" "*" "*" "*" " " "*" "*"
## 8  ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
#hacemos las mismas comparaciones que antes 
nPredictoresR2<-which.max(summary(modelosDetras)$adjr2)
nPredictoresAIC<-which.min(summary(modelosDetras)$cp)
nPredictoresBIC<-which.min(summary(modelosDetras)$bic)

nPredictoresR2
## [1] 4
nPredictoresAIC
## [1] 4
nPredictoresBIC
## [1] 3

En este caso el numero de estimadores optimo es el mismo; sin embargo, M4 cambia y ahora usa X8 como predictor

#creamos el nuevo modelo
M4 <- lm(Y ~ X1 + X2 + X3 + X8, data = df_entrenamiento)

mat_testM4 <- model.matrix(Y ~ X1 + X2 + X3 + X8, data = df_test)
prediccionesM4 <- mat_testM4[, names(coef(M4))] %*% coef(M4)

#Comparamos errores cuadraticos medios
ECM_M4 <- mean((df_test$Y - prediccionesM4)^2)

# Mostrar los resultados
ECM_M4 
## [1] 0.01061142
ECM_M3 
## [1] 0.01064443

Con este método tambien obtenemos que el mejor modelo es de 4 predictores, pero en lugar de X4 tenemos X8

Conclusiones (f)

En este caso los metodos de selección de predictores arrojaron resultados similares. El mejor metodo de seleccion y el metodo hacia adelante nos dieron un modelo de 4 predictores usando X1,X2,X3, y X4, mientras que el metodo hacia atras nos dio uno con el mismo numero de predictores pero con X1,X2,X3, y X8.

En teoria, comparando los ECM de los dos ultimos modelos, el mejor modelo DADA LA MUESTRA TEST es el M4 dado por el mejor metodo de selección.