Model Selection For Regression
Seleccion de Modelos Exhaustiva Forward y Backward (BIC,CP,R2ajustado)
Adrián Alejandro Rodríguez Villarreal
12 de septiembre de 2018
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).
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.
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
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.
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 10Xmat_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.001701273coef_ex_bic## (Intercept) V1 V2 V3
## 2.8193841 2.8152432 5.1686084 0.5914385coef_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
p <- ggplot(data = data.frame(n_predictores = 1:10,
adjr2_fw = summ_fw$adjr2[1:10]),
aes(x = n_predictores, y = adjr2_fw)) +
geom_line(color="blue") +
geom_point(color="blue")+
xlab("# Predictores")+
ylab("Valor R^2 ajustado")
p
#Matriz de modelos seleccionados por # de variables
#summ_fw$which
summ_fw$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 TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE
## 6 TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE TRUE FALSE TRUE
## 7 TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE
## 8 TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE
## 9 TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE
## 10 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE # Mejores modelos
best_fw_adjr2<-which.max(summ_fw$adjr2)
best_fw_bic<-which.min(summ_fw$bic)
best_fw_cp<-which.min(summ_fw$cp)
#Encontarar coeficientes de los modelos
coef_fw_adjr2<-coef(mejores_modelos_fw,best_fw_adjr2)
coef_fw_bic<-coef(mejores_modelos_fw,best_fw_bic)
coef_fw_cp<-coef(mejores_modelos_fw,best_fw_cp)
coef_fw_adjr2## (Intercept) V1 V2 V3 V4 V5
## 2.59181719 2.79061681 5.93389195 -1.28454793 -0.47883758 2.69219967
## V7 V8 V9 V10
## -1.02436211 0.06396386 0.11319110 -0.01167515 coef_fw_bic## (Intercept) V1 V2 V3
## 2.8193841 2.8152432 5.1686084 0.5914385 coef_fw_cp## (Intercept) V1 V2 V3
## 2.8193841 2.8152432 5.1686084 0.5914385 #### BACKWARD ####
mejores_modelos_bw<-regsubsets(Y~.,data =Xmat_df,nvmax = 10,method = "backward")
#summary(mejores_modelos_bw)
summ_bw<- summary(mejores_modelos_bw)
# summ_bw$adjr2
# summ_bw$cp
# summ_bw$bic
data.frame(bw_adjr2=summ_bw$adjr2,bw_cp=summ_bw$cp,bw_bic=summ_bw$bic) p <- ggplot(data = data.frame(n_predictores = 1:10,
bic_bw = summ_bw$bic[1:10]),
aes(x = n_predictores, y = bic_bw)) +
geom_line(color="green") +
geom_point(color="green")+
xlab("# Predictores")+
ylab("Valor BIC")
p
p <- ggplot(data = data.frame(n_predictores = 1:10,
cp_bw = summ_bw$cp[1:10]),
aes(x = n_predictores, y = cp_bw)) +
geom_line(color="red") +
geom_point(color="red")+
xlab("# Predictores")+
ylab("Valor CP")
p
p <- ggplot(data = data.frame(n_predictores = 1:10,
adjr2_bw = summ_bw$adjr2[1:10]),
aes(x = n_predictores, y = adjr2_bw)) +
geom_line(color="blue") +
geom_point(color="blue")+
xlab("# Predictores")+
ylab("Valor R^2 ajustado")
p
#Matriz de modelos seleccionados por # de variables
summ_bw$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 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE
## 4 TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE
## 5 TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE TRUE FALSE
## 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_bw_adjr2<-which.max(summ_bw$adjr2)
best_bw_bic<-which.min(summ_bw$bic)
best_bw_cp<-which.min(summ_bw$cp)
#Encontarar coeficientes de los modelos
coef_bw_adjr2<-coef(mejores_modelos_bw,best_bw_adjr2)
coef_bw_bic<-coef(mejores_modelos_bw,best_bw_bic)
coef_bw_cp<-coef(mejores_modelos_bw,best_bw_cp)
coef_bw_adjr2## (Intercept) V1 V2 V5 V7
## 2.626821468 2.406967613 5.454770870 1.509440964 -0.648478914
## V9 V10
## 0.074987462 -0.001701273 coef_bw_bic## (Intercept) V1 V2 V5
## 2.88679470 3.43147143 5.12739087 0.09577905 coef_bw_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 = 1:10,
adjr2_bw = summ_bw$adjr2[1:10],
cp_bw = summ_bw$cp[1:10],
bic_bw = summ_bw$bic[1:10])) +
geom_line(aes(x = n_predictores, y = adjr2_bw,color="green")) +
xlab("# Predictores")+
ylab("Valor Metrica")+
geom_point(aes(x = n_predictores, y = adjr2_bw,color="green"))+
geom_line(aes(x = n_predictores, y = cp_bw,color="red")) +
geom_point(aes(x = n_predictores, y = cp_bw,color="red"))+
geom_line(aes(x = n_predictores, y = bic_bw,color="blue")) +
geom_point(aes(x = n_predictores, y = bic_bw,color="blue"))+
scale_color_discrete(name = "Accuracy Metrics", labels = c("CP", "adj_r2","BIC"))
p
#When picking from several models, the one with the lowest BIC is preferred.
# The highest adj_r2 is preferred
# the lowest CP is preferred
###https://rpubs.com/Joaquin_AR/242707 (Usar regsubsets y graficas de regsubsets)
#http://jadianes.me/best-subset-model-selection-with-R