Más sobre ciencia de datos: cienciadedatos.net

Introducción

Los árboles de decisión son modelos predictivos formados por reglas binarias (si/no) con las que se consigue repartir las observaciones en función de sus atributos y predecir así el valor de la variable respuesta.

Muchos métodos predictivos generan modelos globales en los que una única ecuación se aplica a todo el espacio muestral. Cuando el caso de uso implica múltiples predictores, que interaccionan entre ellos de forma compleja y no lineal, es muy difícil encontrar un único modelo global que sea capaz de reflejar la relación entre las variables. Los métodos estadísticos y de machine learning basados en árboles engloban a un conjunto de técnicas supervisadas no paramétricas que consiguen segmentar el espacio de los predictores en regiones simples, dentro de las cuales es más sencillo manejar las interacciones. Es esta característica la que les proporciona gran parte de su potencial.

Los modelos Random Forest están formados por un conjunto de árboles de decisión individuales, cada uno entrenado con una muestra ligeramente distinta de los datos de entrenamiento generada mediante bootstrapping. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.

Los modelos Gradient Boosting están formados por un conjunto de árboles de decisión individuales, entrenados de forma secuencial, de forma que cada nuevo árbol trata de mejorar los errores de los árboles anteriores. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.

Los métodos basados en árboles se han convertido en uno de los referentes dentro del ámbito predictivo debido a los buenos resultados que generan en problemas muy diversos. A lo largo de este documento se explora la forma en que se construyen y predicen los árboles de decisión (clasificación y regresión), elementos fundamentales de modelos predictivos más complejos como Random Forest y Gradient Boosting.

Ventajas

Los árboles son fáciles de interpretar aun cuando las relaciones entre predictores son complejas.
Los modelos basados en un solo árbol (no es el caso de Random Forest y Boosting) se pueden representar gráficamente aun cuando el número de predictores es mayor de 3.
Los árboles pueden, en teoría, manejar tanto predictores numéricos como categóricos sin tener que crear variables dummy o one-hot-encoding. En la práctica, esto depende de la implementación del algoritmo que tenga cada librería.
Al tratarse de métodos no paramétricos, no es necesario que se cumpla ningún tipo de distribución específica.
Por lo general, requieren mucha menos limpieza y preprocesado de los datos en comparación con otros métodos de aprendizaje estadístico (por ejemplo, no requieren estandarización).
No se ven muy influenciados por outliers.
Si para alguna observación, el valor de un predictor no está disponible, a pesar de no poder llegar a ningún nodo terminal, se puede conseguir una predicción empleando todas las observaciones que pertenecen al último nodo alcanzado. La precisión de la predicción se verá reducida pero al menos podrá obtenerse.
Son muy útiles en la exploración de datos, permiten identificar de forma rápida y eficiente las variables (predictores) más importantes.
Son capaces de seleccionar predictores de forma automática.
Pueden aplicarse a problemas de regresión y clasificación.

Desventajas

La capacidad predictiva de los modelos basados en un único árbol es bastante inferior a la conseguida con otros modelos. Esto es debido a su tendencia al overfitting y alta varianza. Sin embargo, existen técnicas más complejas que, haciendo uso de la combinación de múltiples árboles (bagging, random forest, boosting), consiguen mejorar en gran medida este problema.
Son sensibles a datos de entrenamiento desbalanceados (una de las clases domina sobre las demás).
Cuando tratan con predictores continuos, pierden parte de su información al categorizarlos en el momento de la división de los nodos.
Tal y como se describe más adelante, la creación de las ramificaciones de los árboles se consigue mediante el algoritmo de recursive binary splitting. Este algoritmo identifica y evalúa las posibles divisiones de cada predictor acorde a una determinada medida (RSS, Gini, entropía…). Los predictores continuos tienen mayor probabilidad de contener, solo por azar, algún punto de corte óptimo, por lo que suelen verse favorecidos en la creación de los árboles.
No son capaces de extrapolar fuera del rango de los predictores observado en los datos de entrenamiento.

Paquetes R

Debido a sus buenos resultados, Random Forest y Gradient Boosting, se han convertido en los algoritmos de referencia cuando se trata con datos tabulares, de ahí que se hayan desarrollado múltiples implementaciones. Cada una tiene unas características que las hacen más adecuadas dependiendo del caso de uso. Algunos de los paquetes más empleados en R son:

tree y rpart: permiten crear y representar árboles de regresión y clasificación.
randomForest: dispone de los principales algoritmos para crear modelos random forest. Destaca por su fácil uso, pero no por su rapidez.
ranger: algoritmos para crear modelos random forest. Es similar a randomForest pero mucho más rápido. Además, incorpora extremely randomized trees y quantile regression forests.
gbm: dispone de los principales algoritmos de boosting. Este paquete ya no está mantenido, aunque es útil para explicar los conceptos, no se recomienda su uso en producción.
XGBoost: esta librería permite acceder al algoritmo XGboost (Extra Gradient boosting). Una adaptación de gradient boosting que destaca por su eficiencia y rapidez.
H2O: implementaciones muy optimizadas de los principales algoritmos de machine learning, entre ellos random forest, gradient boosting y XGBoost. Para conocer más detalles sobre cómo utilizar esta librería consultar Machine Learning con H2O y R.
C50: implementación de los algoritmos C5.0 para árboles de clasificación.

Árboles de regresión

Idea intuitiva

Los árboles de regresión son el subtipo de árboles de predicción que se aplica cuando la variable respuesta es continua. En términos generales, en el entrenamiento de un árbol de regresión, las observaciones se van distribuyendo por bifurcaciones (nodos) generando la estructura del árbol hasta alcanzar un nodo terminal. Cuando se quiere predecir una nueva observación, se recorre el árbol acorde al valor de sus predictores hasta alcanzar uno de los nodos terminales. La predicción del árbol es la media de la variable respuesta de las observaciones de entrenamiento que están en ese mismo nodo terminal.

La forma más sencilla de entender la idea detrás de los árboles de regresión es mediante de un ejemplo simplificado. El set de datos Hitter contiene información sobre 322 jugadores de béisbol de la liga profesional. Entre las variables registradas para cada jugador se encuentran: el salario (Salary), años de experiencia (Years) y el número de bateos durante los últimos años (Hits). Utilizando estos datos, se quiere predecir el salario (en unidades logarítmicas) de un jugador en base a su experiencia y número de bateos. El árbol resultante se muestra en la siguiente imagen:

La interpretación del árbol se hace en sentido descendente. La primera división es la que separa a los jugadores en función de si superan o no los 4.5 años de experiencia. En la rama izquierda quedan aquellos con menos de 4.5 años, y su salario predicho se corresponde con el salario promedio de todos los jugadores de este grupo. La rama derecha se subdivide en función de si superan o no 117.5 bateos. Como resultado de las bifurcaciones, se han generado 3 regiones que pueden identificarse con la siguiente nomenclatura:

\(R_1 = \{X| Year < 4.5\}\): jugadores que han jugado menos de 4.5 años.
\(R_2 = \{X| Year \geq 4.5, Hits < 117.5\}\): jugadores que han jugado 4.5 años o más y que han conseguido menos de 117.5 bateos.
\(R_3 = \{X| Year \geq 4.5, Hits \geq 117.5\}\): jugadores que han jugado 4.5 años o más y que han conseguido 117.5 o más bateos.

A las regiones \(R_1\), \(R_2\) y \(R_3\) se les conoce como nodos terminales (terminal nodes o leaves), a los puntos en los que el espacio de los predictores sufre una división como nodos internos (internal nodes o splits) y a los segmentos que conectan dos nodos como ramas (branches).

Viendo el árbol de la imagen anterior, la interpretación del modelo es: la variable más importante a la hora de determinar el salario de un jugador es el número de años jugados, los jugadores con más experiencia ganan más. Entre los jugadores con pocos años de experiencia, el número de bateos logrados en los años previos no tiene mucho impacto en el salario, sin embargo, sí lo tiene para jugadores con cuatro años y medio o más de experiencia. Para estos últimos, a mayor número de bateos logrados, mayor salario. Con este ejemplo queda patente que, la principal ventaja de los árboles de decisión frente a otros métodos de regresión, es su fácil interpretación y la gran utilidad de su representación gráfica.

Entrenamiento del árbol

El proceso de entrenamiento de un árbol de predicción (regresión o clasificación) se divide en dos etapas:

División sucesiva del espacio de los predictores generando regiones no solapantes (nodos terminales) \(R_1\), \(R_2\), \(R_3\), …, \(R_j\). Aunque, desde el punto de vista teórico las regiones podrían tener cualquier forma, si se limitan a regiones rectangulares (de múltiples dimensiones), se simplifica en gran medida el proceso de construcción y se facilita la interpretación.
Predicción de la variable respuesta en cada región.

A pesar de la sencillez con la que se puede resumir el proceso de construcción de un árbol, es necesario establecer una metodología que permita crear las regiones \(R_1\), \(R_2\), \(R_3\), …, \(R_j\), o lo que es equivalente, decidir donde se introducen las divisiones: en qué predictores y en qué valores de los mismos. Es en este punto donde se diferencian los algoritmos de árboles de regresión y clasificación.

En el caso de los árboles de regresión, el criterio empleado con más frecuencia para identificar las divisiones es el Residual Sum of Squares (RSS). El objetivo es encontrar las \(J\) regiones \((R_1\),…, \(R_j)\) que minimizan el Residual Sum of Squares (RSS) total:

\[RSS=\displaystyle\sum_{j=1}^J \displaystyle\sum_{i \in R_j}(y_i -\hat{y}_{R_j})^2 ,\]

donde \(\hat{y}_{R_j}\) es la media de la variable respuesta en la región \(R_j\). Una descripción menos técnica equivale a decir que: se busca una distribución de regiones tal que, el sumatorio de las desviaciones al cuadrado entre las observaciones y la media de la región a la que pertenecen sea lo menor posible.

Desafortunadamente, no es posible considerar todas las posibles particiones del espacio de los predictores. Por esta razón, se recurre a lo que se conoce como recursive binary splitting (división binaria recursiva). Esta solución sigue la misma idea que la selección de predictores stepwise (backward o fordward) en regresión lineal múltiple, no evalúa todas las posibles regiones pero, alcanza un buen balance computación-resultado.

Recursive binary splitting

El objetivo del método recursive binary splitting es encontrar en cada iteración el predictor \(X_j\) y el punto de corte (umbral) \(s\) tal que, si se distribuyen las observaciones en las regiones \(\{X|X_j < s \}\) y \(\{X|X_j \geq s \}\), se consigue la mayor reducción posible en el RSS. El algoritmo seguido es:

El proceso se inicia en lo más alto del árbol, donde todas las observaciones pertenecen a la misma región.
Se identifican todos los posibles puntos de corte (umbrales) \(s\) para cada uno de los predictores (\(X_1\), \(X_2\),…, \(X_p\)). En el caso de predictores cualitativos, los posibles puntos de corte son cada uno de sus niveles. Para predictores continuos, se ordenan de menor a mayor sus valores, el punto intermedio entre cada par de valores se emplea como punto de corte.
Se calcula el RSS total que se consigue con cada posible división identificada en el paso 2.

\[\sum_{i: x_i \in R_1(j,s)} (y_i -\hat{y}_{R_1})^2 + \sum_{i: x_i \in R_2(j,s)} (y_i -\hat{y}_{R_2})^2\]

donde el primer término es el RSS de la región 1 y el segundo término es el RSS de la región 2, siendo cada una de las regiones el resultado de separar las observaciones acorde al predictor \(j\) y valor \(s\).

Se selecciona el predictor \(X_j\) y el punto de corte \(S\) que resulta en el menor RSS total, es decir, que da lugar a las divisiones más homogéneas posibles. Si existen dos o más divisiones que consiguen la misma mejora, la elección entre ellas es aleatoria.
Se repiten de forma iterativa los pasos 1 a 4 para cada una de las regiones que se han creado en la iteración anterior hasta que se alcanza alguna norma de stop. Algunas de las más empleadas son: que ninguna región contenga un mínimo de n observaciones, que el árbol tenga un máximo de nodos terminales o que la incorporación del nodo reduzca el error en al menos un % mínimo.

Está metodología conlleva dos hechos:

Que cada división óptima se identifica acorde al impacto que tiene en ese momento. No se tiene en cuenta si es la división que dará lugar a mejores árboles en futuras divisiones.
En cada división se evalúa un único predictor haciendo preguntas binarias (si, no), lo que genera dos nuevas ramas del árbol por división. A pesar de que es posible evaluar divisiones más complejas, hacer una pregunta sobre múltiples variables a la vez es equivalente a hacer múltiples preguntas sobre variables individuales.

Nota: Los algoritmos que implementan recursive binary splitting suelen incorporar estrategias para evitar evaluar todos los posibles puntos de corte. Por ejemplo, para predictores continuos, primero se crea un histograma que agrupa los valores y luego se evalúan los puntos de corte de cada región del histograma (bin).

Predicción del árbol

Tras la creación de un árbol, las observaciones de entrenamiento quedan agrupadas en los nodos terminales. Para predecir una nueva observación, se recorre el árbol en función de los valores que tienen sus predictores hasta llegar a uno de los nodos terminales. En el caso de regresión, el valor predicho suele ser la media de la variable respuesta de las observaciones de entrenamiento que están en ese mismo nodo. Si bien la media es el valor más empleado, se puede utilizar cualquier otro (mediana, cuantil…).

Supóngase que se dispone de 10 observaciones, cada una con un valor de variable respuesta \(Y\) y unos predictores \(X\).

La siguiente imagen muestra cómo sería la predicción del árbol para una nueva observación. El camino hasta llegar al nodo final está resaltado. En cada nodo terminal se detalla el índice de las observaciones de entrenamiento que forman parte.

Predicción con un árbol de regresión: el camino hasta llegar al nodo final está resaltado. En cada nodo terminal se detalla el índice de las observaciones de entrenamiento que forman parte.

El valor predicho por el árbol es la media de la variable respuesta \(Y\) de las observaciones con \(id\): 2, 3, 5, 9.

\[\hat{\mu} = \frac{18+24+2+20}{4} = 16\]

Aunque la anterior es la forma más común de obtener las predicciones de un árbol de regresión, existe otra aproximación. La predicción de un árbol de regresión puede verse como una variante de vecinos cercanos en la que, solo las observaciones que forman parte del mismo nodo terminal que la observación predicha, tienen influencia. Siguiendo esta aproximación, la predicción del árbol se define como la media ponderada de todas las observaciones de entrenamiento, donde el peso de cada observación depende únicamente de si forma parte o no del mismo nodo terminal.

\[\hat{\mu} = \sum^n_{i =1} \textbf{w}_iY_i\]

El valor de las posiciones del vector de pesos \(\textbf{w}\) es \(1\) para las observaciones que están en el mismo nodo y \(0\) para el resto. En este ejemplo sería:

\[\textbf{w} =(0, 1, 1, 0, 1, 0, 0, 0, 1, 0)\]

Para que la suma de todos los pesos sea \(1\), se dividen por el número total de observaciones en el nodo terminal seleccionado, en este caso 4.

\[\textbf{w} = (0, \frac{1}{4}, \frac{1}{4}, 0, \frac{1}{4}, 0, 0, 0, \frac{1}{4}, 0)\]

Así pues, el valor predicho es:

\[\hat{\mu} = (0\times10) + (\frac{1}{4}\times18) + (\frac{1}{4}\times24) + (0\times8) + (\frac{1}{4}\times2) + \\ (0\times9) + (0\times16) + (0\times10) + (\frac{1}{4}\times20) + (0\times14) = 16\]

Evitar el overfitting

El proceso de construcción de árboles descrito en las secciones anteriores tiende a reducir rápidamente el error de entrenamiento, es decir, el modelo se ajusta muy bien a las observaciones empleadas como entrenamiento. Como consecuencia, se genera un overfitting que reduce su capacidad predictiva al aplicarlo a nuevos datos. La razón de este comportamiento radica en la facilidad con la que los árboles se ramifican adquiriendo estructuras complejas. De hecho, si no se limitan las divisiones, todo árbol termina ajustándose perfectamente a las observaciones de entrenamiento creando un nodo terminal por observación. Existen dos estrategias para prevenir el problema de overfitting de los árboles: limitar el tamaño del árbol (parada temprana o early stopping) y el proceso de podado (pruning).

Controlar el tamaño del árbol (parada temprana)

El tamaño final que adquiere un árbol puede controlarse mediante reglas de parada que detengan la división de los nodos dependiendo de si se cumplen o no determinadas condiciones. El nombre de estas condiciones puede variar dependiendo del software o librería empleada, pero suelen estar presentes en todos ellos.

Observaciones mínimas para división: define el número mínimo de observaciones que debe tener un nodo para poder ser dividido. Cuanto mayor el valor, menos flexible es el modelo.
Observaciones mínimas de nodo terminal: define el número mínimo de observaciones que deben tener los nodos terminales. Su efecto es muy similar al de observaciones mínimas para división.
Profundidad máxima del árbol: define la profundidad máxima del árbol, entendiendo por profundidad máxima el número de divisiones de la rama más larga (en sentido descendente) del árbol.
Número máximo de nodos terminales: define el número máximo de nodos terminales que puede tener el árbol. Una vez alcanzado el límite, se detienen las divisiones. Su efecto es similar al de controlar la profundidad máxima del árbol.
Reducción mínima de error: define la reducción mínima de error que tiene que conseguir una división para que se lleve a cabo.

A todos estos parámetros se les conoce como hiperparámetros porque no se aprenden durante el entrenamiento del modelo. Su valor tiene que ser especificado por el usuario en base a su conocimiento del problema y mediante el uso de validación cruzada.

Tree pruning

La estrategia de controlar el tamaño del árbol mediante reglas de parada tiene un inconveniente, el árbol se crece seleccionando la mejor división en cada momento. Al evaluar las divisiones sin tener en cuenta las que vendrán después, nunca se elige la opción que resulta en el mejor árbol final, a no ser que también sea la que genera en ese momento la mejor división. A este tipo de estrategias se les conoce como greedy. Un ejemplo que ilustra el problema de este tipo de estrategia es el siguiente: supóngase que un coche circula por el carril izquierdo de una carretera de dos carriles en la misma dirección. En el carril que se encuentra hay muchos coches circulando a 100 km/h, mientras que el otro carril se encuentra vacío. A cierta distancia se observa que hay un vehículo circulando por el carril derecho a 20 km/h. Si el objetivo del conductor es llegar a su destino lo antes posible tiene dos opciones: cambiarse de carril o mantenerse en el que está. Una aproximación de tipo greedy evaluaría la situación en ese instante y determinaría que la mejor opción es cambiarse de carril y acelerar a más de 100 km/h, sin embargo, a largo plazo, esta no es la mejor solución, ya que una vez alcance al vehículo lento, tendrá que reducir mucho su velocidad.

Una alternativa no greedy que consigue evitar el overfitting consiste en generar árboles grandes, sin condiciones de parada más allá de las necesarias por las limitaciones computacionales, para después podarlos (pruning), manteniendo únicamente la estructura robusta que consigue un test error bajo. La selección del sub-árbol óptimo puede hacerse mediante cross-validation, sin embargo, dado que los árboles se crecen lo máximo posible (tienen muchos nodos terminales) no suele ser viable estimar el test error de todas las posibles sub-estructuras que se pueden generar. En su lugar, se recurre al cost complexity pruning o weakest link pruning.

Cost complexity pruning es un método de penalización de tipo Loss + Penalty, similar al empleado en ridge regression o lasso. En este caso, se busca el sub-árbol \(T\) que minimiza la ecuación:

\[\displaystyle\sum_{j=1}^{|T|} \displaystyle\sum_{i \in R_j}(y_i -\hat{y}_{R_j})^2 + \alpha|T| \]

donde \(|T|\) es el número de nodos terminales del árbol.

