Forward Stepwise Selection
Este algoritmo es computacionalmente mas eficiente para seleccionar un subset. Mientras otros algoritmos consideran todas las combinaciones posibles de los predictores, forward stepwise considera mucho menos modelos. Comienza con un modelo sin predictores y luego agrega predictores uno a la vez hasta que todos los predictores esten en el modelo. En particular, la variable que es agregada al modelo es aquella que aporte una mejora en el ajuste del modelo.
library(dplyr)
library(MASS)
library(ISLR)
library(caret)
Se crean los folds, en la variable kfolds hay una lista que contiene listas de un train y un test.
set.seed(1)
#devuelve un dataset con la columna fold dependiendo de en cuanto se quiere dividir
kfolds_tag <- function(k, dataset){
folds <- 1:k*rep(1,nrow(dataset))
dataset <- dataset[sample(1:nrow(dataset)),]
dataset$fold <- folds
return(dataset)
}
#Devuelve una lista que contiene el test y el train ya con el fold que se quiere usar como train.
kfoldcv_list <- function(x,dataset){
test <- dataset %>%
filter(fold == x)
train <- dataset %>%
filter(fold != x)
return(list(train,test))
}
k = 5
kfolds <- lapply(1:k,
kfoldcv_list,
dataset = kfolds_tag(k, Boston))
A continuacion se presenta un ejemplo de cross-validation usando kfolds en un modelo de una sola variable.
Vemos que los distintos modelos de una variable que se pueden formar son 13, que es la cantidad de features que tiene el dataset menos \(1\) que es la variable que vamos a predecir.
Se obtienen las variables para los modelos. Se va a predecir medv:
variables <- colnames(Boston)
variables <- setdiff(variables, "medv")
features_subset <-
combn(variables, 1, simplify = FALSE)
features_subset %>% unlist()
[1] "crim" "zn" "indus" "chas" "nox"
[6] "rm" "age" "dis" "rad" "tax"
[11] "ptratio" "black" "lstat"
Funcion que devuelve el Squared Error del modelo.
kfolds_metric <- function(datasets_list, label, feature){
train <- datasets_list[[1]]
test <- datasets_list[[2]]
formula <- paste0(feature, collapse = "+")
formula <- paste0(label, "~",
formula,
collapse = "")
fit <- lm(formula, data = train)
test_pred <- predict(fit, test)
se <- sum((test_pred - test[,label])^2)
return(se)
}
Aplicamos kfolds_metric y obtenemos el error para predecir la variable medv usando la variable crim.
ejemplo <-
lapply(kfolds,
kfolds_metric,
label = "medv",
feature = (features_subset %>% unlist())[1])
Y el MSE es el siguiente:
ejemplo
[[1]]
[1] 6483.559
[[2]]
[1] 7896.243
[[3]]
[1] 6335.429
[[4]]
[1] 7257.351
[[5]]
[1] 8519.34
mean(ejemplo %>% unlist())
[1] 7298.384
La siguiente funcion obtiene los SE de cada fold para un feature en especifico y luego devuelve el MSE.
mse_in_model <- function(variables_subset, kfolds_list, label, i){
specific_feature <- (variables_subset %>% unlist())[i]
se_model_list <-
lapply(kfolds_list,
kfolds_metric,
label = label,
feature = specific_feature)
mse_model <- mean(se_model_list %>% unlist())
return(mse_model)
}
Como podemos ver el MSE usando el feature crim es de \(MSE=7298.38\):
mse_in_model(features_subset, kfolds, "medv", 1)
[1] 7298.384
A continuacion se muestran los MSE para cada uno de los modelos con los distintos \(13\) features:
lapply(features_subset,
mse_in_model,
kfolds_list = kfolds,
label = "medv",
i = 1)
[[1]]
[1] 7298.384
[[2]]
[1] 7446.966
[[3]]
[1] 6579.919
[[4]]
[1] 8500.715
[[5]]
[1] 7019.304
[[6]]
[1] 4505.497
[[7]]
[1] 7372.615
[[8]]
[1] 8051.34
[[9]]
[1] 7330.823
[[10]]
[1] 6693.436
[[11]]
[1] 6376.768
[[12]]
[1] 7607.52
[[13]]
[1] 3917.973
A continuacion elegimos el modelo de 1 variable que menor MSE tiene:
choose_model <- function(variables_subset, kfolds_list, label){
models_mse <-
lapply(variables_subset,
mse_in_model,
kfolds_list = kfolds_list,
label = label,
i = 1)
return(which.min(models_mse %>% unlist()))
}
modelo_menor_mse <- choose_model(features_subset, kfolds, "medv")
modelo_menor_mse
[1] 13
features_subset[modelo_menor_mse]
[[1]]
[1] "lstat"
Como podemos ver es el modelo que contiene el feature lstat. A partir de este modelo se deberian empezar a adicionar variables de una en una y fijar el modelo que menor MSE posea.
LS0tDQp0aXRsZTogIkhvamEgZGUgVHJhYmFqbyA2OiBGb3J3YXJkIGFuZCBCYWNrd2FyZCBTdGVwd2lzZSBTZWxlY3Rpb24iDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQojIyBQYXVsYSBDYXphbGkNCiMjIyBGaWFiaWxpZGFkDQoNCiMjIEZvcndhcmQgU3RlcHdpc2UgU2VsZWN0aW9uDQoNCkVzdGUgYWxnb3JpdG1vIGVzIGNvbXB1dGFjaW9uYWxtZW50ZSBtYXMgZWZpY2llbnRlIHBhcmEgc2VsZWNjaW9uYXIgdW4gc3Vic2V0LiBNaWVudHJhcyBvdHJvcyBhbGdvcml0bW9zIGNvbnNpZGVyYW4gdG9kYXMgbGFzIGNvbWJpbmFjaW9uZXMgcG9zaWJsZXMgZGUgbG9zIHByZWRpY3RvcmVzLCBmb3J3YXJkIHN0ZXB3aXNlIGNvbnNpZGVyYSBtdWNobyBtZW5vcyBtb2RlbG9zLiBDb21pZW56YSBjb24gdW4gbW9kZWxvIHNpbiBwcmVkaWN0b3JlcyB5IGx1ZWdvIGFncmVnYSBwcmVkaWN0b3JlcyB1bm8gYSBsYSB2ZXogaGFzdGEgcXVlIHRvZG9zIGxvcyBwcmVkaWN0b3JlcyBlc3RlbiBlbiBlbCBtb2RlbG8uIEVuIHBhcnRpY3VsYXIsIGxhIHZhcmlhYmxlIHF1ZSBlcyBhZ3JlZ2FkYSBhbCBtb2RlbG8gZXMgYXF1ZWxsYSBxdWUgYXBvcnRlIHVuYSBtZWpvcmEgZW4gZWwgYWp1c3RlIGRlbCBtb2RlbG8uDQoNCmBgYHtyfQ0KbGlicmFyeShkcGx5cikNCmxpYnJhcnkoTUFTUykNCmxpYnJhcnkoSVNMUikNCmxpYnJhcnkoY2FyZXQpDQpgYGANCg0KU2UgY3JlYW4gbG9zIGZvbGRzLCBlbiBsYSB2YXJpYWJsZSBga2ZvbGRzYCBoYXkgdW5hIGxpc3RhIHF1ZSBjb250aWVuZSBsaXN0YXMgZGUgdW4gdHJhaW4geSB1biB0ZXN0LiANCmBgYHtyfQ0Kc2V0LnNlZWQoMSkNCiNkZXZ1ZWx2ZSB1biBkYXRhc2V0IGNvbiBsYSBjb2x1bW5hIGZvbGQgZGVwZW5kaWVuZG8gZGUgZW4gY3VhbnRvIHNlIHF1aWVyZSBkaXZpZGlyDQprZm9sZHNfdGFnIDwtIGZ1bmN0aW9uKGssIGRhdGFzZXQpew0KICBmb2xkcyA8LSAxOmsqcmVwKDEsbnJvdyhkYXRhc2V0KSkNCiAgZGF0YXNldCA8LSBkYXRhc2V0W3NhbXBsZSgxOm5yb3coZGF0YXNldCkpLF0NCiAgZGF0YXNldCRmb2xkIDwtIGZvbGRzDQogIHJldHVybihkYXRhc2V0KQ0KfQ0KDQojRGV2dWVsdmUgdW5hIGxpc3RhIHF1ZSBjb250aWVuZSBlbCB0ZXN0IHkgZWwgdHJhaW4geWEgY29uIGVsIGZvbGQgcXVlIHNlIHF1aWVyZSB1c2FyIGNvbW8gdHJhaW4uDQprZm9sZGN2X2xpc3QgPC0gZnVuY3Rpb24oeCxkYXRhc2V0KXsNCiAgdGVzdCA8LSBkYXRhc2V0ICU+JQ0KICAgIGZpbHRlcihmb2xkID09IHgpDQogIHRyYWluIDwtIGRhdGFzZXQgJT4lIA0KICAgIGZpbHRlcihmb2xkICE9IHgpDQogIHJldHVybihsaXN0KHRyYWluLHRlc3QpKQ0KfQ0KDQprID0gNQ0Ka2ZvbGRzIDwtIGxhcHBseSgxOmssDQogICAgICAgICAgICAgICAgICBrZm9sZGN2X2xpc3QsDQogICAgICAgICAgICAgICAgICBkYXRhc2V0ID0ga2ZvbGRzX3RhZyhrLCBCb3N0b24pKQ0KYGBgDQoNCkEgY29udGludWFjaW9uIHNlIHByZXNlbnRhIHVuIGVqZW1wbG8gZGUgY3Jvc3MtdmFsaWRhdGlvbiB1c2FuZG8ga2ZvbGRzIGVuIHVuIG1vZGVsbyBkZSB1bmEgc29sYSB2YXJpYWJsZS4gDQoNClZlbW9zIHF1ZSBsb3MgZGlzdGludG9zIG1vZGVsb3MgZGUgdW5hIHZhcmlhYmxlIHF1ZSBzZSBwdWVkZW4gZm9ybWFyIHNvbiAxMywgcXVlIGVzIGxhIGNhbnRpZGFkIGRlIGZlYXR1cmVzIHF1ZSB0aWVuZSBlbCBkYXRhc2V0IG1lbm9zICQxJCBxdWUgZXMgbGEgdmFyaWFibGUgcXVlIHZhbW9zIGEgcHJlZGVjaXIuDQoNClNlIG9idGllbmVuIGxhcyB2YXJpYWJsZXMgcGFyYSBsb3MgbW9kZWxvcy4gU2UgdmEgYSBwcmVkZWNpciBgbWVkdmA6DQpgYGB7cn0NCnZhcmlhYmxlcyA8LSBjb2xuYW1lcyhCb3N0b24pDQp2YXJpYWJsZXMgPC0gc2V0ZGlmZih2YXJpYWJsZXMsICJtZWR2IikNCmZlYXR1cmVzX3N1YnNldCA8LSANCiAgY29tYm4odmFyaWFibGVzLCAxLCBzaW1wbGlmeSA9IEZBTFNFKQ0KZmVhdHVyZXNfc3Vic2V0ICU+JSB1bmxpc3QoKQ0KYGBgDQoNCkZ1bmNpb24gcXVlIGRldnVlbHZlIGVsIFNxdWFyZWQgRXJyb3IgZGVsIG1vZGVsby4NCmBgYHtyfQ0Ka2ZvbGRzX21ldHJpYyA8LSBmdW5jdGlvbihkYXRhc2V0c19saXN0LCBsYWJlbCwgZmVhdHVyZSl7DQogIHRyYWluIDwtIGRhdGFzZXRzX2xpc3RbWzFdXQ0KICB0ZXN0IDwtIGRhdGFzZXRzX2xpc3RbWzJdXQ0KICBmb3JtdWxhIDwtIHBhc3RlMChmZWF0dXJlLCBjb2xsYXBzZSA9ICIrIikNCiAgZm9ybXVsYSA8LSBwYXN0ZTAobGFiZWwsICJ+IiwNCiAgICAgICAgICAgICAgICAgICAgZm9ybXVsYSwNCiAgICAgICAgICAgICAgICAgICAgY29sbGFwc2UgPSAiIikNCiAgZml0IDwtIGxtKGZvcm11bGEsIGRhdGEgPSB0cmFpbikNCiAgdGVzdF9wcmVkIDwtIHByZWRpY3QoZml0LCB0ZXN0KQ0KICBzZSA8LSBzdW0oKHRlc3RfcHJlZCAtIHRlc3RbLGxhYmVsXSleMikNCiAgcmV0dXJuKHNlKQ0KfQ0KYGBgDQoNCkFwbGljYW1vcyBga2ZvbGRzX21ldHJpY2AgeSBvYnRlbmVtb3MgZWwgZXJyb3IgcGFyYSBwcmVkZWNpciBsYSB2YXJpYWJsZSBgbWVkdmAgdXNhbmRvIGxhIHZhcmlhYmxlIGBjcmltYC4NCmBgYHtyfQ0KZWplbXBsbyA8LSANCiAgbGFwcGx5KGtmb2xkcywNCiAgICAgICAgIGtmb2xkc19tZXRyaWMsDQogICAgICAgICBsYWJlbCA9ICJtZWR2IiwNCiAgICAgICAgIGZlYXR1cmUgPSAoZmVhdHVyZXNfc3Vic2V0ICU+JSB1bmxpc3QoKSlbMV0pDQpgYGANCg0KWSBlbCBNU0UgZXMgZWwgc2lndWllbnRlOg0KYGBge3J9DQplamVtcGxvDQptZWFuKGVqZW1wbG8gJT4lIHVubGlzdCgpKQ0KYGBgDQoNCkxhIHNpZ3VpZW50ZSBmdW5jaW9uIG9idGllbmUgbG9zIFNFIGRlIGNhZGEgZm9sZCBwYXJhIHVuIGZlYXR1cmUgZW4gZXNwZWNpZmljbyB5IGx1ZWdvIGRldnVlbHZlIGVsIE1TRS4NCmBgYHtyfQ0KbXNlX2luX21vZGVsIDwtIGZ1bmN0aW9uKHZhcmlhYmxlc19zdWJzZXQsIGtmb2xkc19saXN0LCBsYWJlbCwgaSl7DQogIHNwZWNpZmljX2ZlYXR1cmUgPC0gKHZhcmlhYmxlc19zdWJzZXQgJT4lIHVubGlzdCgpKVtpXQ0KICAgIHNlX21vZGVsX2xpc3QgPC0gDQogICAgICBsYXBwbHkoa2ZvbGRzX2xpc3QsDQogICAgICAgICAgICAga2ZvbGRzX21ldHJpYywNCiAgICAgICAgICAgICBsYWJlbCA9IGxhYmVsLA0KICAgICAgICAgICAgIGZlYXR1cmUgPSBzcGVjaWZpY19mZWF0dXJlKQ0KICAgIG1zZV9tb2RlbCA8LSBtZWFuKHNlX21vZGVsX2xpc3QgJT4lIHVubGlzdCgpKQ0KICAgIHJldHVybihtc2VfbW9kZWwpDQp9DQpgYGANCg0KQ29tbyBwb2RlbW9zIHZlciBlbCBNU0UgdXNhbmRvIGVsIGZlYXR1cmUgYGNyaW1gIGVzIGRlICRNU0U9NzI5OC4zOCQ6DQpgYGB7cn0NCm1zZV9pbl9tb2RlbChmZWF0dXJlc19zdWJzZXQsIGtmb2xkcywgIm1lZHYiLCAxKQ0KYGBgDQoNCkEgY29udGludWFjaW9uIHNlIG11ZXN0cmFuIGxvcyBNU0UgcGFyYSBjYWRhIHVubyBkZSBsb3MgbW9kZWxvcyBjb24gbG9zIGRpc3RpbnRvcyAkMTMkIGZlYXR1cmVzOg0KYGBge3J9DQpsYXBwbHkoZmVhdHVyZXNfc3Vic2V0LA0KICAgICAgICAgbXNlX2luX21vZGVsLA0KICAgICAgICAga2ZvbGRzX2xpc3QgPSBrZm9sZHMsDQogICAgICAgICBsYWJlbCA9ICJtZWR2IiwNCiAgICAgICAgIGkgPSAxKQ0KYGBgDQoNCg0KQSBjb250aW51YWNpb24gZWxlZ2ltb3MgZWwgbW9kZWxvIGRlIDEgdmFyaWFibGUgcXVlIG1lbm9yIE1TRSB0aWVuZToNCmBgYHtyfQ0KY2hvb3NlX21vZGVsIDwtIGZ1bmN0aW9uKHZhcmlhYmxlc19zdWJzZXQsIGtmb2xkc19saXN0LCBsYWJlbCl7DQogIG1vZGVsc19tc2UgPC0NCiAgICBsYXBwbHkodmFyaWFibGVzX3N1YnNldCwNCiAgICAgICAgIG1zZV9pbl9tb2RlbCwNCiAgICAgICAgIGtmb2xkc19saXN0ID0ga2ZvbGRzX2xpc3QsDQogICAgICAgICBsYWJlbCA9IGxhYmVsLA0KICAgICAgICAgaSA9IDEpDQogIHJldHVybih3aGljaC5taW4obW9kZWxzX21zZSAlPiUgdW5saXN0KCkpKQ0KfQ0KDQptb2RlbG9fbWVub3JfbXNlIDwtIGNob29zZV9tb2RlbChmZWF0dXJlc19zdWJzZXQsIGtmb2xkcywgIm1lZHYiKQ0KbW9kZWxvX21lbm9yX21zZQ0KZmVhdHVyZXNfc3Vic2V0W21vZGVsb19tZW5vcl9tc2VdDQoNCmBgYA0KDQpDb21vIHBvZGVtb3MgdmVyIGVzIGVsIG1vZGVsbyBxdWUgY29udGllbmUgZWwgZmVhdHVyZSBgbHN0YXRgLg0KQSBwYXJ0aXIgZGUgZXN0ZSBtb2RlbG8gc2UgZGViZXJpYW4gZW1wZXphciBhIGFkaWNpb25hciB2YXJpYWJsZXMgZGUgdW5hIGVuIHVuYSB5IGZpamFyIGVsIG1vZGVsbyBxdWUgbWVub3IgTVNFIHBvc2VhLg0KDQoNCg0KDQo=