Model Selection For Regression

6 minute read

Seleccion de Modelos Exhaustiva Forward y Backward (BIC,CP,R2ajustado)

Ejercicios

1.- Generación de datos simulados y aplicación de los métodos de selección de subconjuntos (selección de los mejores subconjuntos y selección paso a paso).

  1. Usa una función en R para generar una variable predictora X de longitud n = 100 , asi como un vector de ruido \(\epsilon\) de tamaño n = 100.

  2. Genera un vector de respuestas Y de longitud n = 100 de acuerdo al modelo \[Y = \beta_0 + \beta_1X + \beta_2X^2 +\beta_3X^3 +\epsilon\] donde \(\beta_0 \beta_1 \beta_2 \beta_3\) son constantes de tu elección

  3. Utiliza la función regsubsets() para realizar la selección de los mejores subconjuntos con el n de elegir el mejor modelo que contenga los predictores \(X\), \(X^2\),..\(X^{10}\) . ¿Cuál es el mejor modelo obtenido según el C P , BIC y el \(R^2\) ajustado?, Muestra algunas gráficas que proporcionen evidencia de tu respuesta y reporta los coefcientes del mejor modelo obtenido.

  4. Repite (c), usando la selección paso a paso adelante y la selección paso paso atrás. ¿ Cómo se compara tu respuesta con los resultados obtenidos en (c)?

### a)
rm(list = ls())
# t=2
# Z<-function(x)
# {
#   2*x*x+(1/x)+1
#   #rnorm(100)
# }

####(a)####
#X<-sapply(1:100, function(x){Z(x/5)})
set.seed(5)
X<-rnorm(100)
set.seed(6)
epsilon<-rnorm(100,mean = 2,sd = 2)

####(b)####
B0<-1;B1<-3;B2<-5;B3<-0.5
#Y = B0 + B1*X*X + B2*X*X +B3*X*X*X + epsilon 
Y = B0 + B1*X + B2*X*X +B3*X*X*X + epsilon 
plot(Y)

####(c)####
library(leaps)
?regsubsets()
Xmat<-matrix(ncol = 10,nrow = 100)
for (i in 1:10){
  Xmat[,i]<-X**i
}
dim(Xmat)
## [1] 100  10
Xmat_df<-as.data.frame(Xmat)
#cbind(Xmat[,1],X)
#cbind(Xmat[,2],X)

mejores_modelos_ex<-regsubsets(Y~.,data =Xmat_df,nvmax = 10,method =  "exhaustive",nbest = 1)
summ_ex<-summary(mejores_modelos_ex)
#summ_ex$adjr2
#summ_ex$cp
#summ_ex$bic
data.frame(ex_adjr2=summ_ex$adjr2,ex_cp=summ_ex$cp,ex_bic=summ_ex$bic)
library(ggplot2)

p <- ggplot(data = data.frame(n_predictores = 1:10,
                              bic_ex = summ_ex$bic[1:10]),
            aes(x = n_predictores, y = bic_ex)) +
  geom_line(color="green") +
  geom_point(color="green")+
  xlab("# Predictores")+
  ylab("Valor BIC")
p

p <- ggplot(data = data.frame(n_predictores = 1:10,
                              cp_ex = summ_ex$cp[1:10]),
            aes(x = n_predictores, y = cp_ex)) +
  geom_line(color="red") +
  geom_point(color="red")+
  xlab("# Predictores")+
  ylab("Valor CP")
p

p <- ggplot(data = data.frame(n_predictores = 1:10,
                              adjr2_ex = summ_ex$adjr2[1:10]),
            aes(x = n_predictores, y = adjr2_ex)) +
  geom_line(color="blue") +
  geom_point(color="blue")+
  xlab("# Predictores")+
  ylab("Valor R^2 ajustado")
p

#Matriz de modelos seleccionados por # de variables
#summ_ex$which
summ_ex$which
##    (Intercept)    V1   V2    V3    V4    V5    V6    V7    V8    V9   V10
## 1         TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2         TRUE  TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 3         TRUE  TRUE TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 4         TRUE  TRUE TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## 5         TRUE FALSE TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE
## 6         TRUE  TRUE TRUE FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE
## 7         TRUE  TRUE TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
## 8         TRUE  TRUE TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## 9         TRUE  TRUE TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## 10        TRUE  TRUE TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
# Mejores modelos
best_ex_adjr2<-which.max(summ_ex$adjr2)
best_ex_bic<-which.min(summ_ex$bic)
best_ex_cp<-which.min(summ_ex$cp)

#Encontarar coeficientes de los modelos 
coef_ex_adjr2<-coef(mejores_modelos_ex,best_ex_adjr2)
coef_ex_bic<-coef(mejores_modelos_ex,best_ex_bic)
coef_ex_cp<-coef(mejores_modelos_ex,best_ex_cp)

coef_ex_adjr2
##  (Intercept)           V1           V2           V5           V7 
##  2.626821468  2.406967613  5.454770870  1.509440964 -0.648478914 
##           V9          V10 
##  0.074987462 -0.001701273
coef_ex_bic
## (Intercept)          V1          V2          V3 
##   2.8193841   2.8152432   5.1686084   0.5914385
coef_ex_cp
##  (Intercept)           V1           V2           V5           V7 
##  2.626821468  2.406967613  5.454770870  1.509440964 -0.648478914 
##           V9          V10 
##  0.074987462 -0.001701273
# p <- ggplot(data = data.frame(n_predictores = 2:10,
#                               adjr2_fw = summ_fw$adjr2[2:10]),
#             aes(x = n_predictores, y = adjr2_fw)) +
#   geom_line() +
#   geom_point()
# p


####(d) Ahora con forward backawrd step####

 mejores_modelos_fw<-regsubsets(Y~.,data =Xmat_df,nvmax = 10,method =  "forward")
 #summary(mejores_modelos_fw)
 summ_fw<- summary(mejores_modelos_fw)
 #summ_fw$adjr2
 #summ_fw$cp
 #summ_fw$bic
 data.frame(fw_adjr2=summ_fw$adjr2,fw_cp=summ_fw$cp,fw_bic=summ_fw$bic)
 p <- ggplot(data = data.frame(n_predictores = 1:10,
                               bic_fw = summ_fw$bic[1:10]),
             aes(x = n_predictores, y = bic_fw)) +
   geom_line(color="green") +
   geom_point(color="green")+
   xlab("# Predictores")+
   ylab("Valor BIC")
 p

 p <- ggplot(data = data.frame(n_predictores = 1:10,
                               cp_fw = summ_fw$cp[1:10]),
             aes(x = n_predictores, y = cp_fw)) +
   geom_line(color="red") +
   geom_point(color="red")+
   xlab("# Predictores")+
   ylab("Valor CP")
 p