El primer término de la ecuación se corresponde con el sumatorio total de los residuos cuadrados RSS. Por definición, cuantos más nodos terminales tenga el modelo menor será esta parte de la ecuación. El segundo término es la restricción, que penaliza al modelo en función del número de nodos terminales (a mayor número, mayor penalización). El grado de penalización se determina mediante el tuning parameter \(\alpha\). Cuando \(\alpha = 0\), la penalización es nula y el árbol resultante es equivalente al árbol original. A medida que se incrementa \(\alpha\) la penalización es mayor y, como consecuencia, los árboles resultantes son de menor tamaño. El valor óptimo de \(\alpha\) puede identificarse mediante cross validation.

Algoritmo para crear un árbol de regresión con pruning

Se emplea recursive binary splitting para crear un árbol grande y complejo (\(T_o\)) empleando los datos de training y reduciendo al máximo posible las condiciones de parada. Normalmente se emplea como única condición de parada el número mínimo de observaciones por nodo terminal.
Se aplica el cost complexity pruning al árbol \(T_o\) para obtener el mejor sub-árbol en función de \(\alpha\). Es decir, se obtiene el mejor sub-árbol para un rango de valores de \(\alpha\).
Identificación del valor óptimo de \(\alpha\) mediante k-cross-validation. Se divide el training data set en \(K\) grupos. Para \(k=1\), …, \(k=K\):
- Repetir pasos 1 y 2 empleando todas las observaciones excepto las del grupo \(k_i\).
- Evaluar el mean squared error para el rango de valores de \(\alpha\) empleando el grupo \(k_i\).
- Obtener el promedio de los \(K\) mean squared error calculados para cada valor \(\alpha\).
Seleccionar el sub-árbol del paso 2 que se corresponde con el valor \(\alpha\) que ha conseguido el menor cross-validation mean squared error en el paso 3.

En el caso de los árboles de clasificación (ver más adelante), en lugar de emplear la suma de residuos cuadrados como criterio de selección, se emplea alguna de las medidas de homogeneidad vistas los siguientes apartados.

Ejemplo regresión

Librerías

Las librerías utilizadas en este ejemplo son:

# Tratamiento de datos
# ==============================================================================
library(MASS)
library(dplyr)
library(tidyr)
library(skimr)

# Gráficos
# ==============================================================================
library(ggplot2)
library(ggpubr)

# Preprocesado y modelado
# ==============================================================================
library(tree)

Datos

El set de datos Boston disponible en el paquete MASS contiene precios de viviendas de la ciudad de Boston, así como información socioeconómica del barrio en el que se encuentran. Se pretende ajustar un modelo de regresión que permita predecir el precio medio de una vivienda (medv) en función de las variables disponibles.

data("Boston")
datos = Boston
head(datos, 3)

skim(datos)

Data summary
Name	datos
Number of rows	506
Number of columns	14
_______________________
Column type frequency:
numeric	14
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
crim	1	3.61	8.60	0.01	0.08	0.26	3.68	88.98	▇▁▁▁▁
zn	1	11.36	23.32	0.00	0.00	0.00	12.50	100.00	▇▁▁▁▁
indus	1	11.14	6.86	0.46	5.19	9.69	18.10	27.74	▇▆▁▇▁
chas	1	0.07	0.25	0.00	0.00	0.00	0.00	1.00	▇▁▁▁▁
nox	1	0.55	0.12	0.38	0.45	0.54	0.62	0.87	▇▇▆▅▁
rm	1	6.28	0.70	3.56	5.89	6.21	6.62	8.78	▁▂▇▂▁
age	1	68.57	28.15	2.90	45.02	77.50	94.07	100.00	▂▂▂▃▇
dis	1	3.80	2.11	1.13	2.10	3.21	5.19	12.13	▇▅▂▁▁
rad	1	9.55	8.71	1.00	4.00	5.00	24.00	24.00	▇▂▁▁▃
tax	1	408.24	168.54	187.00	279.00	330.00	666.00	711.00	▇▇▃▁▇
ptratio	1	18.46	2.16	12.60	17.40	19.05	20.20	22.00	▁▃▅▅▇
black	1	356.67	91.29	0.32	375.38	391.44	396.22	396.90	▁▁▁▁▇
lstat	1	12.65	7.14	1.73	6.95	11.36	16.96	37.97	▇▇▅▂▁
medv	1	22.53	9.20	5.00	17.02	21.20	25.00	50.00	▂▇▅▁▁

Ajuste del modelo

La función tree() del paquete tree permite ajustar árboles de decisión. La elección entre árbol de regresión o árbol de clasificación se hace automáticamente dependiendo de si la variable respuesta es de tipo numeric o factor. Es importante tener en cuenta que solo estos dos tipos de vectores están permitidos, si se pasa uno de tipo character se devuelve un error.

A continuación, se ajusta un árbol de regresión empleando como variable respuesta medv y como predictores todas las variables disponibles. Como en todo estudio de regresión, no solo es importante ajustar el modelo, sino también cuantificar su capacidad para predecir nuevas observaciones. Para poder hacer la posterior evaluación, se dividen los datos en dos grupos, uno de entrenamiento y otro de test.

La función tree() crece el árbol hasta que encuentra una condición de stop. Por defecto, estas condiciones son:

mincut: número mínimo de observaciones que debe de tener al menos uno de los nodos hijos para que se produzca la división.
minsize: número mínimo de observaciones que debe de tener un nodo para que pueda dividirse.

# División de los datos en train y test
# ==============================================================================
set.seed(123)
train <- sample(1:nrow(datos), size = nrow(datos)/2)
datos_train <- datos[train,]
datos_test  <- datos[-train,]

# Creación y entrenamiento del modelo
# ==============================================================================
set.seed(123)
arbol_regresion <- tree::tree(
                    formula = medv ~ .,
                    data    = datos_train,
                    split   = "deviance",
                    mincut  = 20,
                    minsize = 50
                  )

summary(arbol_regresion)

## 
## Regression tree:
## tree::tree(formula = medv ~ ., data = datos_train, split = "deviance", 
##     mincut = 20, minsize = 50)
## Variables actually used in tree construction:
## [1] "rm"    "lstat" "dis"   "tax"  
## Number of terminal nodes:  6 
## Residual mean deviance:  20.56 = 5078 / 247 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -14.5500  -2.8680  -0.3628   0.0000   2.0050  22.1300

La función summary() muestra que, el árbol entrenado, tiene un total de 6 nodos terminales y que se han empleado como predictores las variables rm, lstat, dis y tax. En el contexto de árboles de regresión, el término Residual mean deviance es la suma de cuadrados residuales dividida entre (número de observaciones - número de nodos terminales). Cuanto menor es la deviance, mejor es el ajuste del árbol a las observaciones de entrenamiento.

Una vez creado el árbol, se puede representar mediante la combinación de las funciones plot() y text(). La función plot() dibuja la estructura del árbol, las ramas y los nodos. Mediante su argumento type se puede especificar si se quiere que todas las ramas tengan el mismo tamaño (type = "uniform") o que su longitud sea proporcional a la reducción de impureza (heterogeneidad) de los nodos terminales (type = "proportional"). Esta segunda opción permite identificar visualmente el impacto de cada división en el modelo. La función text() añade la descripción de cada nodo interno y el valor de cada nodo terminal.

# Estructura del árbol creado
# ==============================================================================
par(mar = c(1,1,1,1))
plot(x = arbol_regresion, type = "proportional")
text(x = arbol_regresion, splits = TRUE, pretty = 0, cex = 0.8, col = "firebrick")

Podado del árbol (pruning)

Con la finalidad de reducir la varianza del modelo y así mejorar la capacidad predictiva, se somete al árbol a un proceso de pruning. Tal como se describió con anterioridad, el proceso de pruning intenta encontrar el árbol más sencillo (menor tamaño) que consigue los mejores resultados de predicción.

Para aplicar el proceso de pruning es necesario indicar el grado de penalización por complejidad \(\alpha\). Cuanto mayor es este valor, más agresivo es el podado y menor el tamaño del árbol resultante. Dado que no hay forma de conocer de antemano el valor óptimo de \(\alpha\), se recurre a validación cruzada para identificarlo.

La función cv.tree() emplea cross validation para identificar el valor óptimo de penalización. Por defecto, esta función emplea la deviance para guiar el proceso de pruning.

# Pruning (const complexity pruning) por validación cruzada
# ==============================================================================

# El árbol se crece al máximo posible para luego aplicar el pruning
arbol_regresion <- tree(
                    formula = medv ~ .,
                    data    = datos_train,
                    split   = "deviance",
                    mincut  = 1,
                    minsize = 2,
                    mindev  = 0
                  )

# Búsqueda por validación cruzada
set.seed(123)
cv_arbol <- cv.tree(arbol_regresion, K = 5)

El objeto devuelto por cv.tree() contiene:

size: el tamaño (número de nodos terminales) de cada árbol.
dev: la estimación de cross-validation test error para cada tamaño de árbol.
k: El rango de valores de penalización \(\alpha\) evaluados.
method: El criterio empleado para seleccionar el mejor árbol.

# Tamaño óptimo encontrado
# ==============================================================================
size_optimo <- rev(cv_arbol$size)[which.min(rev(cv_arbol$dev))]
paste("Tamaño óptimo encontrado:", size_optimo)

## [1] "Tamaño óptimo encontrado: 10"

resultados_cv <- data.frame(
                   n_nodos  = cv_arbol$size,
                   deviance = cv_arbol$dev,
                   alpha    = cv_arbol$k
                 )

p1 <- ggplot(data = resultados_cv, aes(x = n_nodos, y = deviance)) +
      geom_line() + 
      geom_point() +
      geom_vline(xintercept = size_optimo, color = "red") +
      labs(title = "Error vs tamaño del árbol") +
      theme_bw() 
  
p2 <- ggplot(data = resultados_cv, aes(x = alpha, y = deviance)) +
      geom_line() + 
      geom_point() +
      labs(title = "Error vs penalización alpha") +
      theme_bw() 

ggarrange(p1, p2)

Una vez identificado el valor óptimo de, con la función prune.tree() se aplica el podado final. Esta función también acepta el valor de \(\alpha\) óptimo en lugar del tamaño.

# Estructura del árbol creado final
# ==============================================================================
arbol_final <- prune.tree(
                  tree = arbol_regresion,
                  best = size_optimo
               )

par(mar = c(1,1,1,1))
plot(x = arbol_final, type = "proportional")
text(x = arbol_final, splits = TRUE, pretty = 0, cex = 0.8, col = "firebrick")

Predicción y evaluación del modelo

Por último, se evalúa la capacidad predictiva del primer árbol y del árbol final empleando el conjunto de test.

# Error de test del modelo inicial
# ==============================================================================
predicciones <- predict(arbol_regresion, newdata = datos_test)
test_rmse    <- sqrt(mean((predicciones - datos_test$medv)^2))
paste("Error de test (rmse) del árbol inicial:", round(test_rmse,2))

## [1] "Error de test (rmse) del árbol inicial: 5.43"

# Error de test del modelo final
# ==============================================================================
predicciones <- predict(arbol_final, newdata = datos_test)
test_rmse    <- sqrt(mean((predicciones - datos_test$medv)^2))
paste("Error de test (rmse) del árbol final:", round(test_rmse,2))

## [1] "Error de test (rmse) del árbol final: 5.13"

El proceso de pruning consigue reducir el error (rmse) del modelo de 5.43 a 5.13. Las predicciones del modelo final se alejan en promedio 5.13 unidades (5130 dólares) del valor real.

Árboles de clasificación

Idea intuitiva

Los árboles de regresión son el subtipo de árboles de predicción que se aplica cuando la variable respuesta es categórica.

Construcción del árbol

Para construir un árbol de clasificación se emplea el mismo método recursive binary splitting descrito en los árboles de regresión. Sin embargo, como la variable respuesta en cualitativa, no es posible emplear el RSS como criterio de selección de las divisiones óptimas. Existen varias alternativas, todas ellas con el objetivo de encontrar nodos lo más puros/homogéneos posible. Las más empleadas son:

Classification Error Rate

Se define como la proporción de observaciones que no pertenecen a la clase más común en el nodo. \[E_m = 1 - max_k(\hat{p}_{mk}), \] donde \(\hat{p}_{mk}\) representa la proporción de observaciones del nodo \(m\) que pertenecen a la clase \(k\). A pesar de la sencillez de esta medida, no es suficientemente sensible para crear buenos árboles, por lo que, en la práctica, suelen emplearse otras medidas.

Gini Index

Es una medida de la varianza total en el conjunto de las \(K\) clases del nodo \(m\). Se considera una medida de pureza del nodo. \[G_m = \sum_{k=1}^K \hat{p}_{mk}(1-\hat{p}_{mk})\] Cuando \(\hat{p}_{mk}\) es cercano a 0 o a 1 (el nodo contiene mayoritariamente observaciones de una clase), el término \(\hat{p}_{mk}(1-\hat{p}_{mk})\) es muy pequeño. Como consecuencia, cuanto mayor sea la pureza del nodo, menor el valor del índice Gini \(G\).

El algoritmo CART (Classification and regression trees) emplea Gini como criterio de división.

Information Gain: Cross Entropy

La entropía es otra forma de cuantificar el desorden de un sistema. En el caso de los nodos, el desorden se corresponde con la impureza. Si un nodo es puro, contiene únicamente observaciones de una clase, su entropía es cero. Por el contrario, si la frecuencia de cada clase es la misma, el valor de la entropía alcanza el valor máximo de 1.

\[D = - \displaystyle\sum_{k=1}^K \hat{p}_{mk} \ log(\hat{p}_{mk})\] Los algoritmos C4.5 y C5.0 emplean information gain como criterio de división.

Chi-Square

Esta aproximación consiste en identificar si existe una diferencia significativa entre los nodos hijos y el nodo parental, es decir, si hay evidencias de que la división consigue una mejora. Para ello, se aplica un test estadístico chi-square goodness of fit empleando como distribución esperada \(H_0\) la frecuencia de cada clase en el nodo parental. Cuanto mayor el estadístico \(\chi^2\), mayor la evidencia estadística de que existe una diferencia.

\[\chi^2 = \sum_{k} \frac{(\text{observado}_{k} - \text{esperado}_{k})^2}{\text{esperado}_{k}}\]

Los árboles generados con este criterio de división reciben el nombre de CHAID (Chi-square Automatic Interaction Detector).

Independientemente de la medida empleada como criterio de selección de las divisiones, el proceso siempre es el mismo:

Para cada posible división se calcula el valor de la medida en cada uno de los dos nodos resultantes.
Se suman los dos valores ponderando cada uno por la fracción de observaciones que contiene cada nodo. Este paso es muy importante, ya que no es lo mismo dos nodos puros con 2 observaciones, que dos nodos puros con 100 observaciones.

\[ \frac{\text{n observaciones nodo A}}{\text{n observaciones totales}} \times \text{pureza A} + \frac{\text{n observaciones nodo B}}{\text{n observaciones totales}} \times\text{pureza B}\]

La división con menor o mayor valor (dependiendo de la medida empleada) se selecciona como división óptima.

Para el proceso de construcción del árbol, acorde al libro Introduction to Statistical Learning, Gini index y cross-entropy son más adecuados que el classification error rate debido a su mayor sensibilidad a la homogeneidad de los nodos. Para el proceso de pruning (descrito en la siguiente sección) los tres son adecuados, aunque, si el objetivo es conseguir la máxima precisión en las predicciones, mejor emplear el classification error rate.

Predicción del árbol

Tras la creación de un árbol, las observaciones de entrenamiento quedan agrupadas en los nodos terminales. Para predecir una nueva observación, se recorre el árbol en función del valor de sus predictores hasta llegar a uno de los nodos terminales. En el caso de clasificación, suele emplearse la moda de la variable respuesta como valor de predicción, es decir, la clase más frecuente del nodo. Además, puede acompañarse con el porcentaje de cada clase en el nodo terminal, lo que aporta información sobre la confianza de la predicción.

Ejemplo clasificación

Librerías

Las librerías utilizadas en este ejemplo son:

# Tratamiento de datos
# ==============================================================================
library(dplyr)
library(tidyr)
library(skimr)
library(ISLR)

# Gráficos
# ==============================================================================
library(ggplot2)
library(ggpubr)

# Preprocesado y modelado
# ==============================================================================
library(tree)

Datos

El set de datos Carseats, disponible en paquete ISLR, contiene información sobre la venta de sillas infantiles en 400 tiendas distintas. Para cada una de las 400 tiendas se han registrado 11 variables. Se pretende generar un modelo de clasificación que permita predecir si una tienda tiene ventas altas (Sales > 8) o bajas (Sales <= 8) en función de todas las variables disponibles.

data("Carseats")
datos = Carseats
head(datos, 3)

skim(datos)

Data summary
Name	datos
Number of rows	400
Number of columns	11
_______________________
Column type frequency:
factor	3
numeric	8
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
ShelveLoc	1	FALSE	3	Med: 219, Bad: 96, Goo: 85
Urban	1	FALSE	2	Yes: 282, No: 118
US	1	FALSE	2	Yes: 258, No: 142

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Sales	1	7.50	2.82	0	5.39	7.49	9.32	16.27	▁▆▇▃▁
CompPrice	1	124.97	15.33	77	115.00	125.00	135.00	175.00	▁▅▇▃▁
Income	1	68.66	27.99	21	42.75	69.00	91.00	120.00	▇▆▇▆▅
Advertising	1	6.64	6.65	0	0.00	5.00	12.00	29.00	▇▃▃▁▁
Population	1	264.84	147.38	10	139.00	272.00	398.50	509.00	▇▇▇▇▇
Price	1	115.80	23.68	24	100.00	117.00	131.00	191.00	▁▂▇▆▁
Age	1	53.32	16.20	25	39.75	54.50	66.00	80.00	▇▆▇▇▇
Education	1	13.90	2.62	10	12.00	14.00	16.00	18.00	▇▇▃▇▇

Como Sales es una variable continua y el objetivo del estudio es clasificar las tiendas según si venden mucho o poco, se crea una nueva variable dicotómica (Si, No) llamada ventas_altas.

datos$ventas_altas <- ifelse(test = datos$Sales > 8, 
                             yes = "Si", no = "No")

# Conversión de la variable respuesta a tipo factor
datos$ventas_altas <- as.factor(datos$ventas_altas)

# Una vez creada la nueva variable respuesta se descarta la original
datos$Sales = NULL

Ajuste del modelo

Se ajusta un árbol de clasificación empleando como variable respuesta ventas_altas y como predictores todas las variables disponibles. Se utilizan en primer lugar minsize=10, el resto de los hiperparámetros se dejan por defecto. Después, se aplica el proceso de pruning y se comparan los resultados frente al modelo inicial.

# División de los datos en train y test
# ==============================================================================
set.seed(123)
train <- sample(1:nrow(datos), nrow(datos) / 2, replace = FALSE)
datos_train <- datos[train,]
datos_test  <- datos[-train,]

# Entrenamiento del modelo
# ==============================================================================
set.seed(123)
arbol_clasificacion <- tree(
                        formula = ventas_altas ~ .,
                        data    = datos_train,
                        minsize = 10
                       )
summary(arbol_clasificacion)

## 
## Classification tree:
## tree(formula = ventas_altas ~ ., data = datos_train, minsize = 10)
## Variables actually used in tree construction:
## [1] "ShelveLoc"   "Price"       "Income"      "CompPrice"   "Advertising"
## [6] "Age"         "Education"   "Population"  "US"         
## Number of terminal nodes:  20 
## Residual mean deviance:  0.5056 = 91.02 / 180 
## Misclassification error rate: 0.115 = 23 / 200

