ÁRBOLES DE DECISIÓN Y MÉTODOS DE ENSEMBLE
Cristina Gil Martínez
Junio, 2018
INTRODUCCIÓN
Los métodos para regresión y clasificación basados en árboles de decisión estratifican o segmentan el espacio del predictor en un número simple de regiones, y para obtener las predicciones se suele usar la media o moda de las observaciones de entrenamiento en la región en la que cada observación a predecir pertenece. Los árboles de decisión son simples y fáciles de interpretar, pero pueden no resultar lo suficientemente competitivos frente a otros métodos de aprendizaje supervisado en cuanto a la precisión de predicción. Los métodos de bagging, random forests y boosting producen múltiples árboles que se combinan para mejorar la predicción, a expensar de una interpretación más complicada.
ÁRBOLES DE DECISIÓN
Los árboles de decisión pueden aplicarse tanto para problemas de regresión como de clasificación, y también pueden contener predictores tanto cuantitativos como cualitativos.
Árboles de regresión
Una de las ventajas de los árboles de decisión en comparación con otros métodos de regresión es su fácil interpretación y útil representación gráfica. A la hora de construir un árbol de regresión, dos son los pasos principales
1. División del espacio del predictor (conjunto de valores posibles para \(\small{X_1, X_2, ..., X_p}\) en J regiones distintas y no solapantes o nodos internos, \(\small{R_1, R_2, ..., R_j}\)). Teóricamente, las regiones podrían tomar cualquier forma, aunque si se opta por regiones rectangulares de múltiples dimensiones se simplifica y facilita la interpretación del modelo final. El objetivo es encontrar las regiones \(\small{R_1, ..., R_j}\) que minimicen el Residual Sum of Squares (RSS)
donde \(\small{y_{Rj}}\) es la media de la variable respuesta en la región j.
2. Para cada una de las observaciones pertenecientes a la región \(\small{R_j}\) se asigna la misma predicción, que es simplemente la media de la variable respuesta de las observaciones de entrenamiento en \(\small{R_j}\).
División binaria recursiva
Desafortunadamente, resulta computacionalmente inviable considerar todas las posibles particiones del predictor, por lo que se opta por una división binaria recursiva (recursive binary splitting), similar al concepto de stepwise selection, consiguiendo también obtener buenos resultados. El proceso comienza en el punto más alto del árbol (donde todas las observaciones perteneces a una misma región), y de manera sucesiva se divide el espacio del predictor, donde cada división genera dos nuevas ramas. Se optará además por la mejor división en cada paso, sin tener en cuenta si dicha división mejorará el árbol en futuros pasos o divisiones.
Para llevar a cabo la división binaria recursiva, lo primero es seleccionar el predictor \(\small{X_j}\) y el punto de corte s de manera que al dividir el espacio del predictor en las regiones {X|\(\small{X_j}\) < s} y {X|\(\small{X_j}\) ≥ s} se consiga la máxima reducción del RSS.
Es decir, se consideran todos los predictores \(\small{X, ..., X_p}\) y todos los posibles puntos de corte para cada uno de los predictores, eligiéndose el predictor \(\small{j}\) y punto de corte \(\small{s}\) que minimicen
correspondiente a la suma del RSS de la región 1 y el RSS de la región 2.
El proceso continúa, pero esta vez en lugar de dividir el espacio completo del predictor, se divide una de las dos regiones creadas en el paso anterior. El proceso continúa hasta que se alcance un criterio de parada, por ejemplo, hasta que ninguna región contenga más de \(\small{n}\) observaciones, hasta alcanzar un máximo de nodos terminales, etc. Una vez las regiones \(\small{R_1, ..., R_j}\) han sido creadas, predecimos la respuesta para una determinada observación de test como la media de las observaciones de entrenamiento en la región en la cual pertenece dicha observación.
Podado del árbol (Pruning)
El proceso de división binaria recursiva puede conseguir buenas predicciones con los datos de entrenamiento, ya que reduce el RSS de entrenamiento, lo que implica un sobreajuste a los datos (derivado de la facilidad de ramificación y posible complejidad del árbol resultante), reduciendo la capacidad predictiva para nuevos datos.
Una posible alternativa supone construir un árbol hasta el punto en el que la reducción del RSS debido a cada división supere un determinado límite (alto). Con esta estrategia obtendremos árboles más pequeños, con menos ramificaciones, con lo que conseguimos reducir la varianza y mejorar la interpretabilidad a expensas de una pequeña reducción del bias. Concretamente, la estrategia consistirá en construir un árbol grande \(\small{T_0}\) y luego podarlo para obtener un sub-árbol que consiga el menor test error, obtenido mediante validación cruzada o validación simple. Sin embargo, suele resultar muy costoso obtener el error de validación para cada posible sub-árbol debido al gran número de posibles sub-árboles, por lo que se opta por seleccionar un grupo pequeño de sub-árboles a considerar. Una manera de conseguir esto es mediante el cost complexity pruning, al que corresponde el siguiente algoritmo:
1. Aplicar el método de división binaria recursiva para obtener un árbol grande (\(\small{T_0}\)) con los datos de entrenamiento, finalizando el proceso con el cumplimiento de una determinada norma de parada.
2. Aplicar cost complexity pruning a \(\small{T_0}\) para obtener una secuencia de mejores sub-árboles, en función del parámetro de penalización o tuning parameter α, que controla el equilibrio entre la complejidad del árbol y su ajuste a los datos de entrenamiento. Conforme α aumenta, más ramas se podan del árbol. A cada valor de α le corresponde un sub-árbol \(\small{T ⊂ T_0}\) tal que se minimice
donde \(\small{|T|}\) es el número de nodos terminales u hojas del árbol \(\small{T}\), y \(\small{R_m}\) es el rectángulo correspondiente al nodo terminal \(\small{m}\). Siendo \(\small{α = 0}\), \(\small{T}\) será igual a \(\small{T_0}\).
3. Usar k-fold cross validation para escoger α. Esto consiste en dividir las observaciones de entrenamiento en \(\small{K}\) grupos de manera que para cada \(\small{k=1}\), …, \(\small{K}\):
Repetir pasos 1 y 2 en los \(\small{k}\) – 1 grupos.
Evaluar el error cuadrático medio (MSE) con el grupo \(\small{k_i}\) de observaciones, en función de α (conveniente graficarlo en función de \(\small{|T|}\)).
Promediar los resultados para cada valor de α, escogiendo el valor que minimice el test error.
4. Seleccionar el sub-árbol del paso 2 correspondiente al valor escogido de α por validación cruzada.
Árboles de clasificación
Los árboles de clasificación son muy similares a los de regresión, con la diferencia de que se usan para predecir una variable respuesta cualitativa, asignando la predicción para cada observación como la clase más común (moda) de observaciones de entrenamiento en la región o nodo terminal al que pertenece dicha observación de test.
Al igual que con los árboles de regresión, se aplica la división binaria recursiva para generar el árbol, pero el RSS no puede usarse como criterio para estas divisiones, sino otras alternativas como:
ERROR DE CLASIFICACIÓN
Fracción de observaciones de entrenamiento en una región o nodo que no pertenece a la clase más frecuente
donde \(\small{\hat{p}_{mk}}\) es la proporción de observaciones de entrenamiento en la región \(\small{m}\) que pertenecen a la clase \(\small{k}\).
El error de clasificación no suele ser lo suficientemente sensible en medir la pureza de los nodos para el crear el árbol (es más recomendable para hacer el pruning posterior), por lo que en la práctica se opta por otras dos medidas, el índice de Gini o el cross-entropy. Aun así, cualquiera de los tres puede usarse para el la poda del árbol. El error de clasificación es preferible si el objetivo es conseguir la máxima predicción en las predicciones del árbol podado final.
ÍNDICE DE GINI
El índice de Gini se define como una medida de la varianza total en el conjunto de las \(\small{K}\) clases, o pureza de los nodos
Este índice será pequeño cuando todos los \(\small{\hat{p}_{mk}}\) estén próximos a 0 o 1. Un valor bajo indica que el nodo contiene predominantemente observaciones de una sola clase.
CROSS-ENTROPY
El cross-entropy viene dado por
Al igual que el índice de Gini, el valor del cross-entropy será pequeño cuando todos los \(\small{\hat{p}_{mk}}\) estén próximos a 0 o 1. Un valor bajo también indica que el nodo contiene predominantemente observaciones de una sola clase.
Árboles de decisión vs modelos lineales
Mientras que la regresión lineal asume un modelo con la forma
los árboles de regresión asumen un modelo con la forma
donde \(\small{R_1, ..., R_M}\) representa la partición del espacio del predictor.
El modelo que pueda dar mejores resultados dependerá del problema en cuestión: si la relación entre los predictores y la variable respuesta es aproximadamente lineal, entonces un modelo de regresión lineal dará buenos resultados.
Si de lo contrario existe una relación altamente no lineal y compleja, los árboles de decisión pueden superar los métodos más clásicos. En ambos casos, el rendimiento de unos modelos u otros puede medirse mediante la estimación del test error obtenido por validación cruzada o validación simple, aunque otras consideraciones pueden tenerse en cuenta, como por ejemplo el nivel de interpretabilidad o visualización.
Ventajas y desventajas de los árboles de decisión
• Fácil interpretación
• Representación gráfica muy intuitiva y posible aun cuando hay más de 3 predictores
• Fácil manejo de predictores cualitativos sin la necesidad de crear variables dummy
• No hay necesidad de asumir a priori ninguna relación entre las variables, pudiendo introducir relaciones muy complejas entre las mismas
• Seleccionan los predictores automáticamente
• Capacidad predictiva superable por otros métodos de regresión/clasificación. Sin embargo, mediante la agregación de varios árboles de decisión, métodos como bagging, random forests y boosting, la capacidad predictiva de los árboles puede mejorarse sustancialmente.
MÉTODOS DE ENSEMBLE
Los métodos de bagging, random forests y boosting nos permiten mejorar sustancialmente el rendimiento predictivo de modelos basados en árboles, aunque con el inconveniente de una considerable reducción en la facilidad de interpretación del modelo final. Estos métodos también se conocen como métodos de ensemble o métodos combinados, que son los que utilizan múltiples algoritmos de aprendizaje para obtener predicciones que mejoren las que se podrían obtener por medio de cualquiera de los algoritmos individuales. Son aplicables a muchos métodos de aprendizaje estadísticos (no solo árboles de decisión) para regresión o clasificación.
Bagging
Los árboles de decisión sufren de alta varianza, lo que significa que, si dividiéramos al azar los datos de entrenamiento en dos grupos y ajustáramos un árbol de decisión a cada mitad, los resultados que obtendríamos podrían ser bastante diferentes. Por el contrario, un procedimiento o método con baja varianza dará resultados parecidos aun aplicándose sobre sets de datos distintos. El método de bagging o bootstrap aggregation es un procedimiento utilizado para reducir la varianza de un método de aprendizaje estadístico, usado muy frecuentemente con árboles de decisión.
REGRESIÓN:
Dado un set de \(\small{n}\) observaciones independientes \(\small{Z_1, ..., Z_n}\), cada una con varianza \(\small{σ^2}\), la varianza del promedio de las \(\small{Z}\) observaciones viene dado por \(\small{σ^2/n}\). En otras palabras, promediar un grupo de observaciones reduce la varianza. De esta forma, un método para reducir la varianza y aumentar por lo tanto la precisión de las predicciones de un método de aprendizaje estadístico implica tomar muchos sets de observaciones de la población, ajustar un modelo de predicción por separado para cada set, y promediar las predicciones resultantes. De manera más formal, calcularíamos \(\small{\hat{f^1}{(x)}}\),\(\small{\hat{f^2}{(x)}}\),…,\(\small{\hat{f^B}{(x)}}\) usando \(\small{B}\) sets de entrenamiento, promediándolos obteniendo
En la práctica, no se suele tener acceso a múltiples sets de observaciones, por lo que podemos optar en este caso por el bootstrap, tomando repetidamente muestras (\(\small{^*b}\)) de nuestro único set de datos de entrenamiento
Por tanto, la aplicación del bagging a árboles de regresión consiste en crear \(\small{B}\) árboles de regresión usando los \(\small{B}\) sets de entrenamiento generados por bootstrapping, promediando finalmente las predicciones resultantes. Estos árboles pueden crecer bastante ya que apenas se aplican restricciones, además de que no son podados. De esta manera cada árbol individual tiene alta varianza y poco bias, pero promediando los \(\small{B}\) árboles se contrarresta la varianza.
CLASIFICACIÓN:
En un problema de clasificación, existen varias posibilidades para aplicar bagging, pero la más simple es la siguiente: dada una observación de test, podemos obtener la clase predicha por cada uno de los \(\small{B}\) árboles, y escoger como predicción final para dicha observación la clase más común de entre las \(\small{B}\) predicciones (predicción de cada árbol).
El número de árboles (\(\small{B}\)) a crear no es un parámetro crítico a la hora de aplicar bagging. Ajustar un gran número de árboles no aumentará el riesgo de overfitting, por lo que usaremos un número lo suficientemente alto como para alcanzar la estabilización en la reducción del test error.
Estimación del error out-of-bag
Una manera directa de estimar el test error de un modelo al que se aplica bagging sin necesidad de aplicar la validación cruzada o simple: en promedio, cada árbol generado por bootstrapping usa en torno a 2/3 de las observaciones. El 1/3 restante de observaciones no usadas para ajustar cada árbol se conocen como out-of-bag (OOB). Por lo tanto, podemos obtener la predicción para la \(\small{i}\)-ésima observación de test usando cada uno de los árboles en los cuales dicha observación sea OOB. Esto resultará en \(\small{B/3}\) predicciones para cada observación, por lo que, para obtener un solo valor, se opta por una medida distinta según nos encontremos en un problema de regresión o clasificación:
REGRESIÓN:
Se promedian los \(\small{B/3}\) valores. La estimación del test error se corresponde al OOB MSE.
CLASIFICACIÓN:
Se escoge la clase mayoritaria, o se clasifica en base al promedio de los valores de probabilidad \(\small{P(Y = y │ X)}\). La estimación del test error se corresponde con el error de clasificación.
Con un número suficientemente elevado de árboles (\(\small{B}\)) el error OOB puede llegar a ser equivalente al error de validación leave-one-out. Además, el método basado en OOB para estimar el test error resulta conveniente cuando se aplica bagging en sets de datos grandes, para los cuales aplicar la validación cruzada sería computacionalmente muy costoso.
Importancia de las variables
Cuando se combinan múltiples árboles mediante bagging, no es posible representar gráficamente el modelo resultante mediante un árbol, y no es identificable de manera inmediata qué variables son las más importantes. Por tanto, bagging mejora la predicción del modelo a expensas de la pérdida de interpretabilidad. Aun así, un modo de poder identificar qué predictores tienen mayor importancia es
REGRESIÓN:
Cantidad total reducida del RSS como resultado de las divisiones sobre cada predictor, promediada sobre todos los \(\small{B}\) árboles. Un predictor importante será el que consiga una reducción promedio mayor del RSS.
CLASIFICACIÓN:
Mismo proceso que en el caso de árboles de regresión, pero empleando el índice de Gini.
Random forests
Supóngase que se cuenta con un set de datos cuenta con un predictor muy importante o influyente que destaca sobre el resto. En este caso, todos o casi todos los árboles generados por bagging usarán este predictor en la primera ramificación, por lo que acabarán siendo similares unos a otros, y las predicciones entre ellos estarán altamente correlacionadas. En este escenario, la aplicación de bagging promediando valores correlacionados no consigue una reducción sustancial de la varianza con respecto a un solo árbol.
El método de random forests proporciona una mejora a los árboles combinados por bagging en cuanto a que los decorrelaciona, teniendo en cuenta solo un subgrupo de predictores en cada división. Al igual que en el bagging, se construyen un número de árboles de decisión a partir de pseudo-muestras generadas por bootstrapping. Esta vez, se escogen de entre todos los p predictores una muestra aleatoria de m predictores como candidatos antes de cada división, generalmente \(\small{m ≈ √p}\) (si \(\small{m=p}\), bagging y random forests darían resultados equivalentes). Solo se aplica la división a uno de los m predictores. Esto hace que, de media, \(\small{(p – m) / p}\) divisiones no tengan en cuenta el predictor más influyente, dando más oportunidades al resto.
Usar un número pequeño de \(\small{m}\) para aplicar random forests puede resultar útil cuando contamos con un gran número de predictores correlacionados.
Al igual que con bagging, aumentar el número de pseudo-árboles no incrementará el riesgo de overfitting, por lo que en la práctica usamos un valor lo suficientemente alto para conseguir una estabilización del test error.
Boosting
Boosting funciona de manera parecida al bagging en cuanto a que combina un gran número de árboles, a excepción de que los árboles se construyen de manera secuencial: cada árbol se genera usando información, concretamente los residuos, de árboles previamente generados, en lugar de utilizar la variable respuesta (por ello suelen ser suficientes árboles más pequeños, en lugar de un gran árbol que pueda sobreajustarse a los datos). Otra diferencia es que boosting no utiliza remuestreo por bootstrapping, sino que cada árbol se genera utilizando una versión modificada del set de datos original. Algunos parámetros principales del algoritmo son:
• Número de árboles (\(\small{B}\)). A diferencia del bagging y random forests, boosting puede sobreajustarse a los datos si el número de árboles es demasiado alto (aunque en caso de darse, aparece lentamente). \(\small{B}\) se selecciona por validación cruzada.
• Número de divisiones (\(\small{d}\)) en cada árbol, que controla el nivel de complejidad. Un valor de \(\small{d = 1}\) (cada árbol contiene una única división, es decir, un único predictor) suele dar buenos resultados.
• Parámetro de penalización (\(\small{λ}\)), que controla el ritmo con el que boosting aprende. Valores comunes para este parámetro suelen ser 0,01 o 0,001, aunque la decisión depende del problema en cuestión. Por ejemplo, un valor muy pequeño de \(\small{λ}\) puede requerir un número elevado de árboles para conseguir buenos resultados.
Algunos de los algoritmos de boosting más utilizados son:
AdaBoost
Gradient Boosting
Stochastic Gradient Boosting
AdaBoost
Adaboost (Adaptive Boosting), aplicado sobretodo a problemas de clasificación, es un algoritmo iterativo que se basa en combinar múltiples weak learners, en un único strong learner a través de una combinación lineal ponderada. En cada iteración se genera un weak learner con un peso asociado. En un escenario de clasificación binaria, primeramente:
Se codifican los niveles de la variable respuesta categórica como +1 y -1.
Se asocia un mismo peso inicial (\(\small{w_i}\)) a todas las observaciones \(\small{i}\) de entrenemiento.
Se elige el tipo de modelo para generar los weak learners.
El proceso iterativo seguido por el algoritmo es el siguiente:
Se comienza con la primera de las \(\small{M}\) iteraciones ajustando un weak learner \(\small{G_m(x)}\) con todas las observaciones de entrenamiento y sus pesos iniciales (\(\small{w_i}\)). Se predicen dichas observaciones con el weak learner obtenido para determinar las observaciones correcta e incorrectamente clasificadas, calculando el error cometido \(\small{err_m}\). Se asocia un peso \(\small{α_m}\) al weak learner en base a sus aciertos: cuantos más aciertos, mayor peso. Además, se actualizan los pesos \(\small{w_i}\) de las observaciones: se disminuye el peso de las bien clasificadas y se aumenta el de las mal clasificadas.
En la segunda iteración, se ajusta de nuevo el weak learner a los datos de entrenamiento, pero utilizando los pesos actualizados. La finalidad es poder predecir las observaciones que se han clasificado mal en la iteración anterior. El proceso continúa con \(\small{M}\) iteraciones con las que se obtienen \(\small{M}\) weak learners que formarán el ensemble final \(\small{G(x)}\). Las clasificaciones de nuevas observaciones se obtienen mediante la agregación ponderada (en función del peso individual de cada weak learner) de los resultados de predicción de los \(\small{M}\) weak learners.
Gradient Boosting
A diferencia de AdaBoost, Gradient Boosting no asigna un peso independiente a cada observación de entrenamiento, sino que hace uso de una función de coste \(\small{L(y_i, f(x))}\) cuyo gradiente o derivada parcial de la función de coste (obtenido a partir de los residuos \(\small{r_{im}}\)) se pretende minimizar. Esto se lleva a cabo en un proceso iterativo. El gradiente se utiliza para encontrar la dirección en la que cambiar los parámetros de los weak learners para reducir el error de predicción en las siguientes iteraciones. Concretamente, las predicciones del weak learner \(\small{m}\) intentan acercarse al gradiente negativo de la función de coste.
Suele iniciarse con la mejor aproximación de la variable respuesta (como la media en un problema de regresión), \(\small{f_0(x)}\): árbol con un solo nodo terminal. Se permite el uso de cualquier función de coste diferenciable (que admita derivadas en cualquier dirección).
En un problema de regresión, el algoritmo es:
Para problemas de clasificación el algoritmo es similar. Los pasos 2(a-d) se repiten \(\small{K}\) veces en cada iteración \(\small{m}\), una por cada clase. El el paso 3 se obtiene un acoplamiento de los \(\small{K}\) clasificadores distintos (\(\small{f_{kM}(x), k=1,2,...,K}\))
Gradientes para funciones de coste comunes:
Así pues, como hiperparámetros relacionados con el descenso de gradiente, tendremos:
loss function: función de coste
learning rate: tamaño de cada paso en el descenso de gradiente. A menor valor, más iteraciones pueden ser requeridas
shrinkage: reducción del learning rate
Stochastic Gradient Boosting
Stochastic Gradient Boosting supone una variación del Gradient Boosting en el sentido que cada iteración hace uso de un subconjunto aleatorio y sin reemplazo de los datos de entrenamiento disponibles para ajustar los weak learners. Esto, además de proporcionar un extra de aleatoreidad, permite obtener el error out-of-bag, lo que permite no tener que recurrir a la validación cruzada para la optimización de hiperparámetros cuando los requerimientos computacionales son limitantes.
EJEMPLOS EN R
Árboles simples: library(tree)
• tree() -> Ajusta árboles de decisión mediante partición binaria recursiva.
• tree.control() -> Selección de parámetros para el árbol. Se aplica al argumento control de la función tree().
Poda:
• cv.tree() -> Lleva a cabo k-fold cross validation para estimar la deviance o error de clasificación en función del parámetro cost-complexity k. Especificar argumento FUN = prune.tree para árboles de regresión y FUN = prune.misclass para árboles de clasificación.
• prune.tree() -> Genera sub-árboles a partir del árbol original mediante cost- complexity pruning.
Bagging y random forests: library(RandomForest)
• randomForest() -> Implementación del algoritmo de bagging (m = p) y random forests (m < p), para regresión y clasificación. Con el argumento mtry se especifica el número de variables, y con ello el método a seguir.
• importance() -> Calcula dos medidas (%IncMSE y IncNodePurity) que identifican la importancia de los predictores del modelo generado por la función ranfomForest().
• varImpPlot() -> Diagrama de representación de la importancia de las variables.
Boosting y AdaBoost
• gbm() -> Implementación del algoritmo de Stochastic Gradient Boosting. Especificar argumento distribution = “gaussian” para árboles de regresión, distribution = “Bernoulli” para árboles de clasificación binaria, distribution = adaboost para implementar AdaBoost para clasiicación binaria. El argumento interaction.depth limita la profundidad (complejidad) de cada árbol. El número de árboles ntrees son por defecto 5000, y el parámetro de penalización λ = 0,01.
Ejemplos de regresión
En este ejemplo utilizaremos el set de datos Hitters del paquete ISLR, que contiene información sobre jugadores de beisbol de la Major League Baseball (MLB) de las temporadas 1986-1987. Aplicando todos los métodos descritos en la sección teórica, intentaremos predecir los salarios de los jugadores.
NOTA: Para ajustar un árbol de regresión, la variable respuesta ha de ser de tipo numeric.
library(ISLR)
str(Hitters)## 'data.frame': 322 obs. of 20 variables:
## $ AtBat : int 293 315 479 496 321 594 185 298 323 401 ...
## $ Hits : int 66 81 130 141 87 169 37 73 81 92 ...
## $ HmRun : int 1 7 18 20 10 4 1 0 6 17 ...
## $ Runs : int 30 24 66 65 39 74 23 24 26 49 ...
## $ RBI : int 29 38 72 78 42 51 8 24 32 66 ...
## $ Walks : int 14 39 76 37 30 35 21 7 8 65 ...
## $ Years : int 1 14 3 11 2 11 2 3 2 13 ...
## $ CAtBat : int 293 3449 1624 5628 396 4408 214 509 341 5206 ...
## $ CHits : int 66 835 457 1575 101 1133 42 108 86 1332 ...
## $ CHmRun : int 1 69 63 225 12 19 1 0 6 253 ...
## $ CRuns : int 30 321 224 828 48 501 30 41 32 784 ...
## $ CRBI : int 29 414 266 838 46 336 9 37 34 890 ...
## $ CWalks : int 14 375 263 354 33 194 24 12 8 866 ...
## $ League : Factor w/ 2 levels "A","N": 1 2 1 2 2 1 2 1 2 1 ...
## $ Division : Factor w/ 2 levels "E","W": 1 2 2 1 1 2 1 2 2 1 ...
## $ PutOuts : int 446 632 880 200 805 282 76 121 143 0 ...
## $ Assists : int 33 43 82 11 40 421 127 283 290 0 ...
## $ Errors : int 20 10 14 3 4 25 7 9 19 0 ...
## $ Salary : num NA 475 480 500 91.5 750 70 100 75 1100 ...
## $ NewLeague: Factor w/ 2 levels "A","N": 1 2 1 2 2 1 1 1 2 1 ...Primero eliminaremos las observaciones para la que no hay información sobre la variable respuesta Salary. Hay otras técnicas para manejar los valores ausentes, pero para este ejemplo se decide eliminarlos:
sum(is.na(Hitters$Salary))## [1] 59# Omisión de NAs y dimensión final de los datos
datos_Hitters <- na.omit(Hitters)
dim(datos_Hitters)## [1] 263 20# Distribución variable respuesta
library(ggplot2)
ggplot(data = datos_Hitters, aes(x = Salary)) +
geom_histogram(color = "white", fill = "black") +
labs(title = "Distribución Salary") +
theme(plot.title = element_text(hjust = 0.5))
La variable respuesta posee una distribución sesgada. Aplicaremos una transformación logarítmica para hacer su distribución más normal.
# Transformación logaritmica de la variable respuesta Salary
datos_Hitters$Salary <- log(datos_Hitters$Salary)
ggplot(data = datos_Hitters, aes(x = Salary)) +
geom_histogram(color = "white", fill = "black") +
labs(title = "Distribución log(Salary)") +
theme(plot.title = element_text(hjust = 0.5))
Antes de pasar a generar los modelos, dividimos el set de datos en un grupo de entrenamiento (para el ajuste de los modelos) y otro de test (para la evaluación de los mismos). Esta división dependerá de la cantidad de observaciones con las que contemos y la seguridad con la que queramos obtener la estimación del test error.
# Índice observaciones de entrenamiento
train <- 1:200
# Datos entrenamiento
datosH_train <- datos_Hitters[train, ]
# Datos test
datosH_test <- datos_Hitters[-train, ]Árbol de regresión simple
Ajuste del modelo
Por defecto, las condiciones de stop que regulan el crecimiento del árbol, y que hacen que no todas las variables incluidas en la fórmula tengan por qué ser utilizadas en el árbol final, son:
• mincut = 5: Cantidad ponderada del nº de observaciones mínimas de al menos uno de los nodos hijos resultantes para permitir la división.
• minsize = 10: Cantidad ponderada del nº de observaciones mínimas en un nodo para permitir su división.
• mindev = 0,01: Deviance mínima para proceder con una división. Para obtener un árbol más pequeño, reducir el valor.
Estos valores pueden ser cambiados por el usuario. A continuación, se muestra el ajuste, pero usando los valores por defecto:
library(tree)
# Selección de parámetros para el árbol
setup <- tree.control(nobs = nrow(datosH_train),
mincut = 5,
minsize = 10,
mindev = 0.01)
# Ajuste del árbol de regresión
modelo_arbolR <- tree(Salary ~ ., data = datosH_train,
split = "deviance",
control = setup)
summary(modelo_arbolR)##
## Regression tree:
## tree(formula = Salary ~ ., data = datosH_train, control = setup,
## split = "deviance")
## Variables actually used in tree construction:
## [1] "CAtBat" "CRuns" "AtBat" "Walks" "Hits" "CHits" "CRBI"
## [8] "PutOuts"
## Number of terminal nodes: 10
## Residual mean deviance: 0.1665 = 31.64 / 190
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.003000 -0.221300 0.009429 0.000000 0.246800 2.342000Según la información devuelta, del conjunto de las 19 variables, el modelo ha utilizado para generar el árbol solo 8 (nodos internos): CAtBat, CRuns, AtBat, Walks, Hits, CHits, CRBI y PutOuts. Con estas variables el número de nodos terminales es 10. El residual mean deviance se corresponde con el training RSS (31,64) dividido entre el nº de observaciones – nº de nodos terminales (200 – 10 = 190). Cuanto menor es este valor, mejor se ajusta el modelo a los datos de entrenamiento.
plot(modelo_arbolR, type = "proportional")
# pretty = 0 incluye los nombres de los niveles para las variables cualitativas, en lugar de mostrar solo una letra
text(modelo_arbolR, splits = TRUE, pretty = 0, cex = 0.8, col = "firebrick")
La variable que se encuentra en lo más alto del árbol (primer nodo), y por tanto, la más importante, es CAtBat. Esta primera rama diferencia a los jugadores según el número de bateos a lo largo de la trayectoria profesional del jugador. En este caso, los que superan 1322 bateos tienen un salario más alto. De los que tienen menor salario, la variable CRuns es más determinante como segunda variable más importante que en el grupo de jugadores que ganan más, donde Walks juega un papel más determinante a este nivel. Además, la altura de las ramas indica la efectividad de la división en cuanto a reducción del RSS.
Si además inspeccionamos el objeto tree, R nos muestra información correspondiente a cada rama del árbol: criterio de división, número de observaciones en cada rama antes de la división (n), el error o deviance (RSS), y la media de la predicción final para cada rama (yval). Con asteriscos se indican los nodos terminales.
modelo_arbolR## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 200 166.4000 5.940
## 2) CAtBat < 1322 70 23.3000 4.971
## 4) CRuns < 92.5 43 12.6000 4.696
## 8) AtBat < 173 5 7.2760 5.321 *
## 9) AtBat > 173 38 3.1220 4.614 *
## 5) CRuns > 92.5 27 2.2610 5.409 *
## 3) CAtBat > 1322 130 42.0400 6.462
## 6) Walks < 48.5 80 20.2600 6.232
## 12) Hits < 98.5 41 9.8330 6.019
## 24) AtBat < 335 36 6.8960 6.097 *
## 25) AtBat > 335 5 1.1580 5.461 *
## 13) Hits > 98.5 39 6.6220 6.456
## 26) CHits < 590.5 12 1.8470 6.118 *
## 27) CHits > 590.5 27 2.7950 6.606 *
## 7) Walks > 48.5 50 10.8100 6.829
## 14) CRBI < 369.5 16 0.9873 6.498 *
## 15) CRBI > 369.5 34 7.2280 6.985
## 30) PutOuts < 286 20 2.6210 6.786 *
## 31) PutOuts > 286 14 2.6770 7.270 *Por ejemplo, el primer nodo contiene todas las observaciones (200), con un RSS de 166,4 y una predicción promedio del salario de 5,94 en unidades log Salary. La primera división (\(\small{R_1}\)) está definida por la condición CAtBat < 1322, habiendo 70 jugadores que cumplen esta condición (sub-árbol por debajo de este nodo). El RSS para este grupo es de 23,3, siendo el salario medio de estos jugadores de \(\small{e^{4,971} = 144,17}\) dólares.
Evaluación del modelo
# Predicciones del modelo sobre los datos de test
predicciones <- predict(modelo_arbolR, newdata = datosH_test)
plot(x = predicciones, y = datosH_test$Salary,
main = "Predicciones modelo vs valor real",
xlab = "Predicción",
ylab = "Salary real",
col = "darkgrey", pch = 19)
abline(a = 0, b = 1, col = "blue")
# Error MSE en test
testMSE <- mean((predicciones - datosH_test$Salary)^2)
testMSE## [1] 0.3116231Poda
Ajuste del modelo
El siguiente paso es considerar si la reducción de la varianza y complejidad mediante la poda del árbol (mediante cost complexity pruning) puede mejorar las predicciones del modelo. Para ello hacemos uso de la función cv.tree() para encontrar le parámetro de penalización α óptimo por validación cruzada. Por defecto la función utiliza la deviance o RSS (FUN = prune.tree) como método para guiar el proceso de validación cruzada y poda (puede cambiarse con el argumento FUN ).
# 10-fold cross validation
set.seed(356)
cv_arboles <- cv.tree(modelo_arbolR, K = 10, FUN = prune.tree)
cv_arboles## $size
## [1] 10 9 8 7 6 5 4 3 2 1
##
## $dev
## [1] 64.97950 64.21269 63.76574 62.63029 61.24361 59.65822 57.95912
## [8] 62.48435 70.11162 166.80136
##
## $k
## [1] -Inf 1.779047 1.930350 1.980293 2.206202 2.590201
## [7] 3.801387 8.432809 10.977019 101.065841
##
## $method
## [1] "deviance"
##
## attr(,"class")
## [1] "prune" "tree.sequence"Entre la información devuelta por el objeto cv.tree se encuentra: número de nodos terminales de cada árbol considerado (size), con su correspondiente error de validación (dev) y parámetro de penalización (k). Con esta información, podemos representar el error de validación (deviance) en función del tamaño del árbol y el parámetro de penalización para identificar el mejor sub-árbol:
par(mfrow = c(1, 2))
# Cost complexity pruning
plot(cv_arboles$size, cv_arboles$dev, xlab = "nodos terminales", ylab = "RSS",
type = "b", pch = 19)
plot(cv_arboles$k, cv_arboles$dev, xlab = "alpha", ylab = "RSS", type = "b")
n_mejores_nodos <- cv_arboles$size[which.min(cv_arboles$dev)]
n_mejores_nodos## [1] 4El número de nodos terminales que consigue minimizar el error de validación es 4, por lo que se reducen a la mitad con respecto al árbol inicial.
NOTA: Puede haber casos en los que la poda no consiga mejorar el modelo, escogiéndose por validación cruzada el mayor tamaño de modelo.
# Poda del arbol
modelo_arbolRpodado <- prune.tree(tree = modelo_arbolR, best = n_mejores_nodos)
par(mfrow = c(1,1))
plot(x = modelo_arbolRpodado, type = "proportional")
text(modelo_arbolRpodado, splits = TRUE, pretty = 0, cex = 0.8, col = "firebrick")
Evaluación del modelo
Por último, hacemos uso del set de datos de test para evaluar la capacidad predictiva del sub-árbol obtenido.
# Predicciones del modelo sobre los datos de test
predicciones <- predict(modelo_arbolRpodado, newdata = datosH_test)
plot(x = predicciones, y = datosH_test$Salary,
main = "Predicciones modelo pruned vs valor real",
xlab = "Predicción",
ylab = "Salary real",
col = "darkgrey",
pch = 19)
abline(a = 0, b = 1, col = "blue")
# Error MSE con datos de test
testMSE <- mean((predicciones - datosH_test$Salary)^2)
testMSE## [1] 0.3978137Bagging
Ajuste del modelo
Para intentar mejorar la capacidad predictiva del árbol obtenido en el paso inicial, aplicaremos el método de bagging. El set de datos cuenta con 20 variables, 19 predictores y 1 variable respuesta. Por tanto, para aplicar bagging todos los predictores han de considerarse en cada partición. El número de pseudo-árboles que se generan por defecto son 500:
library(randomForest)
# Bagging
set.seed(356)
modelo_bagging <- randomForest(Salary ~ ., data = datosH_train,
mtry = 19,
importance = TRUE, #evaluar importancia predictores
ntree = 500)
modelo_bagging##
## Call:
## randomForest(formula = Salary ~ ., data = datosH_train, mtry = 19, importance = TRUE, ntree = 500)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 19
##
## Mean of squared residuals: 0.218353
## % Var explained: 73.76El porcentaje de la varianza explicada por el modelo aplicando bagging es del 73,76%. El Mean of squared residuals se corresponde con el out-of-bag-MSE, que puede entenderse como un test error. Su evolución respecto al número de árboles puede representarse de la siguiente manera:
plot(modelo_bagging, col = "firebrick")
El error se estabiliza al alcanzar los 100 árboles aproximadamente.
Importancia de variables
Al combinar múltiples árboles perdemos la posibilidad de obtener una representación gráfica como en el caso de un árbol único. Sin embargo, podemos utilizar la función importance() para obtener la importancia de cada predictor:
# Importancia de las variables en el modelo
importance(modelo_bagging)## %IncMSE IncNodePurity
## AtBat 9.8025542 7.39507746
## Hits 7.0348420 4.12897254
## HmRun 0.7290035 1.58618258
## Runs 5.0863829 2.88099141
## RBI 2.3253375 3.30797998
## Walks 9.6632238 7.27289997
## Years 12.3307475 1.88681219
## CAtBat 35.0961363 84.18605178
## CHits 7.9038625 12.39389250
## CHmRun 9.6036059 4.73271172
## CRuns 12.2386029 13.29016163
## CRBI 9.1685615 9.23786396
## CWalks 5.7962703 5.87407636
## League -0.6993126 0.09175593
## Division -3.1237905 0.17006076
## PutOuts 1.8834923 2.76595582
## Assists -0.1377997 1.52320320
## Errors 0.2149212 1.27567210
## NewLeague 0.6814064 0.19992497Para cada predictor se devuelven dos valores:
• %IncMSE: disminución media de la precisión de las predicciones sobre las muestras OOB cuando la variable dada se excluye del modelo.
• IncNodePurity: medida de la disminución total de impureza de los nodos (medida por el training RSS) que resulta de la división de la variable dada.
Los predictores más importantes se corresponderán a aquellos con mayor %IncMSE y IncNodePurity.
varImpPlot(modelo_bagging)
En este ejemplo la variable más importante es CAtBat, al igual que en el árbol simple.
Evaluación del modelo
# Predicciones del modelo sobre los datos de test
predicciones <- predict(modelo_bagging, newdata = datosH_test)
plot(x = predicciones, y = datosH_test$Salary,
main = "Predicciones modelo vs valor real",
xlab = "predicciones",
col = "darkgrey",
pch = 19)
abline(a = 0, b = 1, col = "blue")
# Error MSE en test
testMSE <- mean((predicciones - datosH_test$Salary)^2)
testMSE## [1] 0.2307334Random forests
El método de random forests se aplica de la misma forma que anteriormente con bagging, pero especificando en el argumento mtry un número menor de predictores respecto al total (\(\small{m<p}\)). Con este método, la función randomForest() usa por defecto \(\small{p/3}\) variables para árboles de regresión, y \(\small{√p}\) con árboles de clasificación. En lugar de utilizar \(\small{p/3}\) en este caso, es mejor determinar el valor óptimo en función de un error de validación. Con el paquete caret podemos llevar a cabo esta búsqueda. Contamos con \(\small{p=19}\) variables, por lo que como máximo tenemos que evaluar \(\small{p-1}\) para que el proceso sea considerado random forest.
Ajuste del modelo
library(caret)
library(parallel)
library(doParallel)
cluster <- makeCluster(detectCores() - 1)
registerDoParallel(cluster) # procesamiento paralelo
# Método de validación cruzada (10-fold)
fitControl <- trainControl(method = "cv",
number = 10,
search = "grid",
allowParallel = TRUE)
# Hiperparámetro a optimizar: número de predictores aleatorios en cada ramificación.
grid_mtry <- expand.grid(mtry = c(2:18))
# Ajuste del modelo random forest
set.seed(356)
modelo_rf <- train(Salary ~ ., data = datosH_train,
method = "rf",
metric = "RMSE",
ntree = 500,
tuneGrid = grid_mtry,
trControl = fitControl)
modelo_rf## Random Forest
##
## 200 samples
## 19 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 180, 180, 180, 180, 179, 180, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 2 0.4343907 0.7721497 0.3118140
## 3 0.4327753 0.7724073 0.3073994
## 4 0.4358981 0.7690836 0.3119219
## 5 0.4349379 0.7694211 0.3090616
## 6 0.4370962 0.7678961 0.3118758
## 7 0.4364547 0.7683587 0.3101804
## 8 0.4374803 0.7673074 0.3103542
## 9 0.4401669 0.7652127 0.3140100
## 10 0.4400828 0.7642849 0.3103472
## 11 0.4416511 0.7637998 0.3132640
## 12 0.4426732 0.7620780 0.3156464
## 13 0.4414202 0.7629171 0.3129724
## 14 0.4422273 0.7628450 0.3137601
## 15 0.4414179 0.7632810 0.3134858
## 16 0.4455226 0.7600228 0.3176640
## 17 0.4446994 0.7606683 0.3158113
## 18 0.4438593 0.7615774 0.3162322
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 3.plot(modelo_rf)
El valor más óptimo de mtry en cuanto a menor RMSE ha sido de 3 predictores.
Importancia de variables
# Top 10 variables mas importantes en el modelo
plot(varImp(modelo_rf, scale = FALSE), top = 10) 
En este caso, las dos variables con mayor importancia son CHits y CAtBat.
Evaluación del modelo
# Predicciones del modelo sobre los datos de test
predicciones <- predict(modelo_rf, newdata = datosH_test)
plot(x = predicciones, y = datosH_test$Salary,
main = "Predicciones modelo_rf vs valor real",
xlab = "Predicción",
ylab = "Salary real",
col = "darkgrey",
pch = 19)
abline(a = 0, b = 1, col = "blue")
testMSE <- mean((predicciones - datosH_test$Salary)^2)
testMSE## [1] 0.2103545El test error se ha reducido levemente respecto al modelo ajustado por bagging. Por lo general, random forest ofrece mejores resultados que bagging al decorrelacionar los árboles.
Stochastic Gradient Boosting
Ajuste del modelo
De entre los hiperparámetros principales, tenemos:
shrinkage: tasa de aprendizajen.trees: número de iteraciones (weak learners)n.minobsinnnode: nº mínimo de observaciones en nodo para su divisióninteraction.depth: complejidad del árbol
# Método de validación cruzada (10-fold)
fitControl <- trainControl(method = "cv",
number = 10,
search = "grid",
allowParallel = TRUE)
# Combinación de hiperparámetro a evaluar: 5*3*4*3*=180
grid_hiperparametros <- expand.grid(shrinkage = c(0.002, 0.005, 0.01, 0.012, 0.02),
n.trees = c(200, 300, 400),
n.minobsinnode = c(5, 7, 10, 20),
interaction.depth = c(1, 2, 3))
# Ajustamos el modelo obteniendo la combinacion de mejores hiperparametros
set.seed(356)
modelo_gboost <- train(Salary ~ ., data = datosH_train,
method = "gbm",
metric = "RMSE",
tuneGrid = grid_hiperparametros,
trControl = fitControl,
verbose = FALSE)
# Evolución del RMSE en función de los hiperparámetros
plot(modelo_gboost)
Importancia de variables
library(erboost)
# Importancia de las variables en el modelo (top 10)
plot(varImp(modelo_gboost, scale = FALSE), top = 10)
CAtBat es, con diferencia, el predictor más influyente del modelo ajustado por boosting, seguido de CRBI.
Evaluación del modelo
# Predicciones del modelo sobre los datos de test
predicciones <- predict(modelo_gboost, newdata = datosH_test)
plot(x = predicciones, y = datosH_test$Salary,
main = "Predicciones modelo boosting vs valor real",
xlab = "Predicción",
ylab = "Salary real",
col = "darkgrey", pch = 19)
abline(a = 0, b = 1, col = "blue")
# Error MSE en test
testMSE <- mean((predicciones - datosH_test$Salary)^2)
testMSE## [1] 0.2853416Ejemplos de clasificación
Para este ejemplo de clasificación utilizaremos el set de datos OJ, también parte del paquete ISLR. Contiene información sobre compra de dos tipos de bebida (Citrus Hill y Minute Maid Orange Juice) por parte de 1070 clientes (las variables registran distintas características del cliente y el producto). En este caso, generaremos árboles de clasificación que predigan que tipo de bebida (Purchase) se compra, en función del conjunto de predictores.
str(OJ)## 'data.frame': 1070 obs. of 18 variables:
## $ Purchase : Factor w/ 2 levels "CH","MM": 1 1 1 2 1 1 1 1 1 1 ...
## $ WeekofPurchase: num 237 239 245 227 228 230 232 234 235 238 ...
## $ StoreID : num 1 1 1 1 7 7 7 7 7 7 ...
## $ PriceCH : num 1.75 1.75 1.86 1.69 1.69 1.69 1.69 1.75 1.75 1.75 ...
## $ PriceMM : num 1.99 1.99 2.09 1.69 1.69 1.99 1.99 1.99 1.99 1.99 ...
## $ DiscCH : num 0 0 0.17 0 0 0 0 0 0 0 ...
## $ DiscMM : num 0 0.3 0 0 0 0 0.4 0.4 0.4 0.4 ...
## $ SpecialCH : num 0 0 0 0 0 0 1 1 0 0 ...
## $ SpecialMM : num 0 1 0 0 0 1 1 0 0 0 ...
## $ LoyalCH : num 0.5 0.6 0.68 0.4 0.957 ...
## $ SalePriceMM : num 1.99 1.69 2.09 1.69 1.69 1.99 1.59 1.59 1.59 1.59 ...
## $ SalePriceCH : num 1.75 1.75 1.69 1.69 1.69 1.69 1.69 1.75 1.75 1.75 ...
## $ PriceDiff : num 0.24 -0.06 0.4 0 0 0.3 -0.1 -0.16 -0.16 -0.16 ...
## $ Store7 : Factor w/ 2 levels "No","Yes": 1 1 1 1 2 2 2 2 2 2 ...
## $ PctDiscMM : num 0 0.151 0 0 0 ...
## $ PctDiscCH : num 0 0 0.0914 0 0 ...
## $ ListPriceDiff : num 0.24 0.24 0.23 0 0 0.3 0.3 0.24 0.24 0.24 ...
## $ STORE : num 1 1 1 1 0 0 0 0 0 0 ...sum(is.na(OJ$Purchase))## [1] 0El set de datos no contiene observaciones ausentes para la variable respuesta Purchase.
# Distribución variable respuesta
library(ggplot2)
ggplot(data = OJ, aes(x = Purchase, y = ..count.., fill = Purchase)) +
geom_bar() +
labs(title = "Distribución Purchase") +
scale_fill_manual(values = c("darkgreen", "orangered2"),
labels = c("Citrus Hill", "Orange Juice")) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
# Tabla frecuencias variable respuesta
table(OJ$Purchase)##
## CH MM
## 653 417# Tabla proporciones variable respuesta
library(dplyr)
prop.table(table(OJ$Purchase)) %>% round(digits = 2)##
## CH MM
## 0.61 0.39La clase mayoritaria (moda) en este caso es la bebida CH con el 61% de las compras. Este será el nivel basal a superar por el modelo (este es el porcentaje mínimo de aciertos si siempre se predijera CH).
Esta división dependerá de la cantidad de observaciones con las que contemos y la seguridad con la que queramos obtener la estimación del test error. Antes de pasar a generar los modelos, dividimos el set de datos en un grupo de entrenamiento (para el ajuste de los modelos) y otro de test (para la evaluación de los mismos). En este caso se opta por una división 80%-20%.
# Índices observaciones de entrenamiento
set.seed(524)
train <- createDataPartition(y = OJ$Purchase, p = 0.8, list = FALSE, times = 1)
# Datos entrenamiento
datosOJ_train <- OJ[train, ]
dim(datosOJ_train)## [1] 857 18# Datos test
datosOJ_test <- OJ[-train, ]
dim(datosOJ_test)## [1] 213 18Árbol de clasificación simple
Ajuste del modelo
NOTA: Para ajustar un árbol de clasificación, la variable respuesta tiene que ser de tipo factor.
# Modelo árbol de clasificación
modelo_arbolC <- tree(Purchase ~ ., data = datosOJ_train)
summary(modelo_arbolC)##
## Classification tree:
## tree(formula = Purchase ~ ., data = datosOJ_train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "SalePriceMM" "SpecialCH" "PriceDiff"
## Number of terminal nodes: 8
## Residual mean deviance: 0.7314 = 621 / 849
## Misclassification error rate: 0.1669 = 143 / 857El modelo ha utilizado solo 4 variables del total de 17, con un número final de 8 nodos terminales. El training error (observaciones mal clasificadas) es del 16,69%. El residual mean deviance se corresponde con la deviance (621) dividida entre el nº de observaciones – nº de nodos terminales (857 – 18 = 849). Para árboles de clasificación, la deviance (índice Gini o entropy) viene dada por
donde \(\small{n_{mk}}\) es el número de observaciones en el nodo terminal que pertenece a la clase \(\small{k_i}\).
# Representación gráfica del árbol
plot(modelo_arbolC)
# pretty = 0 incluye los nombres de los niveles para las variables cualitativas, en lugar de mostrar solo una letra
text(modelo_arbolC, splits = TRUE, pretty = 0, cex = 0.8, col = "firebrick")
modelo_arbolC## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 857 1146.00 CH ( 0.61027 0.38973 )
## 2) LoyalCH < 0.48285 319 343.20 MM ( 0.22884 0.77116 )
## 4) LoyalCH < 0.0616725 65 0.00 MM ( 0.00000 1.00000 ) *
## 5) LoyalCH > 0.0616725 254 304.70 MM ( 0.28740 0.71260 )
## 10) SalePriceMM < 2.04 143 135.60 MM ( 0.18182 0.81818 )
## 20) SpecialCH < 0.5 122 94.81 MM ( 0.13115 0.86885 ) *
## 21) SpecialCH > 0.5 21 29.06 MM ( 0.47619 0.52381 ) *
## 11) SalePriceMM > 2.04 111 151.30 MM ( 0.42342 0.57658 ) *
## 3) LoyalCH > 0.48285 538 479.40 CH ( 0.83643 0.16357 )
## 6) LoyalCH < 0.764572 255 314.00 CH ( 0.69412 0.30588 )
## 12) PriceDiff < -0.165 40 47.05 MM ( 0.27500 0.72500 ) *
## 13) PriceDiff > -0.165 215 230.80 CH ( 0.77209 0.22791 )
## 26) PriceDiff < 0.265 120 152.80 CH ( 0.66667 0.33333 ) *
## 27) PriceDiff > 0.265 95 59.54 CH ( 0.90526 0.09474 ) *
## 7) LoyalCH > 0.764572 283 86.50 CH ( 0.96466 0.03534 ) *El predictor más importante es la lealtad a la marca del consumidor al producto CH (LoyalCH). En función del LoyalCH, la decisión también depende del PriceDiff: Si LoyalCH es < 0,48, la predicción depende de SalePriceMM y SpecialCH. Si es mayor, depende de PriceDiff.
La interpretación el output anterior es la siguiente: por ejemplo, si escogemos el nodo terminal 7), vemos que la división viene dada por la condición LoyalCH > 0,76. Este nodo posee 283 observaciones, con un RSS de 86,5 y como clase predicha CH. Más del 96% de observaciones de entrenamiento en este nodo pertenecen a dicha clase y poco más del 3% a la clase MM.
Evaluación del modelo
# Predicción y metricas del modelo con los datos de test
predicciones <- predict(modelo_arbolC, datosOJ_test, type="class")
confusionMatrix(predicciones, datosOJ_test[ , 1])## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 108 17
## MM 22 66
##
## Accuracy : 0.8169
## 95% CI : (0.7583, 0.8664)
## No Information Rate : 0.6103
## P-Value [Acc > NIR] : 6.226e-11
##
## Kappa : 0.6192
##
## Mcnemar's Test P-Value : 0.5218
##
## Sensitivity : 0.8308
## Specificity : 0.7952
## Pos Pred Value : 0.8640
## Neg Pred Value : 0.7500
## Prevalence : 0.6103
## Detection Rate : 0.5070
## Detection Prevalence : 0.5869
## Balanced Accuracy : 0.8130
##
## 'Positive' Class : CH
## Poda
Ajuste del modelo
Para la poda de árboles de clasificación, especificamos en la función cv.tree() el argumento FUN = prune.misclass para indicar que queremos el error de clasificación como método para guiar el proceso.
# 10-fold cross validation
set.seed(356)
cv_arboles <- cv.tree(modelo_arbolC, K = 10, FUN = prune.misclass)
cv_arboles## $size
## [1] 8 4 2 1
##
## $dev
## [1] 170 170 176 334
##
## $k
## [1] -Inf 0 9 173
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"En este output, la deviance se corresponde con el error de validación.
par(mfrow = c(1, 2))
# Cost complexity pruning
plot(cv_arboles$size, cv_arboles$dev, xlab = "nodos terminales",
ylab = "error clasificación", type = "b", pch = 19)
plot(x = cv_arboles$k, y = cv_arboles$dev, xlab = "alpha",
ylab = "error clasificación", type = "b")
Podemos ver en la evolución del gráfico que con 4 nodos terminales conseguimos una reducción del error prácticamente idéntica que con 8.
Para la poda de un árbol de clasificación se utiliza la función prune.misclass(). En este caso, podaremos la mitad de los nodos (4):
# Poda del arbol a cuatro nodos terminales
modelo_arbolCpodado <- prune.misclass(tree = modelo_arbolC, best = 4)
par(mfrow = c(1,1))
plot(x = modelo_arbolCpodado, type = "proportional")
text(modelo_arbolCpodado, splits = TRUE, pretty = 0, cex = 0.8, col = "firebrick")
Evaluación del modelo
# Predicción y metricas del modelo con los datos de test
predicciones <- predict(modelo_arbolCpodado, newdata = datosOJ_test, type = "class")
confusionMatrix(predicciones, datosOJ_test[ , 1])## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 108 17
## MM 22 66
##
## Accuracy : 0.8169
## 95% CI : (0.7583, 0.8664)
## No Information Rate : 0.6103
## P-Value [Acc > NIR] : 6.226e-11
##
## Kappa : 0.6192
##
## Mcnemar's Test P-Value : 0.5218
##
## Sensitivity : 0.8308
## Specificity : 0.7952
## Pos Pred Value : 0.8640
## Neg Pred Value : 0.7500
## Prevalence : 0.6103
## Detection Rate : 0.5070
## Detection Prevalence : 0.5869
## Balanced Accuracy : 0.8130
##
## 'Positive' Class : CH
## AdaBoost
Ajuste del modelo
Para aplicar AdaBoost con caret especificamos el argumento method = "adaboost" en la función train.
# Método de validacion cruzada (10-fold)
fitControl <- trainControl(method = "cv",
number = 10,
search = "grid",
allowParallel = TRUE)
# Hiperparámetros a evaluar
grid <- expand.grid(nIter = c(100, 400, 600),
method = "adaboost") #nº de iteraciones (clasificadores)
# Evaluación del mejor conjunto de hiperparámetros
set.seed(356)
modelo_adaboost <- train(Purchase ~ ., data = datosOJ_train,
method = "adaboost",
tuneGrid = grid,
trControl = fitControl,
#metric = "Accuracy",
verbose = FALSE)
modelo_adaboost## AdaBoost Classification Trees
##
## 857 samples
## 17 predictor
## 2 classes: 'CH', 'MM'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 772, 771, 770, 771, 771, 772, ...
## Resampling results across tuning parameters:
##
## nIter Accuracy Kappa
## 100 0.7969304 0.5711528
## 400 0.7969304 0.5711528
## 600 0.7969304 0.5711528
##
## Tuning parameter 'method' was held constant at a value of adaboost
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were nIter = 100 and method = adaboost.Evaluación del modelo
# Metricas con los datos de test
confusionMatrix(predict(modelo_adaboost, datosOJ_test), datosOJ_test[ , 1])## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 106 18
## MM 24 65
##
## Accuracy : 0.8028
## 95% CI : (0.743, 0.854)
## No Information Rate : 0.6103
## P-Value [Acc > NIR] : 1.281e-09
##
## Kappa : 0.5908
##
## Mcnemar's Test P-Value : 0.4404
##
## Sensitivity : 0.8154
## Specificity : 0.7831
## Pos Pred Value : 0.8548
## Neg Pred Value : 0.7303
## Prevalence : 0.6103
## Detection Rate : 0.4977
## Detection Prevalence : 0.5822
## Balanced Accuracy : 0.7993
##
## 'Positive' Class : CH
## stopCluster(cluster)
registerDoSEQ()EJEMPLOS EN PYTHON
Para ver la aplicación de estos algoritmos sobre los mismos datos utilizando código en lenguaje Python y la libreríascikit-learn, visitar mi repositorio GitHub.
BIBLIOGRAFÍA
An Introduction to Statistical Learning: with Applications in R (2013). James G., Witten D., Hastie T., Tibshirani R.
The Elements of Statistical Learning (2009). Hastie, Trevor, Tibshirani, Robert, Friedman, Jerome
Stochastic Gradient Boosting (1999). Jerome H. Friedman
https://topepo.github.io/caret/model-training-and-tuning.html
Shield: 
This work by Cristina Gil Martínez is licensed under a Creative Commons Attribution 4.0 International License.
