Arboles de Regresión
Miguel Arquez Abdala
30/3/2020
Arboles de Regresión
Los arboles de regresión invoucran métodos de segmentación de los datos en grupos más pequeños y luego estimar un modelo simple para cada subgrupo. Aunque los arboles de decisión simples suelen tener un poder predictivo pobre y ser muy inestables (alta varianza). Para ellos hay soluciones como los modelos de agregación por bootstrap(Bagging) o los bosques aleatorios
library("rsample") # data splitting
library("dplyr") # data wrangling
library("rpart") # performing regression trees
library("rpart.plot") # plotting regression trees
library("ipred") # bagging
library("caret") # baggingUsamos para los ejercicios la base de datos AmesHousing para predecir precios de la vivienda basado en multiples características.
# Creamos la partición usual de datos
set.seed(123)
ames_split <- initial_split(AmesHousing::make_ames(), prop = .7)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)Hay multiples metodologías para construir arboles de regresión, sin embargo uno de los más populares es el CART (Classification and regression tree). Arboles de regresión básicos son construidos cuando se realiza una partición en grupos más pequeños y se estima una constante para cada grupo. La partición alcanzada por particiones sucesivas binarias (Partición recursiva) está basada en mutiples variables de predicción. El valor predicho no es más que el promedio de la variable independiente de los elementos que se encuetran en ese grupo.
La idea similar a la regresión por MCO. Buscamos un número de splits que reduzcan la suma de los errores al cuadrado del arbol:
\[ min\{SSE = \sum_{i\in R_1}{(y_i- c_1)^2} + \sum_{i\in R_2}{(y_i- c_2)^2} \} \]
Aqui hay que tomar varias desiciones antes de proceder a contruir un arbol de forma automatica:
- ¿Qué profundida tendrá el arbol?
- ¿Qué tan complejo debe ser el arbol?
Hay una forma de encontrar un modelo que controle la complejidad del arbol sin aumentar la varianza del mismo:
\[ min\{SSE + \alpha|T| \} \] Usamos un parametro de complejidad que penaliza la función objetivo cuando aumenta el número de nodos terminales u hojas (subgrupos finales generados por nuevas particiones) que son agregados innecesariamente (no aumentan el poder predictivo del modelo ).
Para un parametro dado de alfa, encontramos el arbol “podado” más pequeño que tiene el mínimo error penalizado. Pequeñas penalizaciones tienen de producir arboles más comlejos, y por lo tanto más grandes. Alfa puede ser identificado con cross-validation para no subestimar su valor, pero tampoco tan grande que el arboles tenga poca varianza pero alto sesgo.
Ventajas de los arboles de decisión
- Faciles de interpretar y entender
- Hace predicciones muy rápido
- Facil de entender qué variables son más importantes para hacer la predicción
- Funciona bien cuando hay no linealidades en los datos
Desventas arboles de decisión
- tiene alta varianza y poco poder predictivo
m1 <- rpart(
formula = Sale_Price ~ ., # Ecuación var dependiente vs. independientes
data = ames_train, # Dataset
method = "anova" # Anova para especificar que es un arbol de regresión
)
rpart.plot(m1) 
plotcp(m1)
m2 <- rpart(
formula = Sale_Price ~ .,
data = ames_train,
method = "anova",
control = list(cp = 0, xval = 15) # cp = cost complexity - xval = cross validation
)
plotcp(m2) # Error relativo: rescalamiento del error entre 0 y 1
abline(v = 12, lty = "dashed")
Se puede observar que a medida que la penalización es más alta el tamaño del arbol es menor, mientras que con penalizaciones muy bajas el arbol puede crecer en ramas incluso tanto como observaciones hay en el dataset.
También se puede observar que las disminuciones en el error son casi nulas después de un cierto número de arboles. La función rpart() seleccina automaticamente estos valores arrojando el mejor arbol basado en la penalización.
m1$cptable## CP nsplit rel error xerror xstd
## 1 0.46690132 0 1.0000000 1.0009222 0.05855161
## 2 0.11961409 1 0.5330987 0.5347929 0.03116217
## 3 0.06955813 2 0.4134846 0.4151417 0.03058554
## 4 0.02559992 3 0.3439265 0.3461258 0.02207839
## 5 0.02196620 4 0.3183265 0.3242197 0.02182111
## 6 0.02023390 5 0.2963603 0.3074877 0.02129292
## 7 0.01674138 6 0.2761264 0.2963372 0.02106996
## 8 0.01188709 7 0.2593850 0.2795199 0.01903482
## 9 0.01127889 8 0.2474980 0.2762666 0.01936472
## 10 0.01109955 9 0.2362191 0.2699895 0.01902217
## 11 0.01060346 11 0.2140200 0.2672133 0.01883219
## 12 0.01000000 12 0.2034165 0.2635207 0.01881691Aunque aún podemos tunear más el modelo para mejorar el desempeño (después de todo en algunos casos para nosotros podría ser deseable seguir trabajando con modelos sencillos y faciles de intrerpretar).
Además del parametro de complejidad (penalización) podemos controlar el número mínimo de observaciones necesarias para realizar un nuevo split del arbol y la máxima profundidad del mismo. Los valores por defecto para rpart son 20 y 30 respectivamente.
m3 <- rpart(
formula = Sale_Price ~ .,
data = ames_train,
method = "anova",
control = list(minsplit = 10, maxdepth = 12, xval = 10)
)
m3$cptable## CP nsplit rel error xerror xstd
## 1 0.46690132 0 1.0000000 1.0004448 0.05850012
## 2 0.11961409 1 0.5330987 0.5343156 0.03093134
## 3 0.06955813 2 0.4134846 0.4148699 0.03035832
## 4 0.02559992 3 0.3439265 0.3455539 0.02190359
## 5 0.02196620 4 0.3183265 0.3259151 0.02168056
## 6 0.02023390 5 0.2963603 0.3062045 0.02114604
## 7 0.01674138 6 0.2761264 0.3061135 0.02176061
## 8 0.01188709 7 0.2593850 0.2917534 0.02058535
## 9 0.01127889 8 0.2474980 0.2872380 0.02441006
## 10 0.01109955 9 0.2362191 0.2850234 0.02440721
## 11 0.01060346 11 0.2140200 0.2829790 0.02334151
## 12 0.01000000 12 0.2034165 0.2735069 0.02260957Ahora, encontrar los valores optimos podría ser bastante complicado, ya que habría que hacerlo a mano. Para ello podemos usar el enfoque grid search para crear multiples valores de los hiperparametros ( asi son conocidos los valores que no son estimados por el modelo pero necesarios para optimizarlo)
# Creamos un dataframe con las posibles combinaciones de hiperparametros
hyper_grid <- expand.grid(
minsplit = seq(5, 20, 1), # poblaciones entre 5 y 20
maxdepth = seq(8, 15, 1) #profundidad del arbol entre 8 y 15
)
models <- vector(mode = "list", nrow(hyper_grid))
for (i in 1:nrow(hyper_grid)) {
# Obtenemos los hiperparametros de la estimación i
minsplit <- hyper_grid$minsplit[i]
maxdepth <- hyper_grid$maxdepth[i]
# Entrenamos el modelo
models[[i]] <- rpart(
formula = Sale_Price ~ .,
data = ames_train,
method = "anova",
control = list(minsplit = minsplit, maxdepth = maxdepth)
)
}Podemos ver una mejoría con respecto al modelo sin optimización
# Creamos una función para extraer el cost complexity de la lis de modelos
get_cp <- function(x) {
min <- which.min(x$cptable[, "xerror"])
cp <- x$cptable[min, "CP"]
}
# Función para extraer el modelo con el error minimo
get_min_error <- function(x) {
min <- which.min(x$cptable[, "xerror"])
xerror <- x$cptable[min, "xerror"]
}
# Miramos el top 5 de valore minimos
hyper_grid %>%
mutate(
cp = purrr::map_dbl(models, get_cp),
error = purrr::map_dbl(models, get_min_error)
) %>%
arrange(error) %>%
top_n(-5, wt = error)## minsplit maxdepth cp error
## 1 16 12 0.01060346 0.2628987
## 2 6 11 0.01000000 0.2645615
## 3 11 11 0.01000000 0.2650862
## 4 10 10 0.01000000 0.2655860
## 5 7 15 0.01000000 0.2656602Finalmente aplicamos los valores del arbol optimo y estiamos sobre los datos de testeo.
optimal_tree <- rpart(
formula = Sale_Price ~ .,
data = ames_train,
method = "anova",
control = list(minsplit = 7, maxdepth = 12, cp = 0.01060346)
)
pred <- predict(optimal_tree, newdata = ames_test)
RMSE(pred = pred, obs = ames_test$Sale_Price)## [1] 39852.01Bagging
Aunque el modelo de arbol de regresión podría ser bueno, recordemos que sufre de alta varianza. incluso cuando el modelo está optimizado. En tal caso podemos ir un paso más allá y mejorar el desempeño del modelo con bagging (boostrap aggregating).
Esta metodología consiste en crear m muestras booststrap de los datos de entrenamiento. A cada muestrale estimamos un arbol de regresión sin ninguna optimización ni parametro de control y finalmente predecimos con los promedios de cada arbol para crear una predicción definitiva.
Un beneficio del bagging es que, en promedio, una muestra bootstrap contendrá unicamente 63% de los datos de entrenamiento, esto deja otro 37% de los datos fuera la muestra boostsrappeada, que se denomina out-of-bag (OOB) y sobre ellas hacemos los análisis de desempeño del modelo
# Semilla
set.seed(123)
# Entrenamos el modelo
bagged_m1 <- bagging(
formula = Sale_Price ~ .,
data = ames_train,
coob = TRUE # Errores sobre el OOB
)
bagged_m1##
## Bagging regression trees with 25 bootstrap replications
##
## Call: bagging.data.frame(formula = Sale_Price ~ ., data = ames_train,
## coob = TRUE)
##
## Out-of-bag estimate of root mean squared error: 36991.67Por defecto la función bagging crea 25 muestrar para el modelo, pero este número está sujeto a optimizaciones conforme a el comportamiento del error (Aunque por lo general entre más arboles menor error y varianza, hasta cierto punto)
# Evaluamos el modelo entre 10 y 50 arboles
ntree <- 10:50
# create empty vector to store OOB RMSE values
rmse <- vector(mode = "numeric", length = length(ntree))
for (i in seq_along(ntree)) {
# reproducibility
set.seed(123)
# perform bagged model
model <- bagging(
formula = Sale_Price ~ .,
data = ames_train,
coob = TRUE,
nbagg = ntree[i]
)
# get OOB error
rmse[i] <- model$err
}
plot(ntree, rmse, type = 'l', lwd = 2)
abline(v = 25, col = "red", lty = "dashed")
Aunque realizar bagging con caret es mucho mejor:
- Calcula la importancia de las variables: se evalua calculando el monto total se RSS que disminiye por hacer un split sobre una variable en particular, promediado por el número total de arboles. La variable predictora con el mayor impacto sobre el RSS es considerada la más importante. El valor de importancias es simplemente el decrecimiento promedio del RSS comparado con el el valor más alto (Reescalamiento de 0 a 100)
# Specify 10-fold cross validation
ctrl <- trainControl(method = "cv", number = 10)
# CV bagged model
bagged_cv <- train(
Sale_Price ~ .,
data = ames_train,
method = "treebag",
trControl = ctrl,
importance = TRUE
)
# assess results
bagged_cv## Bagged CART
##
## 2051 samples
## 80 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 1846, 1845, 1846, 1845, 1847, 1847, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 35854.02 0.8009063 23785.85## Bagged CART
##
## 2051 samples
## 80 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 1846, 1845, 1847, 1845, 1846, 1847, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 36477.25 0.8001783 24059.85
# plot most important variables
plot(varImp(bagged_cv), 20) 
pred <- predict(bagged_cv, ames_test)
RMSE(pred, ames_test$Sale_Price)## [1] 35357.89Este documento es una copia traducida del artículo publicado por la universidad de Cincinnnati: http://uc-r.github.io/regression_trees