La función summary() muestra que el árbol ajustado tiene un total de 20 nodos terminales y un classification error rate de entrenamiento del 11.5%. La deviance de los árboles de clasificación se calcula como \(-2 \displaystyle\sum_m \displaystyle\sum_k n_{mk}log(\hat{p}_{mk})\), donde \(n_{mk}\) es el número de observaciones en el nodo terminal \(m\) que pertenecen a la clase \(k\). El término Residual mean deviance mostrado en el summary es simplemente la deviance dividida entre (número de observaciones - número de nodos terminales). Cuanto menor es la deviance mejor es el ajuste del árbol a las observaciones de entrenamiento.

# Estructura del árbol creado
# ==============================================================================
par(mar = c(1,1,1,1))
plot(x = arbol_clasificacion, type = "proportional")
text(x = arbol_clasificacion, splits = TRUE, pretty = 0, cex = 0.7, col = "firebrick")

Predicción y evaluación del modelo

Se evalúa la capacidad predictiva del árbol inicial calculando el accuracy en el conjunto de test.

# Error de test del modelo
# ==============================================================================
predicciones <- predict(
                  arbol_clasificacion,
                  newdata = datos_test,
                  type    = "class"
                )
table(predicciones, datos_test$ventas_altas)

##             
## predicciones No Si
##           No 87 13
##           Si 35 65

El modelo inicial es capaz de predecir correctamente un 76 % de las observaciones del conjunto de test.

Podado del árbol (pruning)

Con la finalidad de reducir la varianza del modelo y así disminuir el test error, se somete al árbol a un proceso de pruning. A diferencia del ejemplo anterior, al ser este un árbol de clasificación se indica en la función cv.tree() que FUN = prune.misclass. De esta manera se emplea el classification error rate para guiar el proceso de pruning.

# Pruning (const complexity pruning) por validación cruzada
# ==============================================================================

# El árbol se crece al máximo posible para luego aplicar el pruning
arbol_clasificacion <- tree(
                          formula = ventas_altas ~ .,
                          data    = datos_train,
                          mincut  = 1,
                          minsize = 2,
                          mindev  = 0
                       )

# Búsqueda por validación cruzada
set.seed(123)
cv_arbol <- cv.tree(arbol_clasificacion, FUN = prune.misclass, K = 5)

# Tamaño óptimo encontrado
# ==============================================================================
size_optimo <- rev(cv_arbol$size)[which.min(rev(cv_arbol$dev))]
size_optimo

## [1] 33

resultados_cv <- data.frame(n_nodos = cv_arbol$size, clas_error = cv_arbol$dev,
                            alpha = cv_arbol$k)

p1 <- ggplot(data = resultados_cv, aes(x = n_nodos, y = clas_error)) +
      geom_line() + 
      geom_point() +
      geom_vline(xintercept = size_optimo, color = "red") +
      labs(title = " Error de clasificación vs \n tamaño del árbol") +
      theme_bw() 
  
p2 <- ggplot(data = resultados_cv, aes(x = alpha, y = clas_error)) +
      geom_line() + 
      geom_point() +
      labs(title = " Error de clasificación vs \n penalización alpha") +
      theme_bw() 

ggarrange(p1, p2)

En este caso, un árbol con 33 nodos terminales minimiza el test error.

Con la función prune.missclas() se obtiene el mejor árbol de clasificación del tamaño identificado como óptimo (no confundir con la función prune.tree() para árboles de regresión). A esta función también se le puede indicar el valor de \(\alpha\) óptimo en lugar del tamaño.

arbol_final <- prune.misclass(
                  tree = arbol_clasificacion,
                  best = size_optimo
               )

# Error de test del modelo final
#-------------------------------------------------------------------------------
predicciones <- predict(arbol_clasificacion, newdata = datos_test, type = "class")
table(predicciones, datos_test$ventas_altas)

##             
## predicciones No Si
##           No 95 16
##           Si 27 62

paste("El porcentaje de acierto es de",
      100 * ((95 + 62) / (95 + 62 + 16 + 27)), "%")

## [1] "El porcentaje de acierto es de 78.5 %"

Gracias al proceso de pruning el porcentaje de acierto ha pasado de 76% a 78.5%.

Random Forest

Idea intuitiva

Un modelo Random Forest está formado por un conjunto (ensemble) de árboles de decisión individuales, cada uno entrenado con una muestra aleatoria extraída de los datos de entrenamiento originales mediante bootstrapping. Esto implica que cada árbol se entrena con unos datos ligeramente distintos. En cada árbol individual, las observaciones se van distribuyendo por bifurcaciones (nodos) generando la estructura del árbol hasta alcanzar un nodo terminal. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.

Para entender cómo funcionan los modelos Random Forest es necesario conocer primero los conceptos de ensemble y bagging.

Métodos de ensemble

Todos los modelos de aprendizaje estadístico y machine learning sufren el problema de equilibrio entre bias y varianza.

El término bias (sesgo) hace referencia a cuánto se alejan en promedio las predicciones de un modelo respecto a los valores reales. Refleja cómo de capaz es el modelo de aprender la relación real que existe entre los predictores y la variable respuesta. Por ejemplo, si la relación sigue un patrón no lineal, por muchos datos de los que se disponga, un modelo de regresión lineal no podrá modelar correctamente la relación, por lo que tendrá un bias alto.

El término varianza hace referencia a cuánto cambia el modelo dependiendo de los datos utilizados en su entrenamiento. Idealmente, un modelo no debería modificarse demasiado por pequeñas variaciones en los datos de entrenamiento, si esto ocurre, es porque el modelo está memorizando los datos en lugar de aprender la verdadera relación entre los predictores y la variable respuesta. Por ejemplo, un modelo de árbol con muchos nodos, suele variar su estructura con que apenas cambien unos pocos datos de entrenamiento, tiene mucha varianza.

A medida que aumenta la complejidad de un modelo, este dispone de mayor flexibilidad para adaptarse a las observaciones, reduciendo así el bias y mejorando su capacidad predictiva. Sin embargo, alcanzado un determinado grado de flexibilidad, aparece el problema de overfitting, el modelo se ajusta tanto a los datos de entrenamiento que es incapaz de predecir correctamente nuevas observaciones. El mejor modelo es aquel que consigue un equilibrio óptimo entre bias y varianza.

¿Cómo se controlan el bias y varianza en los modelos basados en árboles? Por lo general, los árboles pequeños (pocas ramificaciones) tienen poca varianza pero no consiguen representar bien la relación entre las variables, es decir, tienen bias alto. En contraposición, los árboles grandes se ajustan mucho a los datos de entrenamiento, por lo que tienen muy poco bias pero mucha varianza. Una forma de solucionar este problema son los métodos de ensemble.

Los métodos de ensemble combinan múltiples modelos en uno nuevo con el objetivo de lograr un equilibro entre bias y varianza, consiguiendo así mejores predicciones que cualquiera de los modelos individuales originales. Dos de los tipos de ensemble más utilizados son:

Bagging: Se ajustan múltiples modelos, cada uno con un subconjunto distinto de los datos de entrenamiento. Para predecir, todos los modelos que forman el agregado participan aportando su predicción. Como valor final, se toma la media de todas las predicciones (variables continuas) o la clase más frecuente (variables categóricas). Los modelos Random Forest están dentro de esta categoría.
Boosting: Se ajustan secuencialmente múltiples modelos sencillos, llamados weak learners, de forma que cada modelo aprende de los errores del anterior. Como valor final, al igual que en bagging, se toma la media de todas las predicciones (variables continuas) o la clase más frecuente (variables cualitativas). Tres de los métodos de boosting más empleados son AdaBoost, Gradient Boosting y Stochastic Gradient Boosting.

Aunque el objetivo final es el mismo, lograr un balance óptimo entre bias y varianza, existen dos diferencias importantes:

Forma en que consiguen reducir el error total. El error total de un modelo puede descomponerse como \(bias + varianza + \epsilon\). En bagging, se emplean modelos con muy poco bias pero mucha varianza, agregándolos se consigue reducir la varianza sin apenas inflar el bias. En boosting, se emplean modelos con muy poca varianza pero mucho bias, ajustando secuencialmente los modelos se reduce el bias. Por lo tanto, cada una de las estrategias reduce una parte del error total.
Forma en que se introducen variaciones en los modelos que forman el ensemble. En bagging, cada modelo es distinto del resto porque cada uno se entrena con una muestra distinta obtenida mediante bootstrapping. En boosting, los modelos se ajustan secuencialmente y la importancia (peso) de las observaciones va cambiando en cada iteración, dando lugar a diferentes ajustes.

La clave para que los métodos de ensemble consigan mejores resultados que cualquiera de sus modelos individuales es que, los modelos que los forman, sean lo más diversos posibles (sus errores no estén correlacionados). Una analogía que refleja este concepto es la siguiente: supóngase un juego como el trivial en el que los equipos tienen que acertar preguntas sobre temáticas diversas. Un equipo formado por muchos jugadores, cada uno experto en un tema distinto, tendrá más posibilidades de ganar que un equipo formado por jugadores expertos en un único tema o por un único jugador que sepa un poco de todos los temas.

A continuación, se describe con más detalle la estrategia de bagging, sobre la que se fundamenta el modelo Random Forest.

Bagging

El término bagging es el diminutivo de bootstrap aggregation, y hace referencia al empleo del muestreo repetido con reposición bootstrapping con el fin de reducir la varianza de algunos modelos de aprendizaje estadístico, entre ellos los basados en árboles.

Dadas \(n\) muestras de observaciones independientes \(Z_1\), …, \(Z_n\), cada una con varianza \(\sigma^2\), la varianza de la media de las observaciones \(\overline{Z}\) es \(\sigma^2/n\). En otras palabras, promediando un conjunto de observaciones se reduce la varianza. Basándose en esta idea, una forma de reducir la varianza y aumentar la precisión de un método predictivo es obtener múltiples muestras de la población, ajustar un modelo distinto con cada una de ellas, y hacer la media (la moda en el caso de variables cualitativas) de las predicciones resultantes. Como en la práctica no se suele tener acceso a múltiples muestras, se puede simular el proceso recurriendo a bootstrapping, generando así pseudo-muestras con los que ajustar diferentes modelos y después agregarlos. A este proceso se le conoce como bagging y es aplicable a una gran variedad de métodos de regresión.

En el caso particular de los árboles de decisión, dada su naturaleza de bajo bias y alta varianza, bagging ha demostrado tener muy buenos resultados. La forma de aplicarlo es:

Generar \(B\) pseudo-training sets mediante bootstrapping a partir de la muestra de entrenamiento original.
Entrenar un árbol con cada una de las \(B\) muestras del paso 1. Cada árbol se crea sin apenas restricciones y no se somete a pruning, por lo que tiene varianza alta pero poco bias. En la mayoría de casos, la única regla de parada es el número mínimo de observaciones que deben tener los nodos terminales. El valor óptimo de este hiperparámetro puede obtenerse comparando el out of bag error o por validación cruzada.
Para cada nueva observación, obtener la predicción de cada uno de los \(B\) árboles. El valor final de la predicción se obtiene como la media de las \(B\) predicciones en el caso de variables cuantitativas y como la clase predicha más frecuente (moda) para variables cualitativas.

En el proceso de bagging, el número de árboles creados no es un hiperparámetro crítico en cuanto a que, por mucho que se incremente el número, no se aumenta el riesgo de overfitting. Alcanzado un determinado número de árboles, la reducción de test error se estabiliza. A pesar de ello, cada árbol ocupa memoria, por lo que no conviene almacenar más de los necesarios.

Entrenamiento de Random Forest

El algoritmo de Random Forest es una modificación del proceso de bagging que consigue mejorar los resultados gracias a que decorrelaciona aún más los árboles generados en el proceso.

Recordando el apartado anterior, los beneficios del bagging se basan en el hecho de que, promediando un conjunto de modelos, se consigue reducir la varianza. Esto es cierto siempre y cuando los modelos agregados no estén correlacionados. Si la correlación es alta, la reducción de varianza que se puede lograr es pequeña.

Supóngase un set de datos en el que hay un predictor muy influyente, junto con otros moderadamente influyentes. En este escenario, todos o casi todos los árboles creados en el proceso de bagging estarán dominados por el mismo predictor y serán muy parecidos entre ellos. Como consecuencia de la alta correlación entre los árboles, el proceso de bagging apenas conseguirá disminuir la varianza y, por lo tanto, tampoco mejorar el modelo. Random forest evita este problema haciendo una selección aleatoria de \(m\) predictores antes de evaluar cada división. De esta forma, un promedio de \((p-m)/p\) divisiones no contemplarán el predictor influyente, permitiendo que otros predictores puedan ser seleccionados. Añadiendo este paso extra se consigue decorrelacionar los árboles todavía más, con lo que su agregación consigue una mayor reducción de la varianza.

Los métodos de random forest y bagging siguen el mismo algoritmo con la única diferencia de que, en random forest, antes de cada división, se seleccionan aleatoriamente \(m\) predictores. La diferencia en el resultado dependerá del valor \(m\) escogido. Si \(m=p\) los resultados de random forest y bagging son equivalentes. Algunas recomendaciones son:

La raíz cuadrada del número total de predictores para problemas de clasificación. \(m \approx \sqrt{p}\)
Un tercio del número de predictores para problemas de regresión. \(m \approx \frac{p}{3}\)
Si los predictores están muy correlacionados, valores pequeños de \(m\) consiguen mejores resultados.

Sin embargo, la mejor forma para encontrar el valor óptimo de \(m\) es evaluar el out-of-bag-error o recurrir a validación cruzada.

Al igual que ocurre con bagging, random forest no sufre problemas de overfit por aumentar el número de árboles creados en el proceso. Alcanzado un determinado número, la reducción del error de test se estabiliza.

Predicción de Random Forest

La predicción de un modelo Random Forest es la media de las predicciones de todos los árboles que lo forman.

Supóngase que se dispone de 10 observaciones, cada una con un valor de variable respuesta \(Y\) y unos predictores \(X\).

La siguiente imagen muestra cómo sería la predicción del modelo para una nueva observación. En cada árbol, el camino hasta llegar al nodo final está resaltado. En cada nodo terminal se detalla el índice de las observaciones de entrenamiento que forman parte.

Predicción con random forest: en cada árbol, el camino hasta llegar al nodo final está resaltado. En cada nodo terminal se detalla el índice de las observaciones de entrenamiento que forman parte de él.

El valor predicho por cada árbol es la media de la variable respuesta \(Y\) en el nodo terminal. Acorde a la imagen, las predicciones de cada uno de los tres árboles (de izquierda a derecha) es:

\[\hat{y}_{arbol_1} = \frac{24 + 2 + 20}{3} = 15.33333\]

\[\hat{y}_{arbol_2} = \frac{18 + 24}{2} = 21\]

\[\hat{y}_{arbol_3} = \frac{18 + 24 + 2 + 20}{4} = 16\]

La predicción final del modelo es la media de todas las predicciones individuales:

\[\hat{\mu} = \frac{15.33333+21+16}{3} = 17.4\]

Aunque la anterior es la forma más común de obtener las predicciones de un modelo Random Forest, existe otra aproximación. La predicción de un árbol de regresión puede verse como una variante de vecinos cercanos en la que, solo las observaciones que forman parte del mismo nodo terminal que la observación predicha, tienen influencia. Siguiendo esta aproximación, la predicción del árbol se define como la media ponderada de todas las observaciones de entrenamiento, donde el peso de cada observación depende únicamente de si forma parte o no del mismo nodo terminal. Para Random Forest esto equivale a la media ponderada de todas las observaciones, empleando como pesos \(\textbf{w}\) la media de los vectores de pesos de todos los árboles.

Acorde a la imagen anterior, el vector de pesos para cada uno de los tres árboles (de izquierda a derecha) es:

\[\textbf{w}_{arbol_1} = (0, 0, \frac{1}{3}, 0, \frac{1}{3}, 0, 0, 0, \frac{1}{3}, 0)\]

\[\textbf{w}_{arbol_2} = (0, \frac{1}{2}, \frac{1}{2}, 0, 0, 0, 0, 0, 0, 0)\] \[\textbf{w}_{arbol_3} = (0, \frac{1}{4}, \frac{1}{4}, 0, \frac{1}{4}, 0, 0, 0, \frac{1}{4}, 0)\]

La media de todos los vectores de pesos es:

\[\overline{\textbf{w}} = \frac{1}{3}(\textbf{w}_{arbol_1} + \textbf{w}_{arbol_2} + \textbf{w}_{arbol_3}) = \\ (0, \frac{1}{4}, \frac{13}{36}, 0, \frac{7}{36}, 0, 0, 0, \frac{7}{36}, 0)\]

Una vez obtenido el vector de pesos promedio, se puede calcular la predicción con la media ponderada de todas las observaciones de entrenamiento:

\[\hat{\mu} = \sum^n_{i =1} \overline{\textbf{w}}_i Y_i\]

\[\hat{\mu} = (0\times10) + (\frac{1}{4}\times18) + (\frac{13}{36}\times24) + (0\times8) + (\frac{1}{4}\times2) + \\ (0\times9) + (0\times16) + (0\times10) + (\frac{1}{4}\times20) + (0\times14) = 17.4\]

Out-of-Bag Error

Dada la naturaleza del proceso de bagging, resulta posible estimar el error de test sin necesidad de recurrir a métodos de validación cruzada (cross-validation). El hecho de que los árboles se ajusten empleando muestras generadas por bootstrapping conlleva que, en promedio, cada ajuste use solo aproximadamente dos tercios de las observaciones originales. Al tercio restante se le llama out-of-bag (OOB).

Si para cada árbol ajustado en el proceso de bagging se registran las observaciones empleadas, se puede predecir la respuesta de la observación \(i\) haciendo uso de aquellos árboles en los que esa observación ha sido excluida y promediándolos (la moda en el caso de los árboles de clasificación). Siguiendo este proceso, se pueden obtener las predicciones para las \(n\) observaciones y con ellas calcular el OOB-mean square error (para regresión) o el OOB-classification error (para árboles de clasificación). Como la variable respuesta de cada observación se predice empleando únicamente los árboles en cuyo ajuste no participó dicha observación, el OOB-error sirve como estimación del error de test. De hecho, si el número de árboles es suficientemente alto, el OOB-error es prácticamente equivalente al leave-one-out cross-validation error.

Esta es una ventaja añadida de los métodos de bagging, y por lo tanto de Random Forest ya que evita tener que recurrir al proceso de validación cruzada (computacionalmente costoso) para la optimización de los hiperparámetros.

Dos limitaciones en el uso Out-of-Bag Error:

El Out-of-Bag Error no es adecuado cuando las observaciones tienen una relación temporal (series temporales). Como la selección de las observaciones que participan en cada entrenamiento es aleatoria, no respetan el orden temporal y se estaría introduciendo información a futuro.
El preprocesado de los datos de entrenamiento se hace de forma conjunta, por lo que las observaciones out-of-bag pueden sufrir data leakage. De ser así, las estimaciones del OOB-error son demasiado optimistas.

En un muestreo por bootstrapping, si el tamaño de los datos de entrenamiento es \(n\), cada observación tiene una probabilidad de ser elegida de \(\frac{1}{n}\). Por lo tanto, la probabilidad de no ser elegida en todo el proceso es de \((1-1/n)n\), lo que converge en \(1/\epsilon\), que es aproximadamente un tercio.

Importancia de los predictores

Si bien es cierto que el proceso de bagging (Random Forest) consigue mejorar la capacidad predictiva en comparación a los modelos basados en un único árbol, esto tiene un coste asociado, la interpretabilidad del modelo se reduce. Al tratarse de una combinación de múltiples árboles, no es posible obtener una representación gráfica sencilla del modelo y no es inmediato identificar de forma visual qué predictores son más importantes. Sin embargo, se han desarrollado nuevas estrategias para cuantificar la importancia de los predictores que hacen de los modelos de bagging (Random Forest) una herramienta muy potente, no solo para predecir, sino también para el análisis exploratorio. Dos de estas medidas son: importancia por permutación e impureza de nodos.

Importancia por permutación

Identifica la influencia que tiene cada predictor sobre una determinada métrica de evaluación del modelo (estimada por out-of-bag error o validación cruzada). El valor asociado con cada predictor se obtiene de la siguiente forma:

Crear el conjunto de árboles que forman el modelo.
Calcular una determinada métrica de error (mse, classification error, …). Este es el valor de referencia (\(error_0\)).
Para cada predictor \(j\):
- Permutar en todos los árboles del modelo los valores del predictor \(j\) manteniendo el resto constante.
- Recalcular la métrica tras la permutación, llámese (\(error_j\)).
- Calcular el incremento en la métrica debido a la permutación del predictor \(j\).

\[\text{%Incremento}_j = (error-j - error_0)/ error_0 * 100%\]

Si el predictor permutado estaba contribuyendo al modelo, es de esperar que el modelo aumente su error, ya que se pierde la información que proporcionaba esa variable. El porcentaje en que se incrementa el error debido a la permutación del predictor \(j\) puede interpretarse como la influencia que tiene \(j\) sobre el modelo. Algo que suele llevar a confusiones es el hecho de que este incremento puede resultar negativo. Si la variable no contribuye al modelo, es posible que, al reorganizarla aleatoriamente, solo por azar, se consiga mejorar ligeramente el modelo, por lo que \((error_j - error_0)\) es negativo. A modo general, se puede considerar que estas variables tienen una importancia próxima a cero.

Aunque esta estrategia suele ser la más recomendada, cabe tomar algunas precauciones en su interpretación. Lo que cuantifican es la influencia que tienen los predictores sobre el modelo, no su relación con la variable respuesta. ¿Por qué es esto tan importante? Supóngase un escenario en el que se emplea esta estrategia con la finalidad de identificar qué predictores están relacionados con el peso de una persona, y que dos de los predictores son: el índice de masa corporal (IMC) y la altura. Como IMC y altura están muy correlacionados entre sí (la información que aportan es redundante), cuando se permute uno de ellos, el impacto en el modelo será mínimo, ya que el otro aporta la misma información. Como resultado, estos predictores aparecerán como poco influyentes aun cuando realmente están muy relacionados con la variable respuesta. Una forma de evitar problemas de este tipo es, siempre que se excluyan predictores de un modelo, comprobar el impacto que tiene en su capacidad predictiva.

Incremento de la pureza de nodos

Cuantifica el incremento total en la pureza de los nodos debido a divisiones en las que participa el predictor (promedio de todos los árboles). La forma de calcularlo es la siguiente: en cada división de los árboles, se registra el descenso conseguido en la medida empleada como criterio de división (índice Gini, mse entropía, …). Para cada uno de los predictores, se calcula el descenso medio conseguido en el conjunto de árboles que forman el ensemble. Cuanto mayor sea este valor medio, mayor la contribución del predictor en el modelo.

Ejemplo regresión

Librerías

Las librerías utilizadas en este ejemplo son:

# Tratamiento de datos
# ==============================================================================
library(MASS)
library(dplyr)
library(tidyr)
library(skimr)

# Gráficos
# ==============================================================================
library(ggplot2)
library(ggpubr)

# Preprocesado y modelado
# ==============================================================================
library(tidymodels)
library(ranger)
library(doParallel)

Datos

El set de datos Boston disponible en el paquete MASS contiene precios de viviendas de la ciudad de Boston, así como información socio-económica del barrio en el que se encuentran. Se pretende ajustar un modelo de regresión que permita predecir el precio medio de una vivienda (medv) en función de las variables disponibles.

data("Boston")
datos = Boston
head(datos, 3)

skim(datos)

Data summary
Name	datos
Number of rows	506
Number of columns	14
_______________________
Column type frequency:
numeric	14
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
crim	1	3.61	8.60	0.01	0.08	0.26	3.68	88.98	▇▁▁▁▁
zn	1	11.36	23.32	0.00	0.00	0.00	12.50	100.00	▇▁▁▁▁
indus	1	11.14	6.86	0.46	5.19	9.69	18.10	27.74	▇▆▁▇▁
chas	1	0.07	0.25	0.00	0.00	0.00	0.00	1.00	▇▁▁▁▁
nox	1	0.55	0.12	0.38	0.45	0.54	0.62	0.87	▇▇▆▅▁
rm	1	6.28	0.70	3.56	5.89	6.21	6.62	8.78	▁▂▇▂▁
age	1	68.57	28.15	2.90	45.02	77.50	94.07	100.00	▂▂▂▃▇
dis	1	3.80	2.11	1.13	2.10	3.21	5.19	12.13	▇▅▂▁▁
rad	1	9.55	8.71	1.00	4.00	5.00	24.00	24.00	▇▂▁▁▃
tax	1	408.24	168.54	187.00	279.00	330.00	666.00	711.00	▇▇▃▁▇
ptratio	1	18.46	2.16	12.60	17.40	19.05	20.20	22.00	▁▃▅▅▇
black	1	356.67	91.29	0.32	375.38	391.44	396.22	396.90	▁▁▁▁▇
lstat	1	12.65	7.14	1.73	6.95	11.36	16.96	37.97	▇▇▅▂▁
medv	1	22.53	9.20	5.00	17.02	21.20	25.00	50.00	▂▇▅▁▁

Ajuste del modelo

Se ajusta un modelo empleando como variable respuesta medv y como predictores todas las otras variables disponibles.

La función ranger del paquete ranger permite entrenar modelos random forest para problemas de regresión. Los parámetros e hiperparámetros empleados por defecto son:

formula = NULL
data = NULL
num.trees = 500
mtry = NULL
importance = "none"
write.forest = TRUE
probability = FALSE
min.node.size = NULL
max.depth = NULL
replace = TRUE
sample.fraction = ifelse(replace, 1, 0.632)
case.weights = NULL
class.weights = NULL
splitrule = NULL
num.random.splits = 1
alpha = 0.5
minprop = 0.1
split.select.weights = NULL
always.split.variables = NULL
respect.unordered.factors = NULL
scale.permutation.importance = FALSE
local.importance = FALSE
regularization.factor = 1
regularization.usedepth = FALSE
keep.inbag = FALSE
inbag = NULL
holdout = FALSE
quantreg = FALSE
oob.error = TRUE
num.threads = NULL
save.memory = FALSE
verbose = TRUE
seed = NULL
dependent.variable.name = NULL
status.variable.name = NULL
classification = NULL
x = NULL
y = NULL

De entre todos ellos, destacan aquellos que detienen el crecimiento de los árboles, los que controlan el número de árboles y predictores incluidos, y los que gestionan la paralelización:

num.trees; número de árboles incluidos en el modelo.
max.depth: profundidad máxima que pueden alcanzar los árboles. Por defecto (0) crece al máximo los árboles.
min.node.size: número mínimo de observaciones que debe de tener cada uno de los nodos hijos para que se produzca la división. Por defecto emplea 1 para problemas de clasificación y 10 para problemas de regresión.
mtry: número de predictores considerados en cada división. Por defecto, se emplea como valor la raíz cuadrada del número total de predictores disponible (redondeado a la baja).
oob.error: Si se calcula o no el out-of-bag error (accuracy, mean squared error o \(R^2\)). Por defecto es FALSE ya que aumenta el tiempo de entrenamiento.
num.threads: número de cores empleados para el entrenamiento. En random forest los árboles se ajustan de forma independiente, por lo que la paralelización reduce notablemente el tiempo de entrenamiento. Por defecto se utilizan todos los cores disponibles.
seed: semilla para que los resultados sean reproducibles.

Como en todo estudio predictivo, no solo es importante ajustar el modelo, sino también cuantificar su capacidad para predecir nuevas observaciones. Para poder hacer esta evaluación, se dividen los datos en dos grupos, uno de entrenamiento y otro de test

# División de los datos en train y test
# ==============================================================================
set.seed(123)
train       <- sample(1:nrow(datos), size = nrow(datos)/2)
datos_train <- datos[train,]
datos_test  <- datos[-train,]

# Creación y entrenamiento del modelo
# ==============================================================================
set.seed(123)
modelo  <- ranger(
            formula   = medv ~ .,
            data      = datos_train,
            num.trees = 10,
            seed      = 123
           )

print(modelo)

## Ranger result
## 
## Call:
##  ranger(formula = medv ~ ., data = datos_train, num.trees = 10,      seed = 123) 
## 
## Type:                             Regression 
## Number of trees:                  10 
## Sample size:                      253 
## Number of independent variables:  13 
## Mtry:                             3 
## Target node size:                 5 
## Variable importance mode:         none 
## Splitrule:                        variance 
## OOB prediction error (MSE):       13.98202 
## R squared (OOB):                  0.828173

Predicción y evaluación del modelo

Una vez entrenado el modelo, se evalúa la capacidad predictiva empleando el conjunto de test.

# Error de test del modelo
# ==============================================================================
predicciones <- predict(
                  modelo,
                  data = datos_test
                )

predicciones <- predicciones$predictions
test_rmse    <- sqrt(mean((predicciones - datos_test$medv)^2))
paste("Error de test (rmse) del modelo:", round(test_rmse,2))

## [1] "Error de test (rmse) del modelo: 4.15"

Optimización de hiperparámetros

El modelo inicial se ha entrenado utilizando 10 árboles (num.trees=10) y manteniendo el resto de hiperparámetros con su valor por defecto. Al ser hiperparámetros, no se puede saber de antemano cuál es el valor más adecuado, la forma de identificarlos es mediante el uso de estrategias de validación, por ejemplo validación cruzada.

Los modelos Random Forest tienen la ventaja de disponer del Out-of-Bag error, lo que permite obtener una estimación del error de test sin recurrir a la validación cruzada, que es computacionalmente costosa. En la implementación de ranger, para problemas de regresión, se calculan dos métricas como oob_score: el mean squared error y el \(R^2\). Si se desea otra, se tiene que recurrir a validación cruzada.

Cabe tener en cuenta que, cuando se busca el valor óptimo de un hiperparámetro con dos métricas distintas, el resultado obtenido raramente es el mismo. Lo importante es que ambas métricas identifiquen las mismas regiones de interés. A continuación, se muestran las dos aproximaciones.

Número de árboles

En Random Forest, el número de árboles no es un hiperparámetro crítico en cuanto que, añadir árboles, solo puede hacer que mejorar el resultado. En Random Forest no se produce overfitting por exceso de árboles sin embargo, añadir árboles una vez que la mejora se estabiliza es una pérdida de recursos computacionales.

# Validación empleando el Out-of-Bag error (root mean squared error)
# ==============================================================================

# Valores evaluados
num_trees_range <- seq(1, 400, 20)

# Bucle para entrenar un modelo con cada valor de num_trees y extraer su error
# de entrenamiento y de Out-of-Bag.

train_errors <- rep(NA, times = length(num_trees_range))
oob_errors   <- rep(NA, times = length(num_trees_range))

for (i in seq_along(num_trees_range)){
  modelo  <- ranger(
               formula   = medv ~ .,
               data      = datos_train,
               num.trees = num_trees_range[i],
               oob.error = TRUE,
               seed      = 123
             )
  
  predicciones_train <- predict(
                          modelo,
                          data = datos_train
                        )
  predicciones_train <- predicciones_train$predictions
  
  train_error <- mean((predicciones_train - datos_train$medv)^2)
  oob_error   <- modelo$prediction.error
  
  train_errors[i] <- sqrt(train_error)
  oob_errors[i]   <- sqrt(oob_error)
  
}

# Gráfico con la evolución de los errores
df_resulados <- data.frame(n_arboles = num_trees_range, train_errors, oob_errors)
ggplot(data = df_resulados) +
  geom_line(aes(x = num_trees_range, y = train_errors, color = "train rmse")) + 
  geom_line(aes(x = num_trees_range, y = oob_errors, color = "oob rmse")) +
  geom_vline(xintercept = num_trees_range[which.min(oob_errors)],
             color = "firebrick",
             linetype = "dashed") +
  labs(
    title = "Evolución del out-of-bag-error vs número árboles",
    x     = "número de árboles",
    y     = "out-of-bag-error (rmse)",
    color = ""
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

paste("Valor óptimo de num.trees:", num_trees_range[which.min(oob_errors)])

## [1] "Valor óptimo de num.trees: 281"

El paquete ranger no tiene ninguna funcionalidad propia para realizar validación cruzada, por lo que se recurre a tidymodels. Puede encontrarse una descripción más detallada de este proceso en Machine Learning con R y tidymodels. Otro framework excelente para machine learning en R es mlr3.

# Validación empleando k-cross-validation (root mean squared error)
# ==============================================================================

# Valores evaluados
num_trees_range <- seq(1, 400, 10)

# Bucle para entrenar un modelo con cada valor de num_trees y extraer su error
# de entrenamiento y de Out-of-Bag.

train_errors <- rep(NA, times = length(num_trees_range))
cv_errors    <- rep(NA, times = length(num_trees_range))

for (i in seq_along(num_trees_range)){
  
  # Definición del modelo
  modelo <- rand_forest(
              mode  = "regression",
              trees = num_trees_range[i]
            ) %>%
            set_engine(
              engine = "ranger",
              seed   = 123
            )
  
  # Particiones validación cruzada
  set.seed(1234)
  cv_folds <- vfold_cv(
                data    = datos_train,
                v       = 5,
                repeats = 1
              )
  
  # Ejecución validación cruzada
  validacion_fit <- fit_resamples(
                      preprocessor = medv ~ .,
                      object       = modelo,
                      resamples    = cv_folds,
                      metrics      = metric_set(rmse)
                    )
  
  # Extraer la métrica de validación 
  cv_error <- collect_metrics(validacion_fit)$mean
  
  # Predicción datos train
  modelo_fit <- modelo %>% fit(medv ~ ., data = datos_train)
  predicciones_train <- predict(
                          modelo_fit,
                          new_data = datos_train
                        )
  predicciones_train <- predicciones_train$.pred
  
  train_error <- sqrt(mean((predicciones_train - datos_train$medv)^2))
  
  # Resultados
  train_errors[i] <- train_error
  cv_errors[i]    <- cv_error
  
}

# Gráfico con la evolución de los errores
df_resulados <- data.frame(n_arboles = num_trees_range, train_errors, cv_errors)
ggplot(data = df_resulados) +
  geom_line(aes(x = num_trees_range, y = train_errors, color = "train rmse")) + 
  geom_line(aes(x = num_trees_range, y = cv_errors, color = "cv rmse")) +
  geom_vline(xintercept = num_trees_range[which.min(oob_errors)],
             color = "firebrick",
             linetype = "dashed") +
  labs(
    title = "Evolución del cv-error vs número árboles",
    x     = "número de árboles",
    y     = "cv-error (rmse)",
    color = ""
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

paste("Valor óptimo de num.trees:", num_trees_range[which.min(cv_errors)])

## [1] "Valor óptimo de num.trees: 41"

Ambas métricas indican que, a partir de entre 50 y 100 árboles, el error de validación del modelo se estabiliza.

mtry

El valor de mtry es uno de los hiperparámetros más importantes de random forest, ya que es el que permite controlar cuánto se decorrelacionan los árboles entre sí.

# Validación empleando el Out-of-Bag error (root mean squared error)
# ==============================================================================

# Valores evaluados
mtry_range <- seq(1, ncol(datos_train)-1)

# Bucle para entrenar un modelo con cada valor de mtry y extraer su error
# de entrenamiento y de Out-of-Bag.

train_errors <- rep(NA, times = length(mtry_range))
oob_errors   <- rep(NA, times = length(mtry_range))

for (i in seq_along(mtry_range)){
  modelo  <- ranger(
               formula   = medv ~ .,
               data      = datos_train,
               num.trees = 50,
               mtry      = mtry_range[i],
               oob.error = TRUE,
               seed      = 123
             )
  
  predicciones_train <- predict(
                          modelo,
                          data = datos_train
                        )
  predicciones_train <- predicciones_train$predictions
  
  train_error <- mean((predicciones_train - datos_train$medv)^2)
  oob_error   <- modelo$prediction.error
  
  train_errors[i] <- sqrt(train_error)
  oob_errors[i]   <- sqrt(oob_error)
  
}

# Gráfico con la evolución de los errores
df_resulados <- data.frame(mtry = mtry_range, train_errors, oob_errors)
ggplot(data = df_resulados) +
  geom_line(aes(x = mtry_range, y = train_errors, color = "train rmse")) + 
  geom_line(aes(x = mtry_range, y = oob_errors, color = "oob rmse")) +
  geom_vline(xintercept =  mtry_range[which.min(oob_errors)],
             color = "firebrick",
             linetype = "dashed") +
  labs(
    title = "Evolución del out-of-bag-error vs mtry",
    x     = "mtry",
    y     = "out-of-bag-error (rmse)",
    color = ""
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

paste("Valor óptimo de mtry:", mtry_range[which.min(oob_errors)])

## [1] "Valor óptimo de mtry: 10"

# Validación empleando k-cross-validation (root mean squared error)
# ==============================================================================

# Valores evaluados
mtry_range <- seq(1, ncol(datos_train)-1)

# Bucle para entrenar un modelo con cada valor de mtry y extraer su error
# de entrenamiento y de Out-of-Bag.

train_errors <- rep(NA, times = length(mtry_range))
cv_errors    <- rep(NA, times = length(mtry_range))

for (i in seq_along(mtry_range)){
  
  # Definición del modelo
  modelo <- rand_forest(
              mode  = "regression",
              trees = 100,
              mtry  = mtry_range[i]
            ) %>%
            set_engine(
              engine = "ranger",
              seed   = 123
            )
  
  # Particiones validación cruzada
  set.seed(1234)
  cv_folds <- vfold_cv(
                data    = datos_train,
                v       = 5,
                repeats = 1
              )
  
  # Ejecución validación cruzada
  validacion_fit <- fit_resamples(
                      preprocessor = medv ~ .,
                      object       = modelo,
                      resamples    = cv_folds,
                      metrics      = metric_set(rmse)
                    )
  
  # Extraer datos de validación
  cv_error <- collect_metrics(validacion_fit)$mean
  
  # Predicción datos train
  modelo_fit <- modelo %>% fit(medv ~ ., data = datos_train)
  predicciones_train <- predict(
                          modelo_fit,
                          new_data = datos_train
                        )
  predicciones_train <- predicciones_train$.pred
  
  train_error <- sqrt(mean((predicciones_train - datos_train$medv)^2))
  
  # Resultados
  train_errors[i] <- train_error
  cv_errors[i]    <- cv_error
  
}

# Gráfico con la evolución de los errores
df_resulados <- data.frame(mtry = mtry_range, train_errors, cv_errors)
ggplot(data = df_resulados) +
  geom_line(aes(x = mtry_range, y = train_errors, color = "train error")) + 
  geom_line(aes(x = mtry_range, y = cv_errors, color = "cv error")) +
  geom_vline(xintercept =  mtry_range[which.min(cv_errors)],
             color = "firebrick",
             linetype = "dashed") +
  labs(
    title = "Evolución del out-of-bag-error vs mtry",
    x     = "mtry",
    y     = "cv-error (mse)",
    color = ""
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

paste("Valor óptimo de mtry:", mtry_range[which.min(cv_errors)])

## [1] "Valor óptimo de mtry: 7"

Acorde a las dos métricas utilizadas, el valor óptimo de mtry está entre 7 y 10.

Grid search

Aunque el análisis individual de los hiperparámetros es útil para entender su impacto en el modelo e identificar rangos de interés, la búsqueda final no debe hacerse de forma secuencial, ya que cada hiperparámetro interacciona con los demás. Es preferible recurrir a grid search o random search para analizar varias combinaciones de hiperparámetros. Puede encontrarse más información sobre las estrategias de búsqueda en Tidymodels y en mlr3.

Grid Search basado en out-of-bag error

# Grid de hiperparámetros evaluados
# ==============================================================================
param_grid = expand_grid(
                'num_trees' = c(50, 100, 500, 1000, 5000),
                'mtry'      = c(3, 5, 7, ncol(datos_train)-1),
                'max_depth' = c(1, 3, 10, 20)
                
             )

# Loop para ajustar un modelo con cada combinación de hiperparámetros
# ==============================================================================

oob_error = rep(NA, nrow(param_grid))

for(i in 1:nrow(param_grid)){

  modelo <- ranger(
              formula   = medv ~ .,
              data      = datos_train, 
              num.trees = param_grid$num_trees[i],
              mtry      = param_grid$mtry[i],
              max.depth = param_grid$max_depth[i],
              seed      = 123
            )

 oob_error[i] <- sqrt(modelo$prediction.error)
}


# Resultados
# ==============================================================================
resultados <- param_grid
resultados$oob_error <- oob_error
resultados <- resultados %>% arrange(oob_error)

# Mejores hiperparámetros por out-of-bag error
# ==============================================================================
head(resultados, 1)

Grid Search basado en validación cruzada

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# ==============================================================================
modelo <- rand_forest(
             mode  = "regression",
             mtry  = tune(),
             trees = tune()
          ) %>%
          set_engine(
            engine     = "ranger",
            max.depth  = tune(),
            importance = "none",
            seed       = 123
          )

# DEFINICIÓN DEL PREPROCESADO
# ==============================================================================
# En este caso no hay preprocesado, por lo que el transformer solo contiene
# la definición de la fórmula y los datos de entrenamiento.
transformer <- recipe(
                  formula = medv ~ .,
                  data    =  datos_train
               )

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# ==============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = medv
            )

# WORKFLOW
# ==============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo)
                     

# GRID DE HIPERPARÁMETROS
# ==============================================================================
hiperpar_grid <- expand_grid(
                  'trees'     = c(50, 100, 500, 1000, 5000),
                  'mtry'      = c(3, 5, 7, ncol(datos_train)-1),
                  'max.depth' = c(1, 3, 10, 20)
                 )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# ==============================================================================
cl <- makePSOCKcluster(parallel::detectCores() - 1)
registerDoParallel(cl)

grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              grid      = hiperpar_grid
            )

stopCluster(cl)

# Mejores hiperparámetros por validación cruzada
# ==============================================================================
show_best(grid_fit, metric = "rmse", n = 1)

Una vez identificados los mejores hiperparámetros, se reentrena el modelo indicando los valores óptimos en sus argumentos.

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_final_fit <- finalize_workflow(
                        x = workflow_modelado,
                        parameters = mejores_hiperpar
                    ) %>%
                    fit(
                      data = datos_train
                    ) %>%
                    pull_workflow_fit()

# Error de test del modelo final
# ==============================================================================
predicciones <- modelo_final_fit %>%
                predict(
                  new_data = datos_test,
                  type     = "numeric"
                )

predicciones <- predicciones %>% 
                bind_cols(datos_test %>% dplyr::select(medv))

rmse_test  <- rmse(
                 data     = predicciones,
                 truth    = medv,
                 estimate = .pred,
                 na_rm    = TRUE
              )
rmse_test

Tras optimizar los hiperparámetros, se consigue reducir el error rmse del modelo de 4.15 a 3.57. Las predicciones del modelo final se alejan en promedio 3.57 unidades (3570 dólares) del valor real.

Importancia de predictores

Importancia por pureza de nodos

En los modelos anteriores, el argumento importance se deja por defecto como "none". Esto desactiva el cálculo de importancia de predictores para reducir así el tiempo de entrenamiento. Se entrena de nuevo el modelo, con los mejores hiperparámetros encontrados, pero esta vez indicando importance = "impurity". Los modelos ranger calculan la impureza a partir del índice \(Gini\) en problemas de clasificación y con la varianza en regresión.

# Entrenamiento modelo
modelo <- rand_forest(
             mode  = "regression"
          ) %>%
          set_engine(
            engine     = "ranger",
            importance = "impurity",
            seed       = 123
          )

modelo <- modelo %>% finalize_model(mejores_hiperpar)
modelo <- modelo %>% fit(medv ~., data = datos_train)

# Importancia
importancia_pred <- modelo$fit$variable.importance %>%
                    enframe(name = "predictor", value = "importancia")

# Gráfico
ggplot(
  data = importancia_pred,
  aes(x    = reorder(predictor, importancia),
      y    = importancia,
      fill = importancia)
) +
labs(x = "predictor", title = "Importancia predictores (pureza de nodos)") +
geom_col() +
coord_flip() +
theme_bw() +
theme(legend.position = "none")

Importancia por permutación

Se entrena de nuevo el modelo, con los mejores hiperparámetros encontrados, pero esta vez indicando importance = "permutation".

# Entrenamiento modelo
modelo <- rand_forest(
             mode  = "regression"
          ) %>%
          set_engine(
            engine     = "ranger",
            importance = "permutation",
            seed       = 123
          )

modelo <- modelo %>% finalize_model(mejores_hiperpar)
modelo <- modelo %>% fit(medv ~., data = datos_train)

# Importancia
importancia_pred <- modelo$fit$variable.importance %>%
                    enframe(name = "predictor", value = "importancia")

# Gráfico
ggplot(
  data = importancia_pred,
  aes(x    = reorder(predictor, importancia),
      y    = importancia,
      fill = importancia)
) +
labs(x = "predictor", title = "Importancia predictores (permutación)") +
geom_col() +
coord_flip() +
theme_bw() +
theme(legend.position = "none")

Ambas estrategias identifican lstats y rm como los predictores más influyentes, acorde a los datos de entrenamiento.

Ejemplo clasificación

Librerías

Las librerías utilizadas en este ejemplo son:

# Tratamiento de datos
# ==============================================================================
library(ISLR)
library(dplyr)
library(tidyr)
library(skimr)

# Gráficos
# ==============================================================================
library(ggplot2)
library(ggpubr)

# Preprocesado y modelado
# ==============================================================================
library(tidymodels)
library(ranger)
library(doParallel)

Datos

data("Carseats")
datos = Carseats
head(datos, 3)

Como Sales es una variable continua y el objetivo del estudio es clasificar las tiendas según si venden mucho o poco, se crea una nueva variable dicotómica (Si, No) llamada ventas_altas.

datos$ventas_altas <- ifelse(test = datos$Sales > 8, 
                             yes = "Si", no = "No")

# Conversión de la variable respuesta a tipo factor
datos$ventas_altas <- as.factor(datos$ventas_altas)

# Una vez creada la nueva variable respuesta se descarta la original
datos$Sales = NULL

Ajuste del modelo y optimización de hiperparámetros

Se ajusta un árbol de clasificación empleando como variable respuesta ventas_altas y como predictores todas las variables disponibles. A diferencia del ejemplo de regresión, en estos datos hay variables categóricas. ranger es capaz de utilizar este tipo de variables, pero otras implementaciones requieren aplicar one-hot-encoding.

# División de los datos en train y test
# ==============================================================================
set.seed(123)
train       <- sample(1:nrow(datos), size = nrow(datos)/2)
datos_train <- datos[train,]
datos_test  <- datos[-train,]

Grid Search basado en out-of-bag error

# Grid de hiperparámetros evaluados
# ==============================================================================
param_grid = expand_grid(
                'num_trees' = c(50, 100, 500, 1000, 5000),
                'mtry'      = c(3, 5, 7, ncol(datos_train)-1),
                'max_depth' = c(1, 3, 10, 20)
                
             )

# Loop para ajustar un modelo con cada combinación de hiperparámetros
# ==============================================================================

oob_error = rep(NA, nrow(param_grid))

for(i in 1:nrow(param_grid)){

  modelo <- ranger(
              formula   = ventas_altas ~ .,
              data      = datos_train, 
              num.trees = param_grid$num_trees[i],
              mtry      = param_grid$mtry[i],
              max.depth = param_grid$max_depth[i],
              seed      = 123
            )

 oob_error[i] <- modelo$prediction.error
}


# Resultados
# ==============================================================================
resultados <- param_grid
resultados$oob_error <- oob_error
resultados <- resultados %>% arrange(oob_error)

# Mejores hiperparámetros por out-of-bag error
# ==============================================================================
head(resultados, 1)

Grid Search basado en validación cruzada

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# ==============================================================================
modelo <- rand_forest(
             mode  = "classification",
             mtry  = tune(),
             trees = tune()
          ) %>%
          set_engine(
            engine     = "ranger",
            max.depth  = tune(),
            importance = "none",
            seed       = 123
          )

# DEFINICIÓN DEL PREPROCESADO
# ==============================================================================
# En este caso no hay preprocesado, por lo que el transformer solo contiene
# la definición de la fórmula y los datos de entrenamiento.
transformer <- recipe(
                  formula = ventas_altas ~ .,
                  data    =  datos_train
               )

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# ==============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = ventas_altas
            )

# WORKFLOW
# ==============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo)
                     

# GRID DE HIPERPARÁMETROS
# ==============================================================================
hiperpar_grid <- expand_grid(
                  'trees'     = c(50, 100, 500, 1000, 5000),
                  'mtry'      = c(3, 5, 7, ncol(datos_train)-1),
                  'max.depth' = c(1, 3, 10, 20)
                 )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# ==============================================================================
cl <- makePSOCKcluster(parallel::detectCores() - 1)
registerDoParallel(cl)

grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(accuracy),
              grid      = hiperpar_grid
            )

stopCluster(cl)

# Mejores hiperparámetros por validación cruzada
# ==============================================================================
show_best(grid_fit, metric = "accuracy", n = 1)

Predicción y evaluación del modelo

Una vez identificados los mejores hiperparámetros, se reentrena el modelo indicando los valores óptimos en sus argumentos.

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "accuracy")

modelo_final_fit <- finalize_workflow(
                        x = workflow_modelado,
                        parameters = mejores_hiperpar
                    ) %>%
                    fit(
                      data = datos_train
                    ) %>%
                    pull_workflow_fit()

# Error de test del modelo final
# ==============================================================================
predicciones <- modelo_final_fit %>%
                predict(new_data = datos_test)

predicciones <- predicciones %>% 
                bind_cols(datos_test %>% dplyr::select(ventas_altas))

accuracy_test  <- accuracy(
                     data     = predicciones,
                     truth    = ventas_altas,
                     estimate = .pred_class,
                     na_rm    = TRUE
                  )
accuracy_test

mat_confusion <- predicciones %>%
                 conf_mat(
                   truth     = ventas_altas,
                   estimate  = .pred_class
                 )
mat_confusion

##           Truth
## Prediction  No  Si
##         No 100  15
##         Si  22  63

Tras optimizar los hiperparámetros, se consigue un porcentaje de acierto del 81.5%.

Predicción de probabilidades

La mayoría de implementaciones de Random Forest, entre ellas la de ranger, permiten predecir probabilidades cuando se trata de problemas de clasificación. Es importante entender cómo se calculan estos valores para interpretarlos y utilizarlos correctamente.

En el ejemplo anterior, al aplicar predict() se devuelve Si (ventas elevadas) o No (ventas bajas) para cada observación de test. Sin embargo, no se dispone de ningún tipo de información sobre la seguridad con la que el modelo realiza esta asignación. Con predict(type="prob"), en lugar de una clasificación, se obtiene la probabilidad con la que el modelo considera que cada observación puede pertenecer a cada una de las clases.

# Predicción de probabilidades
# ==============================================================================
predicciones <- modelo_final_fit %>%
                predict(new_data = datos_test, type = "prob")
head(predicciones, 4)

El resultado de predict(type="prob") es un dataframe con una fila por observación y tantas columnas como clases tenga la variable respuesta. El valor de la primera columna se corresponde con la probabilidad, acorde al modelo, de que la observación pertenezca a la clase No, y así sucesivamente. El valor de probabilidad mostrado para cada predicción se corresponde con la fracción de observaciones de cada clase en los nodos terminales a los que ha llegado la observación predicha en el conjunto de los árboles.

Por defecto, predict() asigna cada nueva observación a la clase con mayor probabilidad (en caso de empate se asigna de forma aleatoria). Sin embargo, este no tiene por qué ser el comportamiento deseado en todos los casos.

# Clasificación empleando la clase de mayor probabilidad
# ==============================================================================
predicciones <- predicciones %>%
                mutate(
                  clasificacion_default_0.5 = if_else(.pred_No > .pred_Si, "No", "Si")
                )
head(predicciones, 4)

Supóngase el siguiente escenario: la campaña de navidad se aproxima y los propietarios de la cadena quieren duplicar el stock de artículos en aquellas tiendas de las que se prevé que tengan ventas elevadas. Como el transporte de este material hasta las tiendas supone un coste elevado, el director quiere limitar esta estrategia únicamente a tiendas para las que se tenga mucha seguridad de que van conseguir muchas ventas.

Si se dispone de las probabilidades, se puede establecer un punto de corte concreto, por ejemplo, considerando únicamente como clase Si (ventas altas) aquellas tiendas cuya predicción para esta clase sea superior al 0.9 (90%). De esta forma, la clasificación final se ajusta mejor a las necesidades del caso de uso.

# Clasificación final empleando un threshold de 0.9 para la clase 1.
# ==============================================================================
predicciones <- predicciones %>%
                mutate(
                  clasificacion_default_0.9 = if_else(.pred_Si > 0.9, "Si", "No")
                )
head(predicciones, 4)

¿Hasta qué punto se debe confiar en estas probabilidades?

Es muy importante tener en cuenta la diferencia entre la “visión” que tiene el modelo del mundo y el mundo real. Todo lo que sabe un modelo sobre el mundo real es lo que ha podido aprender de los datos de entrenamiento y, por lo tanto, tiene una “visión” limitada. Por ejemplo, supóngase que, en los datos de entrenamiento, todas las tiendas que están en zona urbana Urban='Yes' tienen ventas altas independientemente del valor que tomen el resto de predictores. Cuando el modelo trate de predecir una nueva observación, si esta está en zona urbana, clasificará a la tienda como ventas elevadas con un 100% de seguridad. Sin embargo, esto no significa que sea inequivocamente cierto, podría haber tiendas en zonas urbanas que no tienen ventas elevadas pero, al no estar presentes en los datos de entrenamiento, el modelo no contempla esta posibilidad.

Teniendo en cuenta todo esto, hay que considerar las probabilidades generadas por el modelo como la seguridad que tiene este, desde su visión limitada, al realizar las predicciones. Pero no como la probabilidad en el mundo real de que así lo sea.

Importancia de predictores

Importancia por pureza de nodos

# Entrenamiento modelo
modelo <- rand_forest(
             mode  = "classification"
          ) %>%
          set_engine(
            engine     = "ranger",
            importance = "impurity",
            seed       = 123
          )

modelo <- modelo %>% finalize_model(mejores_hiperpar)
modelo <- modelo %>% fit(ventas_altas ~., data = datos_train)

# Importancia
importancia_pred <- modelo$fit$variable.importance %>%
                    enframe(name = "predictor", value = "importancia")

# Gráfico
ggplot(
  data = importancia_pred,
  aes(x    = reorder(predictor, importancia),
      y    = importancia,
      fill = importancia)
) +
labs(x = "predictor", title = "Importancia predictores (pureza de nodos)") +
geom_col() +
scale_fill_viridis_c() +
coord_flip() +
theme_bw() +
theme(legend.position = "none")

Importancia por permutación

Se entrena de nuevo el modelo, con los mejores hiperparámetros encontrados, pero esta vez indicando importance = "permutation".

# Entrenamiento modelo
modelo <- rand_forest(
             mode  = "classification"
          ) %>%
          set_engine(
            engine     = "ranger",
            importance = "permutation",
            seed       = 123
          )

modelo <- modelo %>% finalize_model(mejores_hiperpar)
modelo <- modelo %>% fit(ventas_altas ~., data = datos_train)

# Importancia
importancia_pred <- modelo$fit$variable.importance %>%
                    enframe(name = "predictor", value = "importancia")

# Gráfico
ggplot(
  data = importancia_pred,
  aes(x    = reorder(predictor, importancia),
      y    = importancia,
      fill = importancia)
) +
labs(x = "predictor", title = "Importancia predictores (permutación)") +
geom_col() +
scale_fill_viridis_c() +
coord_flip() +
theme_bw() +
theme(legend.position = "none")

Ambas estrategias identifican Price, Age y ShelveLoc como los predictores más influyentes, acorde a los datos de entrenamiento.

Boosting

Idea intuitiva

Boosting es otra estrategia de ensemble que se puede emplear con un amplio grupo de métodos de statistical learning, entre ellos los árboles de decisión. La idea detrás del boosting es ajustar, de forma secuencial, múltiples weak learners (modelos sencillos que predicen solo ligeramente mejor que lo esperado por azar). Cada nuevo modelo emplea información del modelo anterior para aprender de sus errores, mejorando iteración a iteración. En el caso de los árboles de predicción, un weak learners se consigue utilizando árboles con muy pocas ramificaciones. A diferencia del método de bagging (random forest), el boosting no hace uso de muestreo repetido (bootstrapping), la diferencia entre los árboles que forman el ensemble se origina por que la importancia (peso) de las observaciones va cambiando en cada iteración.

Tres de los algoritmos de boosting más empleados son AdaBoost, Gradient Boosting y Stochastic Gradient Boosting. Todos ellos se caracterizan por tener una cantidad considerable de hiperparámetros, cuyo valor óptimo se tiene que identificar mediante validación cruzada. Tres de los más importantes son:

El número de weak learners o número de iteraciones: a diferencia del bagging y random forest, el boosting puede sufrir overfitting si este valor es excesivamente alto. Para evitarlo, se emplea un término de regularización conocido como learning rate.
Learning rate: controla la influencia que tiene cada weak learner en el conjunto del ensemble, es decir, el ritmo al que aprende el modelo. Suelen recomendarse valores de 0.001 o 0.01, aunque la elección correcta puede variar dependiendo del problema. Cuanto menor sea su valor, más árboles se necesitan para alcanzar buenos resultados pero menor es el riesgo de overfitting.
Si los weak learners son árboles, el tamaño máximo permitido de cada árbol. Suelen emplearse valores pequeños, entre 1 y 10.

AdaBoost

La publicación en 1995 del algoritmo AdaBoost (Adaptative Boosting) por parte de Yoav Freund y Robert Schapire supuso un avance muy importante en el campo del aprendizaje estadístico, ya que hizo posible aplicar la estrategia de boosting a multitud de problemas. Para el funcionamiento de AdaBoost (problema de clasificación binaria) es necesario establecer:

Un tipo de modelo, normalmente llamado weak learner o base learner, que sea capaz de predecir la variable respuesta con un porcentaje de acierto ligeramente superior a lo esperado por azar. En el caso de los árboles de regresión, este weak learner suele ser un árbol con apenas unos pocos nodos.
Codificar las dos clases de la variable respuesta como +1 y -1.
Un peso inicial e igual para todas las observaciones que forman el set de entrenamiento.

Una vez que estos tres puntos se han establecido, se inicia un proceso iterativo. En la primera iteración, se ajusta el weak learner empleando los datos de entrenamiento y los pesos iniciales (todos iguales). Con el weak learner ajustado y almacenado, se predicen las observaciones de entrenamiento y se identifican aquellas bien y mal clasificadas. Con esta información:

Se actualizan los pesos de las observaciones, disminuyendo el de las que están bien clasificadas y aumentando el de las mal clasificadas.
Se asigna un peso total al weak learner, proporcional al total de aciertos. Cuantos más aciertos consiga el weak learner, mayor su influencia en el conjunto del ensemble.

En la siguiente iteración, se llama de nuevo al weak learner y se vuelve a ajustar, esta vez, empleando los pesos actualizados en la iteración anterior. El nuevo weak learner se almacena, obteniendo así un nuevo modelo para el ensemble. Este proceso se repite M veces, generando un total de M weak learner. Para clasificar nuevas observaciones, se obtiene la predicción de cada uno de los weak learners que forman el ensemble y se agregan sus resultados, ponderando el peso de cada uno acorde al peso que se le ha asignado en el ajuste. El objetivo detrás de esta estrategia es que cada nuevo weak learner se centra en predecir correctamente las observaciones que los anteriores no han sido capaces.

Algoritmo AdaBoost

\(n\): número de observaciones de entrenamiento

\(M\): número de iteraciones de aprendizaje (número total de weak learners)

\(G_m\): weak learner de la iteración \(m\)

\(w_i\): peso de la observación \(i\)

\(\alpha_m\): peso del weak learner \(m\)

Inicializar los pesos de las observaciones, asignando el mismo valor a todas ellas: \[w_i = \frac{1}{N}, i = 1, 2,..., N\]
Para \(m=1\) hasta \(M\):
- 2.1 Ajustar el weak learner \(G_m\) utilizando las observaciones de entrenamiento y los pesos \(w_i\).
- 2.2 Calcular el error del weak learner como: \[err_m = \frac{\sum^N_{i=1} w_iI(y_i \neq G_m(x_i))}{\sum^N_{i=1}w_i}\]
- 2.3 Calcular el peso asignado al weak learner \(G_m\): \[\alpha_m = log(\frac{1-err_m}{err_m})\]
- 2.4 Actualizar los pesos de las observaciones: \[w_i = w_i exp[\alpha_m I(y_i \neq G_m(x_i))], \ \ \ i = 1, 2,..., N\]
Predicción de nuevas observaciones agregando todos los weak learners y ponderándolos por su peso:

\[G(x) = sign[\sum^M_{m=1} \alpha_mG_m(x)]\]

Gradient Boosting

Gradient Boosting es una generalización del algoritmo AdaBoost que permite emplear cualquier función de coste, siempre que esta sea diferenciable. La flexibilidad de este algoritmo ha hecho posible aplicar boosting a multitud de problemas (regresión, clasificación múltiple…) convirtiéndolo en uno de los métodos de machine learning de mayor éxito. Si bien existen varias adaptaciones, la idea general de todas ellas es la misma: entrenar modelos de forma secuencial, de forma que cada modelo ajusta los residuos (errores) de los modelos anteriores.

Se ajusta un primer weak learner \(f_1\) con el que se predice la variable respuesta \(y\), y se calculan los residuos \(y-f_1(x)\). A continuación, se ajusta un nuevo modelo \(f_2\), que intenta predecir los residuos del modelo anterior, en otras palabras, trata de corregir los errores que ha hecho el modelo \(f_1\).

\[f_1(x) ≈ y\]

\[f_2(x) ≈ y - f_1(x)\]

En la siguiente iteración, se calculan los residuos de los dos modelos de forma conjunta \(y-f_1(x)-f_2(x)\), los errores cometidos por \(f_1\) y que \(f_2\) no ha sido capaz de corregir, y se ajusta un tercer modelo \(f_3\) para tratar de corregirlos.

\[f_3(x) ≈ y - f_1(x) - f_2(x)\]

Este proceso se repite \(M\) veces, de forma que cada nuevo modelo minimiza los residuos (errores) del anterior.

Dado que el objetivo de Gradient Boosting es ir minimizando los residuos iteración a iteración, es susceptible de overfitting. Una forma de evitar este problema es empleando un valor de regularización, también conocido como learning rate (\(\lambda\)), que limite la influencia de cada modelo en el conjunto del ensemble. Como consecuencia de esta regularización, se necesitan más modelos para formar el ensemble pero se consiguen mejores resultados.

\[f_1(x) ≈ y\]

\[f_2(x) ≈ y - \lambda f_1(x)\]

\[f_3(x) ≈ y - \lambda f_1(x) - \lambda f_2(x)\]

\[y ≈ \lambda f_1(x) + \lambda f_2(x) + \lambda f_3(x) + \ ... \ + \lambda f_m(x)\]

Stochastic Gradient Boosting

Tiempo después de la publicación de algoritmo de Gradient Boosting, se le incorporó una de las propiedades de bagging, el muestreo aleatorio de observaciones de entrenamiento. En concreto, en cada iteración del algoritmo, el weak learner se ajusta empleando únicamente una fracción del set de entrenamiento, extraída de forma aleatoria y sin reemplazo (no con bootstrapping). Al resultado de esta modificación se le conoce como Stochastic Gradient Boosting y aporta dos ventajas: consigue mejorar la capacidad predictiva y permite estimar el out-of-bag-error.

A diferencia de Random Forest el out-of-bag-error de Stochastic Gradient Boosting solo es una buena aproximación del error de validación cruzada cuando el número de árboles en bajo. Ver detalles en Gradient Boosting Out-of-Bag estimates.

Binning

El principal cuello de botella que ralentiza al algoritmo de Gradient Boosting es la búsqueda de los puntos de corte (thresholds). Para encontrar el threshold óptimo en cada división, en cada árbol, es necesario iterar sobre todos los valores observados de cada uno de los predictores. El problema no es que haya muchos valores que comparar, sino que requiere ordenar cada vez las observaciones acorde a cada predictor, proceso computacionalmente muy costoso.

La estrategia de binning trata de reducir este cuello de botella agilizando el ordenamiento de las observaciones mediante una discretización de sus valores. Normalmente se emplean los cuantiles para hacer una discretización homogénea en cuanto al número de observaciones que cae en cada intervalo. Como resultado de este proceso, el ordenamiento es varios órdenes de magnitud más rápido. El potencial inconveniente de la discretización es que se pierde parte de la información, ya que únicamente se contemplan como posibles thresholds los límites de cada bin.

Esta es la principal característica que hace que las implementaciones de XGBoost, LightGBM y H2O sean mucho más rápida que la implementación original de gbm.

Parada temprana (early stopping)

Una de las características de los modelos Gradient Boosting es que, con el número suficiente de weak learners, el modelo final tiende a ajustarse perfectamente a los datos de entrenamiento causando overfitting. Este comportamiento implica que el analista tiene que encontrar el número adecuado de árboles y, para ello, suele tener que entrenar el modelo con cientos o miles de árboles hasta identificar el momento en el que empieza el overfitting. Esto suele ser poco eficiente en términos de tiempo, ya que, posiblemente, se estén ajustando muchos árboles innecesarios.

Para evitar este problema, la mayoría de implementaciones incluyen toda una serie de estrategias para detener el proceso de ajuste del modelo a partir del momento en el que este deja de mejorar. Por defecto, la métrica empleada se calcula utilizando un conjunto de validación. En XGBoost, la estrategia de parada está controlada por el argumento early_stopping_rounds.

Importancia de los predictores

Las estrategias de gradient boosting para determinar la importancia de los predictores son las mismas que las explicadas para random forest: pureza de nodos y permutación.

Ejemplo regresión

Librerías

Las librerías utilizadas en este ejemplo son:

# Tratamiento de datos
# ==============================================================================
library(MASS)
library(dplyr)
library(tidyr)
library(skimr)

# Gráficos
# ==============================================================================
library(ggplot2)
library(ggpubr)

# Preprocesado y modelado
# ==============================================================================
library(tidymodels)
library(xgboost)
library(gbm)
library(doParallel)

Datos

data("Boston")
datos = Boston
head(datos, 3)

Ajuste del modelo

Se ajusta un modelo empleando como variable respuesta medv y como predictores todas las otras variables disponibles.

La función xgboost() del paquete xgboost permite entrenar modelos Gradient Boosting para problemas de regresión y clasificación. Los parámetros e hiperparámetros empleados:

params = list()
data
nrounds
watchlist = list()
obj = NULL
feval = NULL
verbose = 1
print_every_n = 1L
early_stopping_rounds = NULL
maximize = NULL
save_period = NULL
save_name = "xgboost.model"
xgb_model = NULL
callbacks = list()

Son muchos los parámetros que regulan el comportamiento de XGBoost, puede encontrarse una descripción detallada de todos ellos en XGBoost Parameters. En la práctica, cabe prestar especial atención a aquellos que controlan el crecimiento de los árboles, la velocidad de aprendizaje del modelo, y los que gestionan la parada temprana para evitar overfitting:

params: lista con los parámetros de entrenamiento:
- booster [default = gbtree ]: tipo de modelo empleado como weak learner, pueden ser árboles (“gbtree”, “dart”) o modelos lineales (“gblinear”)
- eta [default=0.3, alias: learning_rate]: reduce la contribución de cada árbol multiplicando su influencia original por este valor.
- gamma [default=0, alias: min_split_loss]: reducción mínima en la función de coste que debe de conseguir una división para que se lleve a cabo.
- max_depth [default=6]: profundidad máxima que pueden alcanzar los árboles.
- subsample [default=1]: proporción de observaciones utilizadas para el ajsute de cada árbol. Si su valor es inferior a 1, se está aplicando Stochastic Gradient Boosting.
- colsample_bytree: número de predictores considerados a en cada división.
nrounds: número de iteraciones del algoritmo de boosting, es decir, número de modelos que forman el ensemble.
early_stopping_rounds: número de iteraciones consecutivas sin mejora para que el algoritmo se detenga (early stopping). Si su valor es NULL, se desactiva la parada temprana. Para que la parada temprana se pueda llevar a cabo, se necesita un conjunto de validación independiente del de entrenamiento (watchlist).
watchlist: conjunto de validación empleado la la parada temprana.
seed: semilla para que los resultados sean reproducibles. Este parámetro se ignora, utilizar set.seed() en su lugar.

# División de los datos en train y test
# ==============================================================================
set.seed(123)
train <- sample(1:nrow(datos), size = nrow(datos)/2)
datos_train <- datos[train,]
datos_test  <- datos[-train,]

Los modelos de XGBoost pueden trabajar con varios formatos de datos, entre ellos las matrices de R. Sin embargo, es recomendable utilizar xgb.DMatrix, una estructura propia y optimizada esta librería.

datos_train <- xgb.DMatrix(
                data  = datos_train %>% select(-medv) %>% data.matrix(),
                label = datos_train$medv
               )

datos_test <- xgb.DMatrix(
                data  = datos_test %>% select(-medv) %>% data.matrix(),
                label = datos_test$medv
               )

# Creación y entrenamiento del modelo
# ==============================================================================
set.seed(123)
modelo <- xgb.train(
            data    = datos_train,
            params  = list(max_depth = 2),
            nrounds = 10,
            
          )
modelo

## ##### xgb.Booster
## raw: 6.7 Kb 
## call:
##   xgb.train(params = list(max_depth = 2), data = datos_train, nrounds = 10)
## params (as set within xgb.train):
##   max_depth = "2", validate_parameters = "1"
## xgb.attributes:
##   niter
## callbacks:
##   cb.print.evaluation(period = print_every_n)
## # of features: 13 
## niter: 10
## nfeatures : 13

Predicción y evaluación del modelo

Una vez entrenado el modelo, se evalúa la capacidad predictiva empleando el conjunto de test.

# Error de test del modelo
# ==============================================================================
predicciones <- predict(
                  modelo,
                  newdata = datos_test
                )

test_rmse <- sqrt(mean((predicciones - getinfo(datos_test, "label"))^2))
paste("Error de test (rmse) del modelo:", round(test_rmse,2))

## [1] "Error de test (rmse) del modelo: 4.26"

Optimización de hiperparámetros

El modelo inicial se ha entrenado utilizando 10 árboles (nrounds = 10) y manteniendo el resto de hiperparámetros con su valor por defecto. Al ser hiperparámetros, no se puede saber de antemano cuál es el valor más adecuado, la forma de identificarlos es mediante el uso de estrategias de validación, por ejemplo validación cruzada.

Número de árboles

En Gradient Boosting, el número de árboles es un hiperparámetro crítico en cuanto que, con forme se añaden árboles, se incrementa el riesgo de overfitting.

# Validación empleando k-cross-validation (root mean squared error)
# ==============================================================================

resultados_cv <- xgb.cv(
                data      = datos_train,
                params    = list(eta = 0.3, max_depth = 6, subsample = 1),
                nrounds   = 500,
                nfold     = 5,
                metrics   = "rmse", 
                verbose   = 0
             )

resultados_cv <- resultados_cv$evaluation_log

# Gráfico con la evolución de los errores
ggplot(data = resultados_cv) +
  geom_line(aes(x = iter, y = train_rmse_mean, color = "train rmse")) + 
  geom_line(aes(x = iter, y = test_rmse_mean, color = "cv rmse")) +
  geom_point(
    data = slice_min(resultados_cv, order_by = test_rmse_mean, n = 1),
    aes(x = iter, y = test_rmse_mean),
    color = "firebrick"
  ) +
  labs(
    title = "Evolución del cv-error vs número árboles",
    x     = "número de árboles",
    y     = "cv-error (rmse)",
    color = ""
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

paste("Valor óptimo de nrounds:", slice_min(resultados_cv, order_by = test_rmse_mean, n = 1)$iter)

## [1] "Valor óptimo de nrounds: 34"

Learning rate

Junto con el número de árboles, el learning rate (eta) es el hiperparámetro más importantes en Gradient Boosting, ya que es el que permite controlar cómo de rápido aprende el modelo y con ello el riesgo de llegar al overfitting. Estos dos hiperparámetros son interdependientes, cuanto menor es el learning rate, más árboles se necesitan para alcanzar buenos resultados pero menor es el riesgo de overfitting.

eta_range          <- c(0.001, 0.01, 0.1, 0.3)
df_resultados_cv   <- data.frame()

for(i in seq_along(eta_range)){
  set.seed(123)
  resultados_cv <- xgb.cv(
                    data    = datos_train,
                    params  = list(eta = eta_range[i], max_depth = 6, subsample = 1),
                    nrounds = 4000,
                    nfold   = 5,
                    metrics = "rmse", 
                    verbose = 0
                 )
  
  resultados_cv <- resultados_cv$evaluation_log
  resultados_cv <- resultados_cv %>%
                   select(iter, test_rmse_mean) %>%
                   mutate(eta = as.character(eta_range[i]))
  
  df_resultados_cv <- df_resultados_cv %>% bind_rows(resultados_cv)
  
}


# Gráfico con la evolución de los errores
ggplot(data = df_resultados_cv) +
  geom_line(aes(x = iter, y = test_rmse_mean, color = eta)) + 
  labs(
    title = "Evolución del cv-error vs learning rate (eta)",
    x     = "eta",
    y     = "cv-error (rmse)",
    color = ""
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

Los valores estimados por validación cruzada indican que, el mejor modelo se consigue con un eta de 0.3.

Grid search

Grid Search basado en validación cruzada

# División de los datos en train y test
# ==============================================================================
set.seed(123)
train <- sample(1:nrow(datos), size = nrow(datos)/2)
datos_train <- datos[train,]
datos_test  <- datos[-train,]

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# ==============================================================================
set.seed(123)
modelo <- boost_tree(
             mode        = "regression",
             mtry        = tune(),
             trees       = tune(),
             tree_depth  = tune(),
             learn_rate  = tune(),
             sample_size = tune(),
             stop_iter   = NULL
          ) %>%
          set_engine(
            engine     = "xgboost"
          )

# DEFINICIÓN DEL PREPROCESADO
# ==============================================================================
# En este caso no hay preprocesado, por lo que el transformer solo contine
# la definición de la formula y los datos de entrenamiento.
transformer <- recipe(
                  formula = medv ~ .,
                  data    =  datos_train
               )

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# ==============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = medv
            )

# WORKFLOW
# ==============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo)
                     

# GRID DE HIPERPARÁMETROS
# ==============================================================================
hiperpar_grid <- expand_grid(
                  "trees"      = c(100, 500, 1000),
                  "mtry"       = c(ceiling(sqrt(ncol(datos_train))), ncol(datos_train)-1),
                  "tree_depth" = c(1, 3, 10, 20),
                  "learn_rate" = c(0.01, 0.1, 0.3),
                  "sample_size" = c(0.5, 1)
                 )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# ==============================================================================
cl <- makePSOCKcluster(parallel::detectCores() - 1)
registerDoParallel(cl)

set.seed(123)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              grid      = hiperpar_grid
            )

stopCluster(cl)

# Mejores hiperparámetros por validación cruzada
# ==============================================================================
show_best(grid_fit, metric = "rmse", n = 1)

Una vez identificados los mejores hiperparámetros, se reentrena el modelo indicando los valores óptimos en sus argumentos.

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_final_fit <- finalize_workflow(
                        x = workflow_modelado,
                        parameters = mejores_hiperpar
                    ) %>%
                    fit(
                      data = datos_train
                    ) %>%
                    pull_workflow_fit()

# Error de test del modelo final
# ==============================================================================
predicciones <- modelo_final_fit %>%
                predict(
                  new_data = bake(transformer, datos_test %>% select(-medv)),
                  type     = "numeric"
                )

predicciones <- predicciones %>% 
                bind_cols(datos_test %>% dplyr::select(medv))

rmse_test  <- rmse(
                 data     = predicciones,
                 truth    = medv,
                 estimate = .pred,
                 na_rm    = TRUE
              )
rmse_test

Tras optimizar los hiperparámetros, se consigue reducir el error rmse del modelo de 4.26 a 3.63. Las predicciones del modelo final se alejan en promedio 3.63 unidades (3630 dólares) del valor real.

Importancia de predictores

La función xgb.importance() calcula la importancia de los predictores de un modelo XGBoost en base a la pureza de nodos.

importance_matrix <- xgb.importance(model = modelo_final_fit$fit)
xgb.plot.importance(importance_matrix = importance_matrix)

Gradient Boosting Monótono

Los modelos basados en Gradient Boosting son capaces de aprender interacciones complejas entre los predictores y la variable respuesta. Esto puede convertirse en una desventaja cuando, debido al ruido de los datos, el modelo aprende comportamientos que no reflejan la realidad. Una estrategia para evitar este tipo de problemas consiste en forzar a que la relación de cada predictor (numérico) con la variable respuesta sea monótona, es decir tenga siempre la misma dirección (positiva o negativa).

Aunque el algoritmo es más complejo, a grandes rasgos, se consigue una relación monótona respecto a una determinada variable si, además de otros criterios, solo se permiten divisiones de nodos cuando el valor medio del nodo hijo derecho es mayor que el nodo hijo izquierdo (o a la inversa si se quiere una relación negativa).

Los siguientes datos muestran una relación no lineal positiva entre X e Y.

x <- runif(n = 200, min = 0, max = 15)
y <- x^2 + 10
datos <- data.frame(x, y)

ggplot(data = datos, aes(x = x, y = y)) +
  geom_point() +
  labs(title = "Relación Y~X verdadera") + 
  theme_bw()

Supóngase ahora que los datos contienen cierto ruido y que, por lo tanto, la relación observada no es exactamente la verdadera.

datos$y <- datos$y + rnorm(n = 200, mean = 0, sd = 20)
ggplot(data = datos, aes(x = x, y = y)) +
  geom_point() +
  labs(title = "Relación Y~X con ruido") + 
  theme_bw()

Si se ajusta un modelo de gradient boosting a estos datos el modelo asume como verdadero parte del ruido, lo que conlleva que, en determinadas zonas, considere negativa la relación entre X e Y.

library(gbm)
set.seed(123)
modelo_gbm <- gbm(formula = y ~ x,
                  data = datos,
                  n.tree = 500,
                  interaction.depth = 5,
                  shrinkage = 0.1,
                  distribution = "gaussian",
                  verbose = FALSE)

# Se crea un grid de predicciones para representar el modelo.
grid <- data.frame(x = runif(n = 1000,
                             min = min(datos$x),
                             max = max(datos$x)))
grid$prediccion <- predict(object = modelo_gbm,
                           newdata = grid,
                           n.trees = 500)

# Representación del modelo.
ggplot(data = datos, aes(x = x, y = y)) +
  geom_point() +
  geom_line(data = grid,
            aes(x = x, y = prediccion), 
            color = "red", size = 1) +
  labs(title = "Modelo GBM") +
  theme_bw()

Si se dispone de suficiente conocimiento sobre el problema en cuestión como para saber que la relación entre predictor y variable es monótona, forzar el comportamiento puede resultar en un modelo con mayor capacidad de generalizar a nuevas observaciones. La función gbm() tiene un argumento llamado var.monotone que recibe un vector de longitud igual al número de predictores del modelo, con valores +1 para monotonía positiva, -1 para negativa y 0 para no restricciones. Esta misma característica tambien esta dispobible en H2O.

Se repite el ajuste del modelo, esta vez forzando la monotonía positiva.

set.seed(123)
modelo_gbm <- gbm(formula = y ~ x,
                  data = datos,
                  n.tree = 500,
                  interaction.depth = 5,
                  shrinkage = 0.1,
                  distribution = "gaussian",
                  var.monotone = c(1),
                  verbose = FALSE)

# Se crea un grid de predicciones para representar el modelo.
grid <- data.frame(x = runif(n = 1000,
                             min = min(datos$x),
                             max = max(datos$x)))
grid$prediccion <- predict(object = modelo_gbm,
                           newdata = grid,
                           n.trees = 500)

# Representación del modelo.
ggplot(data = datos, aes(x = x, y = y)) +
  geom_point() +
  geom_line(data = grid,
            aes(x = x, y = prediccion), 
            color = "red", size = 1) +
  labs(title = "Modelo GBM",
       subtitle = "Monotonía positiva") +
  theme_bw()

Esta vez, aunque la relación y~x no es constante, si es siempre positiva.

C5.0

Idea intuitiva

C5.0 es el algoritmo sucesor de C4.5, publicado por John Ross Quinlan (1992), con el objetivo de crear árboles de clasificación. Entre sus características, destacan la capacidad para generar árboles de decisión simples, modelos basados en reglas, ensembles basados en boosting y asignación de distintos pesos a los errores. Este algoritmo ha resultado de una gran utilidad a la hora de crear modelos de clasificación y todas sus capacidades son accesibles mediante el paquete C50. Aunque comparte muchas características con los algoritmos de random forest y gradient boosting, cabe tener en cuenta algunas peculiaridades:

La medida de pureza empleada para las divisiones del árbol es la entropía.
El podado de los árboles se realiza por defecto, y el método empleado se conoce como pessimistic pruning.
Los árboles se pueden convertir en modelos basados en reglas.
Emplea un algoritmo de boosting más próximo a AdaBoost que a Gradient Boosting.
Por defecto, el algoritmo de boosting se detiene si la incorporación de nuevos modelos no aporta un mínimo de mejora.
Incorpora una estrategia para la selección de predictores (Winnowing) previo ajuste del modelo.
Permite asignar diferente peso a cada tipo de error.

Ejemplo

Supóngase una entidad bancaria que se dedica a conceder préstamos a sus clientes. La compañía ha detectado que parte de sus pérdidas económicas se deben a que, un porcentaje de clientes a los que se les concede un préstamo, nunca lo devuelven. El banco está interesado en un sistema de aprendizaje estadístico que ayude a identificar a qué clientes denegar el préstamo. Ejemplo obtenido del libro Machine Learning with R Second Edition by Brett Lantz.

Se trata de un problema de clasificación binaria (si/no) que puede resolverse con un número considerable de herramientas (regresión logística, maquinas de vector soporte, redes neuronales…). En este escenario, además de la predicción, es importante poder entender las razones por las que se concede o no un préstamo. En situaciones de este tipo, en las que la transparencia del modelo añade mucho valor, los sistemas basados en árboles suelen ser la elección adecuada. A continuación, se exploran las capacidades del algoritmo C5.0 para resolver problemas de clasificación.

Librerías

Las librerías utilizadas en este ejemplo son:

# Tratamiento de datos
# ==============================================================================
library(tidyverse)
library(skimr)

# Gráficos
# ==============================================================================
library(ggplot2)
library(ggpubr)

# Preprocesado y modelado
# ==============================================================================
library(C50)
library(tidymodels)
library(doParallel)

Datos

El set de datos credit.csv (puede descargarse registrándose de forma gratuita en packtpub), contiene información sobre 1000 clientes. Para cada uno de ellos, se han registrado 17 variables de tipo cualitativo y cuantitativo. La variable default indica si el cliente defraudó o no.

datos <- read_csv(file = "./datos/credit.csv")
head(datos,3)

skim(datos)

Data summary
Name	datos
Number of rows	1000
Number of columns	17
_______________________
Column type frequency:
character	10
numeric	7
________________________
Group variables	None

Variable type: character

skim_variable	complete_rate	min	max	n_unique
checking_balance	1	6	10	4
credit_history	1	4	9	5
purpose	1	3	20	6
savings_balance	1	7	13	5
employment_duration	1	8	11	5
other_credit	1	4	5	3
housing	1	3	5	3
job	1	7	10	4
phone	1	2	3	2
default	1	2	3	2

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
months_loan_duration	1	20.90	12.06	4	12.0	18.0	24.00	72	▇▇▂▁▁
amount	1	3271.26	2822.74	250	1365.5	2319.5	3972.25	18424	▇▂▁▁▁
percent_of_income	1	2.97	1.12	1	2.0	3.0	4.00	4	▂▃▁▂▇
years_at_residence	1	2.85	1.10	1	2.0	3.0	4.00	4	▂▆▁▃▇
age	1	35.55	11.38	19	27.0	33.0	42.00	75	▇▆▃▁▁
existing_loans_count	1	1.41	0.58	1	1.0	1.0	2.00	4	▇▅▁▁▁
dependents	1	1.16	0.36	1	1.0	1.0	1.00	2	▇▁▁▁▂

De los 1000 clientes, 700 devolvieron el préstamo y 300 no. Para que un sistema de clasificación se considere útil, debe de tener una capacidad de acierto mayor que la frecuencia de la clase mayoritaria, en este caso 700/1000 = 0.7.

table(datos$default)

## 
##  no yes 
## 700 300

Separación entrenamiento y test

Para poder determinar la capacidad predictiva del modelo, se seleccionan aleatoriamente un 80% de las observaciones como entrenamiento y el 20% restante como test.

# División de los datos en train y test
# ==============================================================================
set.seed(1234)
indices     <- sample(x = 1:nrow(datos), size = 0.8*nrow(datos))
datos_train <- datos[indices,]
datos_test  <- datos[-indices,]

Cuando se generan particiones de entrenamiento y test, es importante comprobar que las proporciones de los niveles de la variable respuesta son aproximadamente iguales entre ambas particiones.

round(table(datos_train$default)/nrow(datos_train), 3)

## 
##    no   yes 
## 0.699 0.301

round(table(datos_test$default)/nrow(datos_test), 3)

## 
##    no   yes 
## 0.705 0.295

Árbol de clasificación simple

Se procede a crear un modelo de clasificación basado en un único árbol. C5.0 necesita que las variables cualitativas estén almacenadas en forma de factor.

# Preprocesado
# ==============================================================================
# Convertir todas las columnas no numéricas en factor.
datos_train <- datos_train %>% map_if(.p = is_character, .f = as.factor) %>%
               as.data.frame()
datos_test <- datos_test %>% map_if(.p = is_character, .f = as.factor) %>%
              as.data.frame()

# Ajuste del árbol
#===============================================================================
set.seed(123)
arbol_c50 <- C5.0(
              formula = default ~ ., 
              data    = datos_train,
              trials  = 1,
              rules   = FALSE,
              control = C5.0Control(seed = 123)
              
            )

El summary de un modelo C50 muestra información detallada sobre el número de observaciones de entrenamiento, el número de predictores empleados, las divisiones del árbol, los errores de entrenamiento y la importancia de los predictores.

summary(arbol_c50)

## 
## Call:
## C5.0.formula(formula = default ~ ., data = datos_train, trials = 1, rules
##  = FALSE, control = C5.0Control(seed = 123))
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Fri Nov  6 10:58:49 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 800 cases (17 attributes) from undefined.data
## 
## Decision tree:
## 
## checking_balance in {> 200 DM,unknown}:
## :...other_credit in {none,store}: no (326/34)
## :   other_credit = bank:
## :   :...purpose in {car0,furniture/appliances,renovations}: no (20/1)
## :       purpose in {business,car,education}:
## :       :...percent_of_income <= 1: no (4)
## :           percent_of_income > 1:
## :           :...years_at_residence > 3: no (10/3)
## :               years_at_residence <= 3:
## :               :...months_loan_duration <= 12: no (2)
## :                   months_loan_duration > 12: yes (10/1)
## checking_balance in {< 0 DM,1 - 200 DM}:
## :...credit_history in {perfect,very good}:
##     :...housing in {other,rent}: yes (26/2)
##     :   housing = own:
##     :   :...purpose in {car0,renovations}: yes (3)
##     :       purpose = education: no (1)
##     :       purpose = business:
##     :       :...years_at_residence <= 1: no (2)
##     :       :   years_at_residence > 1: yes (3/1)
##     :       purpose = car:
##     :       :...months_loan_duration <= 18: yes (5)
##     :       :   months_loan_duration > 18:
##     :       :   :...job in {management,skilled,unemployed}: no (3)
##     :       :       job = unskilled: yes (1)
##     :       purpose = furniture/appliances:
##     :       :...checking_balance = 1 - 200 DM: no (4)
##     :           checking_balance = < 0 DM:
##     :           :...other_credit = bank: yes (3)
##     :               other_credit in {none,store}: no (2)
##     credit_history in {critical,good,poor}:
##     :...months_loan_duration > 27:
##         :...dependents > 1:
##         :   :...other_credit in {bank,store}: yes (3)
##         :   :   other_credit = none:
##         :   :   :...age <= 45: no (10/1)
##         :   :       age > 45: yes (2)
##         :   dependents <= 1:
##         :   :...savings_balance in {< 100 DM,500 - 1000 DM}: yes (52/11)
##         :       savings_balance = > 1000 DM: no (2/1)
##         :       savings_balance = 100 - 500 DM:
##         :       :...credit_history = critical: no (1)
##         :       :   credit_history = good: yes (7)
##         :       :   credit_history = poor:
##         :       :   :...job = management: yes (2)
##         :       :       job in {skilled,unemployed,unskilled}: no (3)
##         :       savings_balance = unknown:
##         :       :...checking_balance = 1 - 200 DM: no (7/1)
##         :           checking_balance = < 0 DM:
##         :           :...credit_history = critical: no (1)
##         :               credit_history in {good,poor}: yes (4)
##         months_loan_duration <= 27:
##         :...credit_history = poor:
##             :...checking_balance = < 0 DM:
##             :   :...months_loan_duration <= 18: no (2)
##             :   :   months_loan_duration > 18: yes (4)
##             :   checking_balance = 1 - 200 DM:
##             :   :...age <= 45: no (13)
##             :       age > 45: yes (3/1)
##             credit_history = critical:
##             :...other_credit = bank: yes (10/4)
##             :   other_credit = store: no (1)
##             :   other_credit = none:
##             :   :...employment_duration in {< 1 year,4 - 7 years}: no (25)
##             :       employment_duration in {> 7 years,1 - 4 years,unemployed}:
##             :       :...purpose in {business,car0,education,
##             :           :           renovations}: yes (4)
##             :           purpose in {car,furniture/appliances}: no (36/8)
##             credit_history = good:
##             :...purpose in {car0,renovations}: no (8/2)
##                 purpose = business:
##                 :...age <= 37: no (11)
##                 :   age > 37: yes (3/1)
##                 purpose = education:
##                 :...amount <= 1050: yes (3)
##                 :   amount > 1050: no (3)
##                 purpose = furniture/appliances:
##                 :...savings_balance in {< 100 DM,> 1000 DM,
##                 :   :                   500 - 1000 DM}: no (79/22)
##                 :   savings_balance = 100 - 500 DM: yes (7/1)
##                 :   savings_balance = unknown:
##                 :   :...years_at_residence <= 2: yes (4)
##                 :       years_at_residence > 2: no (6/1)
##                 purpose = car:
##                 :...amount > 8471: yes (6)
##                     amount <= 8471:
##                     :...employment_duration = unemployed: no (8)
##                         employment_duration in {< 1 year,> 7 years,1 - 4 years,
##                         :                       4 - 7 years}:
##                         :...amount <= 1577: yes (22/5)
##                             amount > 1577:
##                             :...savings_balance = 500 - 1000 DM: yes (1)
##                                 savings_balance in {> 1000 DM,
##                                 :                   unknown}: no (5)
##                                 savings_balance = 100 - 500 DM:
##                                 :...housing in {other,own}: no (3)
##                                 :   housing = rent: yes (1)
##                                 savings_balance = < 100 DM:
##                                 :...checking_balance = 1 - 200 DM: no (5)
##                                     checking_balance = < 0 DM:
##                                     :...housing = rent: yes (0)
##                                         housing = other: no (1)
##                                         housing = own:
##                                         :...amount <= 2096: no (2)
##                                             amount > 2096: yes (5)
## 
## 
## Evaluation on training data (800 cases):
## 
##      Decision Tree   
##    ----------------  
##    Size      Errors  
## 
##      58  101(12.6%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##     532    27    (a): class no
##      74   167    (b): class yes
## 
## 
##  Attribute usage:
## 
##  100.00% checking_balance
##   58.50% other_credit
##   53.50% credit_history
##   49.50% months_loan_duration
##   37.00% purpose
##   24.75% savings_balance
##   14.75% employment_duration
##   11.75% dependents
##    8.12% amount
##    8.12% housing
##    5.25% age
##    4.62% years_at_residence
##    3.25% percent_of_income
##    1.12% job
## 
## 
## Time: 0.0 secs

En este caso, 101 de los 800 clientes han sido clasificados incorrectamente, el error de entrenamiento es del 12.6%.

Error de test

Una vez ajustado el modelo, se evalúa su capacidad predictiva con el conjunto de test.

# ==============================================================================
predicciones <- predict(arbol_c50, newdata = datos_test, type = "class")
table(predicho = predicciones, real = datos_test$default)

##         real
## predicho  no yes
##      no  113  33
##      yes  28  26

error_clas <- mean(predicciones != datos_test$default)
paste(
  "El error de clasificación es del:", 100 * error_clas, "%.", 
  sum(predicciones == datos_test$default),
  "clasificaciones correctas de un total de",
  length(predicciones)
)

## [1] "El error de clasificación es del: 30.5 %. 139 clasificaciones correctas de un total de 200"

Si se especifica el argumento type = "prob", en lugar de la clase predicha se devuelven las probabilidades de pertenecer a cada clase.

predicciones <- predict(arbol_c50, newdata = datos_test, type = "prob")
head(predicciones, 3)

##          no       yes
## 1 0.8951032 0.1048968
## 2 0.4246875 0.5753125
## 3 0.8951032 0.1048968

El modelo basado en un árbol simple no consigue clasificar correctamente las observaciones mejor de lo que cabría esperar si se asignara a todas a la clase mayoritaria.

Boosting

C5.0 incorpora un método propio de boosting (semejante a AdaBoost) para generar modelos de tipo ensemble basados en árboles.

# Ajuste del modelo
#===============================================================================
modelo_c50_boost <- C5.0(
                      formula = default ~ ., 
                      data    = datos_train,
                      trials  = 100,
                      rules   = FALSE,
                      control = C5.0Control(seed = 123)
                    )

# Error de test
#===============================================================================
predicciones <- predict(modelo_c50_boost, newdata = datos_test, type = "class")
table(predicho = predicciones, real = datos_test$default)

##         real
## predicho  no yes
##      no  128  36
##      yes  13  23

error_clas <- mean(predicciones != datos_test$default)
paste("El error de clasificación es del:", 100 * error_clas, "%.", 
      sum(predicciones == datos_test$default),
      "clasificaciones correctas de un total de", length(predicciones))

## [1] "El error de clasificación es del: 24.5 %. 151 clasificaciones correctas de un total de 200"

Los resultados muestran que, mediante boosting, se mejora la capacidad predictiva del modelo.

Identificación de variables más influyentes

C5.0 incorpora dos formas de cuantificar la influencia de los predictores en el modelo final (simple o ensemble).

Usage: porcentaje de observaciones de entrenamiento que caen en todos los nodos generados tras una división en el que ha participado el predictor. Por ejemplo, en el caso de un árbol simple, el predictor empleado en la primera decisión recibe automáticamente una importancia del 100%, ya que todas las observaciones de entrenamiento caen en uno de los dos nodos hijos generados. Otros predictores puede que participen con frecuencia en divisiones, pero si en los nodos hijos solo caen unas pocas observaciones, su importancia puede ser próxima a cero. En el caso de boosting, se promedia la importancia de cada predictor en todos los árboles que forman el ensemble.
Splits: porcentaje de divisiones en las que participa cada predictor. No tiene en cuenta la importancia de la división.

Las dos medidas cuantifican aspectos muy diferentes, por lo que el ranking de importancia puede variar mucho dependiendo cual se utilice.

importancia_usage <- C5imp(modelo_c50_boost, metric = "usage")
importancia_usage <- importancia_usage %>%
                     rownames_to_column(var = "predictor")

importancia_splits <- C5imp(modelo_c50_boost, metric = "splits")
importancia_splits <- importancia_splits %>%
                     rownames_to_column(var = "predictor")


p1 <- ggplot(data = importancia_usage, aes(x = reorder(predictor, Overall),
                                           y = Overall, fill = Overall)) +
    labs(x = "predictor", title = "usage") +
    geom_col() +
    coord_flip() +
    scale_fill_viridis_c() +
    theme_bw() +
    theme(legend.position = "none")

p2 <- ggplot(data = importancia_splits, aes(x = reorder(predictor, Overall),
                                            y = Overall, fill = Overall)) +
    labs(x = "predictor", title = "splits") +
    geom_col() +
    coord_flip() +
    scale_fill_viridis_c() +
    theme_bw() +
    theme(legend.position = "none")
ggarrange(p1, p2)

Definir el peso de errores

En muchos escenarios, las consecuencias de cometer un error de clasificación son distintas dependiendo el tipo que sean. Por ejemplo, en el ámbito biomédico, no es igual de grave confundir un tumor maligno con uno benigno que viceversa, ya que, en el primer caso, la vida del paciente puede estar en peligro. En el problema de los préstamos de este ejercicio, concederle un préstamo a alguien que no lo va a devolver, es mucho más costoso para el banco que no concedérselo a varios clientes que sí lo devolverían. El algoritmo C5.0 permite asignar diferentes penalizaciones a cada tipo de error, lo que fuerza al modelo (árbol) a minimizar determinados tipos de error.

El peso o penalización de cada error se pasa en forma de matriz al argumento costs. En la matriz se especifica cuánto de costoso es cada error respecto a los demás. En este problema, dado que la variable respuesta es binaria, se necesita una matriz 2x2.

dimensiones <- list(predicho = c("no", "yes"), real = c("no", "yes"))
dimensiones

## $predicho
## [1] "no"  "yes"
## 
## $real
## [1] "no"  "yes"

A continuación, se asigna la penalización de cada tipo de error. Se considera que un préstamo no recuperado es 4 veces más costoso para el banco que un cliente perdido. Las predicciones correctas no tienen coste.

coste_errores <- matrix(data = c(0, 1, 4, 0), nrow = 2, dimnames = dimensiones)
coste_errores

##         real
## predicho no yes
##      no   0   4
##      yes  1   0

Por último, se reajusta un modelo, esta vez empleando los pesos y se comparan los resultados frente al modelo que emplea el mismo peso para todos los errores.

# Mismo peso para todos los errores
modelo_c50_boost <- C5.0(
                      formula = default ~ ., 
                      data    = datos_train,
                      rules   = FALSE,
                      control = C5.0Control(seed = 123)
                    )
predicciones <- predict(modelo_c50_boost, newdata = datos_test, type = "class")
table(predicho = predicciones, real = datos_test$default)

##         real
## predicho  no yes
##      no  113  33
##      yes  28  26

error_clas <- mean(predicciones != datos_test$default)
paste("El error de clasificación es del:", 100 * error_clas,"%")

## [1] "El error de clasificación es del: 30.5 %"

# Distinto peso para cada tipo de error
modelo_c50_boost_2 <- C5.0(
                      formula = default ~ ., 
                      data    = datos_train,
                      rules   = FALSE,
                      costs   = coste_errores,
                      control = C5.0Control(seed = 123)
                    )

predicciones <- predict(modelo_c50_boost_2, newdata = datos_test, type = "class")
table(predicho = predicciones, real = datos_test$default)

##         real
## predicho no yes
##      no  88  12
##      yes 53  47

error_clas <- mean(predicciones != datos_test$default)
paste("El error de clasificación es del:", 100 * error_clas,"%")

## [1] "El error de clasificación es del: 32.5 %"

Como consecuencia de las penalizaciones, el error total del modelo aumenta, pero el tipo de error más costoso se reduce considerablemente. El número de clientes clasificados como “no fraude” pero que realmente sí defraudan se ha reducido de 33 a 12.

Winnowing

Winnowing es una estrategia incorporada en C5.0 seleccionar automáticamente los predictores más relevantes. Consiste en emplear un algoritmo inicial que identifique los predictores que están relacionados con la variable respuesta para, posteriormente, ajustar el modelo final únicamente con estos predictores. Los pasos de Winnowing son:

Se divide aleatoriamente el set de entrenamiento en dos mitades. Con una de ellas se ajusta un árbol (llamado winnowing tree) cuya finalidad es evaluar la utilidad de los predictores.
Los predictores que no se han utilizado en ninguna división del winnowing tree se identifican como no útiles.
La otra mitad de las observaciones de entrenamiento, no utilizadas en el ajuste del winnowing tree, se emplean para estimar el error del árbol. Se inicia un proceso iterativo en el que se eliminan (de uno en uno) cada uno de los predictores que forman el winnowing tree. Si, como resultado de la eliminación, el error del árbol disminuye, se considera que el predictor solo aporta ruido y se identifica como no útil.
Una vez que todos los predictores han sido identificados como útiles o no útiles, se reajusta el winnowing tree utilizando únicamente los predictores útiles. Si al hacerlo, el error del árbol (calculado con la otra mitad del set de entrenamiento) no se reduce, el proceso de winnowing se descarta y no se excluye ningún predictor. Si por lo contrario, el error sí se reduce, se lanza el ajuste C5.0 convencional, empleando todas las observaciones de entrenamiento pero solo los predictores identificados como útiles.

Apuntes varios (miscellaneous)

Comparación árboles frente a modelos lineales

La superioridad de los métodos basados en árboles frente a los métodos lineales depende del problema en cuestión. Cuando la relación entre los predictores y la variable respuesta es aproximadamente lineal, un modelo de tipo regresión lineal funciona bien y supera a los árboles de regresión.

Por contra, si la relación entre los predictores y la variable respuesta es de tipo no lineal y compleja, los métodos basados en árboles suelen superar a las aproximaciones lineales clásicas.

El procedimiento adecuado para encontrar el mejor modelo (lineal, no lineal, árboles…) consiste en estudiar el problema en cuestión, seleccionar los modelos para los que se cumplen las condiciones y comparar el test error.

Comparación Random Forest y Gradient Boosting

Ambos métodos son muy potentes y la superioridad de uno u otro depende del problema al que se apliquen. Algunos aspectos a tener en cuenta son:

Gracias al out-of-bag error, el método de Random Forest no necesita recurrir a validación cruzada para la optimización de sus hiperparámetros. Esto puede ser una ventaja muy importante cuando los requerimientos computacionales son limitantes. Esta característica también está presente en Stochastic Gradient Boosting, sin embargo, suele subestimar la mejora conseguida al añadir más iteraciones.
Random forest tiene menos hiperparámetros, lo que hace más sencillo su correcta implementación.
Si existe una proporción alta de predictores irrelevantes, Random Forest puede verse perjudicado, sobre todo a medida que se reduce el número de predictores (\(m\)) evaluados. Supóngase que, de 100 predictores, 90 de ellos no aportan información (son ruidos). Si el hiperparámetro \(m\) es igual a 3 (en cada división se evalúan 3 predictores aleatorios), es muy probable que los 3 predictores seleccionados sean de los que no aportan nada. Aun así, el algoritmo seleccionará el mejor de ellos, incorporándolo al modelo. Cuanto mayor sea el porcentaje de predictores no relevantes, mayor la frecuencia con la que esto ocurre, por lo que los árboles que forman el Random Forest contendrán predictores irrelevantes. Como consecuencia, se reduce su capacidad predictiva. En el caso de Gradient Boosting Trees, siempre se evalúan todos los predictores, por lo que los no relevantes se ignoran.
En Random Forest, cada modelo que forma el ensemble es independiente del resto, esto permite que el entrenamiento sea paralelizable.
Random Forest no sufre problemas de overfitting por muchos árboles que se agreguen.
Si se realiza una buena optimización de hiperparámetros, Gradient Boosting Trees suele obtener resultados superiores.
El tamaño final de los modelos Gradient Boosting Trees suele ser menor.
Gradient Boosting Trees suele ser más rápido prediciendo.

Creación paso a paso de un árbol CART

Master in Big Data. University of Santiago de Compostela, Manuel Mucientes

Dado el siguiente conjunto de datos con 6 observaciones, 3 variables de entrada continuas y una variable de salida cualitativa:

library(tidyverse)
x1 <- c(4,-3,3,1,-2,-3)
x2 <- c(3,-1,-2,4,3,5)
x3 <- c(-1,-1,0,0,1,5)
y  <- c(1,0,0,1,0,0)
datos <- data.frame(x1,x2,x3,y)
datos

Se construye un árbol de clasificación (sin podar) mediante CART, utilizando como criterio de división la entropía. La condición de parada debe ser que los nodos hoja sean puros (todos los ejemplos son de la misma clase). En cada nodo del árbol se debe indicar:

La variable y su valor umbral.
La entropía correspondiente.
La los nodos hoja, la clase del nodo y los ejemplos que pertenecen al mismo.

Tal como se ha visto en las descripciones teóricas, se trata de un proceso iterativo. Para cada nodo, se evalúan las condiciones de parada, si no se cumplen, se procede con los siguientes pasos:

Para cada uno de los predictores continuos (en este caso todos):
- Se ordenan de menor a mayor los valores observados del predictor.
- Se identifican aquellos puntos en los que se produce un cambio de clase en la variable respuesta.
- Se identifican como posibles puntos de cortes los valores medios entre cada par de observaciones consecutivas en las que hay cambio de clase.
- Se calcula el valor de entropía cruzada para cada uno de los puntos de corte identificados en el paso anterior.
Se selecciona el punto de corte con menor entropía y se divide el nodo acorde a este.
Se repite el proceso con cada uno de los nuevos nodos.

Predictor X1

datos %>% select(x1,y) %>% arrange(x1)

Los posibles puntos de corte para este predictor son: (-2 + 1)/2 = -0.5, (1 + 3)/2 = 2 y (4 + 3)/2 = 3.5

\[D_{(x_1,-0.5)} = \frac{3}{6}(-\frac{3}{3}\log\frac{3}{3} -\frac{0}{3}\log\frac{0}{3}) + \frac{3}{6}(-\frac{1}{3}\log\frac{1}{3} -\frac{2}{3}\log\frac{2}{3})\]

En el cálculo de la entropía, cuando una división tiene como resultado un nodo puro aparece en la ecuación la operación \(0 \cdot log(0)\), lo que equivale a \(0 \cdot - \infty\), cuyo resultado es indeterminación. Dado que se trata de un nodo puro, su entropía es nula y por lo tanto, se considera que el término estero es 0. Para que la operación no genere un NaN al implementarla en R, se calculan primero los términos de cada partición por separado, se sustituyen por cero las indeterminaciones y luego se suman.

d = c((3/6)*(-(3/3) * log(3/3) -(0/3) * log(0/3)), (3/6)*(-(1/3) * log(1/3) - (2/3) * log(2/3)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.3182571

\[D_{(x_1,2)} = \frac{4}{6}(-\frac{3}{4}\log\frac{3}{4} -\frac{1}{4}\log\frac{1}{4}) + \frac{2}{6}(-\frac{1}{2}\log\frac{1}{2} -\frac{1}{2}\log\frac{1}{2})\]

d = c((4/6)*(-(3/4) * log(3/4) -(1/4) * log(1/4)), (2/6)*(-(1/2) * log(1/2) - (1/2) * log(1/2)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.6059392

\[D_{(x_1,3.5)} = \frac{5}{6}(-\frac{4}{5}\log\frac{4}{5} -\frac{1}{5}\log\frac{1}{5}) + \frac{1}{6}(-\frac{0}{1}\log\frac{0}{1} -\frac{1}{1}\log\frac{1}{1})\]

d = c((5/6)*(-(4/5) * log(4/5) -(1/5) * log(1/5)), (12/6)*(-(0/1) * log(0/1) - (1/1) * log(1/1)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.417002

Predictor X2

datos %>% select(x2, y) %>% arrange(x2)

Cuando para un mismo valor del predictor existen dos o más observaciones con clases distintas, se establece un posible punto de corte a cada lado de la misma. En este caso, para \(X_2 = 3\) existe una observación de clase 0 y otra de clase 1. Los posibles puntos de corte para este predictor son: (-1 + 3)/2 = 1, (3 + 4)/2 = 3.5 y (4 + 5)/2 = 4.5.

\[D_{(x_2, 1)} = \frac{2}{6}(-\frac{2}{2}\log\frac{2}{2} -\frac{0}{2}\log\frac{0}{2}) + \frac{4}{6}(-\frac{2}{4}\log\frac{2}{4} -\frac{2}{4}\log\frac{2}{4})\]

d = c((2/6)*(-(2/2) * log(2/2) -(0/2) * log(0/2)), (4/6)*(-(2/4) * log(2/4) - (2/4) * log(2/4)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.4620981

\[D_{(x_2, 3.5)} = \frac{4}{6}(-\frac{3}{4}\log\frac{3}{4} -\frac{1}{4}\log\frac{1}{4}) + \frac{2}{6}(-\frac{1}{2}\log\frac{1}{2} -\frac{1}{2}\log\frac{1}{2})\]

d = c((4/6)*(-(3/4) * log(3/4) -(1/4) * log(1/4)), (2/6)*(-(1/2) * log(1/2) - (1/2) * log(1/2)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.6059392

\[D_{(x_2, 4.5)} = \frac{5}{6}(-\frac{3}{5}\log\frac{3}{5} -\frac{2}{5}\log\frac{2}{5}) + \frac{1}{6}(-\frac{1}{1}\log\frac{1}{1} -\frac{0}{1}\log\frac{0}{1})\]

d = c((5/6)*(-(3/5) * log(3/5) -(2/5) * log(2/5)), (1/6)*(-(1/1) * log(1/1) - (0/1) * log(0/1)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.5608431

Predictor X3

datos %>% select(x3, y) %>% arrange(x3)

Los posibles puntos de corte para este predictor son: (-1 + 0)/2 = -0.5 y (0 + 1)/2 = 0.5.

\[D_{(x_3, -0.5)} = \frac{2}{6}(-\frac{1}{2}\log\frac{1}{2} -\frac{1}{2}\log\frac{1}{2}) + \frac{4}{6}(-\frac{3}{4}\log\frac{3}{4} -\frac{1}{4}\log\frac{1}{4})\]

d = c((2/6)*(-(1/2) * log(1/2) -(1/2) * log(1/2)), (4/6)*(-(3/4) * log(3/4) - (1/4) * log(1/4)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.6059392

\[D_{(x_3, 0.5)} = \frac{4}{6}(-\frac{2}{4}\log\frac{2}{4} -\frac{2}{4}\log\frac{2}{4}) + \frac{2}{6}(-\frac{2}{2}\log\frac{2}{2} -\frac{0}{2}\log\frac{0}{2})\]

d = c((4/6)*(-(2/4) * log(2/4) -(2/4) * log(2/4)), (2/6)*(-(2/2) * log(2/2) - (0/2) * log(0/2)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.4620981

De todos los posibles puntos de corte, el que consigue la menor entropía total es X_1 = -0.5.

predictor <- c("x1", "x1", "x1", "x2", "x2", "x2", "x3", "x3")
valor_division <- c(-0.5, 2, 3.5, 1, 3.5, 4.5, -0.5, 0.5)
entropia  <- c(0.318, 0.606, 0.417, 0.462, 0.606, 0.56, 0.606, 0.462)
data.frame(predictor, valor_division, entropia) %>% arrange(entropia)

Esta división da lugar a dos nuevos nodos:

\(X_1 \leq -0.5\)

nodo_1 <- datos %>% filter(x1 <= -0.5)
nodo_1

\(X_1 > -0.5\)

nodo_2 <- datos %>% filter(x1 > -0.5)
nodo_2

El nodo 1 es puro, todas las observaciones que contiene pertenecen a la clase 0, por lo que se detiene su división, es un nodo terminal. El nodo 2 no cumple la condición de parada, por lo que se repite el proceso de división.

Predictor X1

nodo_2 %>% select(x1,y) %>% arrange(x1)

Los posibles puntos de corte para este predictor son: (1 + 3)/2 = 2 y (3 + 4)/2 = 3.5.

\[D_{(x_1,2)} = \frac{1}{3}(-\frac{0}{1}\log\frac{0}{1} -\frac{1}{1}\log\frac{1}{1}) + \frac{2}{3}(-\frac{1}{2}\log\frac{1}{2} -\frac{1}{2}\log\frac{1}{2})\]

d = c((1/3)*(-(0/1) * log(0/1) -(1/1) * log(1/1)), (2/3)*(-(1/2) * log(1/2) - (1/2) * log(1/2)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.4620981

\[D_{(x_1,3.5)} = \frac{2}{3}(-\frac{1}{2}\log\frac{1}{2} -\frac{1}{2}\log\frac{1}{2}) + \frac{1}{3}(-\frac{0}{1}\log\frac{0}{1} -\frac{1}{1}\log\frac{1}{1})\]

d = c((2/3)*(-(1/2) * log(1/2) - (1/2) * log(1/2)), (1/3)*(-(0/1) * log(0/1) -(1/1) * log(1/1)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.4620981

Predictor X2

nodo_2 %>% select(x2,y) %>% arrange(x2)

Los posibles puntos de corte para este predictor son: (-2 + 3)/2 = 0.5.

\[D_{(x_2,0.5)} = \frac{1}{3}(-\frac{1}{1}\log\frac{1}{1} -\frac{0}{1}\log\frac{0}{1}) + \frac{2}{3}(-\frac{0}{2}\log\frac{0}{2} -\frac{2}{2}\log\frac{2}{2})\]

d = c((1/3)*(-(1/1) * log(1/1) - (0/1) * log(0/1)), (2/3)*(-(0/2) * log(0/2) -(2/2) * log(2/2)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0

Predictor X3

nodo_2 %>% select(x3,y) %>% arrange(x3)

Los posibles puntos de corte para este predictor son: (-1 + 0)/2 = -0.5.

\[D_{(x_3, -0.5)} = \frac{1}{3}(-\frac{0}{1}\log\frac{0}{1} -\frac{1}{1}\log\frac{1}{1}) + \frac{2}{3}(-\frac{1}{2}\log\frac{1}{2} -\frac{1}{2}\log\frac{1}{2})\]

d = c((1/3)*(-(0/1) * log(0/1) - (1/1) * log(1/1)), (2/3)*(-(1/2) * log(1/2) -(1/2) * log(1/2)))
d[is.nan(d)] <- 0
sum(d)

## [1] 0.4620981

De todos los puntos de corte evaluados para el nodo 2, el que consigue la menor entropía total es X_2 = 0.5.

predictor <- c("x1", "x1", "x2", "x3")
valor_division <- c(2, 3.5, 0.5, -0.5)
entropia  <- c(0.462, 0.462, 0, 0.462)
data.frame(predictor, valor_division, entropia) %>% arrange(entropia)

Esta división da lugar a dos nuevos nodos, ambos puros, por lo que se detiene el crecimiento del árbol. Como resultado del proceso se ha generado un árbol con 3 nodos terminales.

Extremely randomized trees

Los métodos de extremely randomized trees llevan la decorrelación de los árboles un paso más allá que random forest. En cada división de los nodos, solo evalúa un subconjunto aleatorio de los predictores disponibles y, además, dentro de cada predictor seleccionado solo se evalúa un subconjunto aleatorio de los posibles puntos de corte. En determinados escenarios, este método consigue reducir un poco más la varianza. El paquete ranger ejecuta este algoritmo cuando se indica splitrule = "extratrees".

Información sesión

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)

Bibliografía

An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics)

Applied Predictive Modeling by Max Kuhn and Kjell Johnson

The Elements of Statistical Learning by T.Hastie, R.Tibshirani, J.Friedman

CMU Statistics

Trevor Hastie - Gradient Boosting Machine Learning - YouTube

Ensemble Learning to Improve Machine Learning Results

Generalized Boosted Models: A guide to the gbm package

https://en.wikipedia.org/wiki/Gradient_boosting

Classification Using C5.0

¿Cómo citar este documento?

Árboles de decisión, random forest, gradient boosting y C5.0 por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/33_arboles_decision_random_forest_gradient_boosting_C50.html

¿Te ha gustado el artículo? Tu ayuda es importante

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊

Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.

Se permite:

Compartir: copiar y redistribuir el material en cualquier medio o formato.
Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
NoComercial: No puedes utilizar el material para fines comerciales.
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.

Árboles de decisión, random forest, gradient boosting y C5.0

Joaquín Amat Rodrigo

Febrero, 2017 (última actualización Octubre, 2020)

Introducción

Paquetes R

Árboles de regresión

Idea intuitiva

Entrenamiento del árbol

Predicción del árbol

Evitar el overfitting

Controlar el tamaño del árbol (parada temprana)

Tree pruning

Ejemplo regresión

Librerías

Datos

Ajuste del modelo

Podado del árbol (pruning)

Predicción y evaluación del modelo

Árboles de clasificación

Idea intuitiva

Construcción del árbol

Predicción del árbol

Ejemplo clasificación

Librerías

Datos

Ajuste del modelo

Predicción y evaluación del modelo

Podado del árbol (pruning)

Random Forest

Idea intuitiva

Métodos de ensemble

Bagging

Entrenamiento de Random Forest

Predicción de Random Forest

Out-of-Bag Error

Importancia de los predictores

Ejemplo regresión

Librerías

Datos

Ajuste del modelo

Predicción y evaluación del modelo

Optimización de hiperparámetros

Número de árboles

mtry

Grid search

Grid Search basado en out-of-bag error

Grid Search basado en validación cruzada

Importancia de predictores

Importancia por pureza de nodos

Importancia por permutación

Ejemplo clasificación

Librerías

Datos

Ajuste del modelo y optimización de hiperparámetros

Grid Search basado en out-of-bag error

Grid Search basado en validación cruzada

Predicción y evaluación del modelo

Predicción de probabilidades

Importancia de predictores

Importancia por pureza de nodos

Importancia por permutación

Boosting

Idea intuitiva

AdaBoost

Gradient Boosting

Stochastic Gradient Boosting

Binning

Parada temprana (early stopping)

Importancia de los predictores

Ejemplo regresión

Librerías

Datos

Ajuste del modelo

Predicción y evaluación del modelo

Optimización de hiperparámetros

Número de árboles

Learning rate

Grid search

Grid Search basado en validación cruzada

Importancia de predictores