Árboles

JPAG1*

1unam

*jp_6@ciencias.unam.mX

Abstract

Un árbol de decisión es una serie de filtros sobre las variables predictoras. La serie de filtros termina en una predicción de clase. Cada filtro es una pregunta binaria de sí / no, que crea bifurcaciones en la serie de filtros, lo que conduce a una estructura en forma de árbol. Los filtros dependen del tipo de variables predictoras. Si las variables son categóricas, como el género, entonces los filtros podrían ser preguntas del tipo “es género femenino.” Si las variables son continuas, como la expresión génica, el filtro podría ser “¿la expresión de PIGX es mayor que 210?” Cada punto en el que filtramos muestras en función de estas preguntas se denomina “nodos de decisión.”El algoritmo de ajuste de árbol encuentra las mejores variables en los nodos de decisión dependiendo de qué tan bien dividieron las muestras en clases después de la aplicación del nodo de decisión. Los árboles de decisión manejan variables predictoras tanto categóricas como numéricas, son fáciles de interpretar y pueden lidiar con variables faltantes.

1 Árboles de decisión

A pesar de sus ventajas, los árboles de decisión tienden a sobreajustarse si crecen a mucha profundidad y pueden aprender patrones irregulares.

Existen muchas variantes de algoritmos de aprendizaje automático basados en árboles. Sin embargo, la mayoría de los algoritmos construyen nodos de decisión de arriba hacia abajo. Seleccionan las mejores variables para usar en los nodos de decisión basándose en cuán homogéneos son los conjuntos de muestras después de la división. Una medida de homogeneidad es la impureza de Gini. Esta medida se calcula para cada subconjunto después de la división y luego se suma como un promedio ponderado.Para un nodo de decisión que divide los datos perfectamente en un problema de dos clases, la impureza de gini será \(0\) , y para un nodo que divide los datos en un subconjunto que tiene \(50%\) de clase A y \(50%\) de clase B, la impureza será \(0.5\). Formalmente, la impureza de gini, \(I_G(p)\) , de un conjunto de muestras con etiquetas de clase conocidas para K clases es la siguiente, donde \(p_i\) es la probabilidad de observar la clase \(i\) en el subconjunto:

\[\begin{equation}I_G(p) \sum_{i=1}^K p_i(1-p_i)=\sum_{i=1}^K p_i-\sum_{i=1}^K p_i^2=1-\sum_{i=1}^K p_i^2 \tag{1} \end{equation}\]

Por ejemplo, si un subconjunto de datos después de la división tiene \(75%\) de clase A y \(25%\) de clase B para ese subconjunto, la impureza sería \(1−(0.752+0.252)=0.375\). Si el otro subconjunto tuviera un \(5%\) de clase A y un \(95%\) de clase B, su impureza sería \(1−(0.952+0.052)=0.095\). Si los tamaños de los subconjuntos después de la división fueran iguales, la impureza ponderada total sería \(0.5∗0.375+0.5∗0.095=0.235\).

Estos cálculos se realizarán para cada variable potencial y la división, y cada nodo se construirá en función de la disminución de impurezas de gini. Si la variable es continua, el valor de corte se decidirá en función de la mejor impureza. Por ejemplo, los valores de expresión génica tendrán divisiones como Cell.size <2.1. Aquí (1) es el valor de corte que produce la mejor impureza. Hay otras medidas de homogeneidad, sin embargo, la impureza de gini es la que se usa para los bosques aleatorios.

Al hacer este análisis, se explora cómo el cambio de hiperparámetros puede ayudar a mejorar el rendimiento. Se utilizará parsnip (Kuhn and Vaughan 2021) para el ajuste del modelo y recipes (Kuhn and Wickham 2021) y workflows para realizar las transformaciones, y tune (Kuhn 2021) y dials (Kuhn 2020) para ajustar los hiperparámetros del modelo. rpart.plot() se usa para visualizar los árboles de decisión creados usando el paquete rpart (Therneau and Atkinson 2019) como motor, y vip (Greenwell and Boehmke 2020) se usa para visualizar la importancia de las variables para modelos posteriores.

Para ajustar un árbol de clasificación y regresión, se crean una especificación de árbol de decisión general usando rpart como motor, con la base de datos de mlbench (Leisch and Dimitriadou 2021; Newman et al. 1998), siguiendo las direcciones de tidymodels (Kuhn and Wickham 2020) y tidyverse [Wickham et al. (2019).

library(mlbench)
library(tidyverse)
library(tidymodels)
library(rpart.plot)
library(vip)
data("BreastCancer")
BCrec <- recipe(Class ~ ., data = BreastCancer) %>%
  step_naomit(all_predictors()) %>%
  prep(BreastCancer, verbose = FALSE) %>%
  bake(new_data = NULL)
BCrec2 <- BCrec %>% select(-Class,-Id) %>% mutate_all(as.numeric)
BCrec4 <- tibble(BCrec2,BCrec %>% select(Class))
tree_spec <- decision_tree() %>%
  set_engine("rpart")

Luego, esta especificación del árbol de decisión se puede utilizar para crear un motor de árbol de decisión de clasificación.

class_tree_spec <- tree_spec %>%
  set_mode("classification")

Con la especificación, se ajusta el modelo:

class_tree_fit <- class_tree_spec %>%
  fit(Class ~ ., data = BCrec4)

El resultado del modelo ofrece un resumen bastante informativo del modelo. Intenta dar una descripción escrita del árbol que se crea.

class_tree_fit$fit

## n= 683 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 683 239 benign (0.65007321 0.34992679)  
##    2) Cell.size< 2.5 418  12 benign (0.97129187 0.02870813)  
##      4) Bare.nuclei< 5.5 410   5 benign (0.98780488 0.01219512) *
##      5) Bare.nuclei>=5.5 8   1 malignant (0.12500000 0.87500000) *
##    3) Cell.size>=2.5 265  38 malignant (0.14339623 0.85660377)  
##      6) Cell.shape< 2.5 23   5 benign (0.78260870 0.21739130)  
##       12) Bl.cromatin< 3.5 16   0 benign (1.00000000 0.00000000) *
##       13) Bl.cromatin>=3.5 7   2 malignant (0.28571429 0.71428571) *
##      7) Cell.shape>=2.5 242  20 malignant (0.08264463 0.91735537)  
##       14) Cell.size< 4.5 68  17 malignant (0.25000000 0.75000000)  
##         28) Bare.nuclei< 2.5 14   4 benign (0.71428571 0.28571429) *
##         29) Bare.nuclei>=2.5 54   7 malignant (0.12962963 0.87037037) *
##       15) Cell.size>=4.5 174   3 malignant (0.01724138 0.98275862) *

Con summary() se obtiene mayor información de utilidad:

summary(class_tree_fit$fit)

## Call:
## rpart::rpart(formula = Class ~ ., data = data)
##   n= 683 
## 
##           CP nsplit  rel error    xerror       xstd
## 1 0.79079498      0 1.00000000 1.0000000 0.05215335
## 2 0.05439331      1 0.20920502 0.2677824 0.03186596
## 3 0.02510460      2 0.15481172 0.1924686 0.02740567
## 4 0.01255230      3 0.12970711 0.1882845 0.02712741
## 5 0.01000000      6 0.09205021 0.1882845 0.02712741
## 
## Variable importance
##       Cell.size      Cell.shape     Bare.nuclei    Epith.c.size     Bl.cromatin 
##              21              18              16              15              15 
## Normal.nucleoli    Cl.thickness 
##              14               1 
## 
## Node number 1: 683 observations,    complexity param=0.790795
##   predicted class=benign     expected loss=0.3499268  P(node) =1
##     class counts:   444   239
##    probabilities: 0.650 0.350 
##   left son=2 (418 obs) right son=3 (265 obs)
##   Primary splits:
##       Cell.size    < 2.5 to the left,  improve=222.3221, (0 missing)
##       Cell.shape   < 3.5 to the left,  improve=216.4111, (0 missing)
##       Bare.nuclei  < 2.5 to the left,  improve=203.7284, (0 missing)
##       Bl.cromatin  < 3.5 to the left,  improve=196.3903, (0 missing)
##       Epith.c.size < 2.5 to the left,  improve=193.1310, (0 missing)
##   Surrogate splits:
##       Cell.shape      < 3.5 to the left,  agree=0.917, adj=0.785, (0 split)
##       Epith.c.size    < 2.5 to the left,  agree=0.900, adj=0.743, (0 split)
##       Bare.nuclei     < 2.5 to the left,  agree=0.880, adj=0.691, (0 split)
##       Normal.nucleoli < 2.5 to the left,  agree=0.877, adj=0.683, (0 split)
##       Bl.cromatin     < 3.5 to the left,  agree=0.876, adj=0.679, (0 split)
## 
## Node number 2: 418 observations,    complexity param=0.0251046
##   predicted class=benign     expected loss=0.02870813  P(node) =0.6120059
##     class counts:   406    12
##    probabilities: 0.971 0.029 
##   left son=4 (410 obs) right son=5 (8 obs)
##   Primary splits:
##       Bare.nuclei     < 5.5 to the left,  improve=11.68296, (0 missing)
##       Cl.thickness    < 6.5 to the left,  improve=10.32214, (0 missing)
##       Normal.nucleoli < 3.5 to the left,  improve=10.32214, (0 missing)
##       Bl.cromatin     < 4.5 to the left,  improve= 8.53307, (0 missing)
##       Epith.c.size    < 3.5 to the left,  improve= 4.63208, (0 missing)
##   Surrogate splits:
##       Cl.thickness    < 8.5 to the left,  agree=0.988, adj=0.375, (0 split)
##       Normal.nucleoli < 3.5 to the left,  agree=0.983, adj=0.125, (0 split)
## 
## Node number 3: 265 observations,    complexity param=0.05439331
##   predicted class=malignant  expected loss=0.1433962  P(node) =0.3879941
##     class counts:    38   227
##    probabilities: 0.143 0.857 
##   left son=6 (23 obs) right son=7 (242 obs)
##   Primary splits:
##       Cell.shape    < 2.5 to the left,  improve=20.58158, (0 missing)
##       Cell.size     < 3.5 to the left,  improve=18.27650, (0 missing)
##       Bare.nuclei   < 1.5 to the left,  improve=16.81493, (0 missing)
##       Bl.cromatin   < 2.5 to the left,  improve=13.91034, (0 missing)
##       Marg.adhesion < 2.5 to the left,  improve=11.17148, (0 missing)
##   Surrogate splits:
##       Bl.cromatin < 1.5 to the left,  agree=0.932, adj=0.217, (0 split)
## 
## Node number 4: 410 observations
##   predicted class=benign     expected loss=0.01219512  P(node) =0.6002928
##     class counts:   405     5
##    probabilities: 0.988 0.012 
## 
## Node number 5: 8 observations
##   predicted class=malignant  expected loss=0.125  P(node) =0.01171303
##     class counts:     1     7
##    probabilities: 0.125 0.875 
## 
## Node number 6: 23 observations,    complexity param=0.0125523
##   predicted class=benign     expected loss=0.2173913  P(node) =0.03367496
##     class counts:    18     5
##    probabilities: 0.783 0.217 
##   left son=12 (16 obs) right son=13 (7 obs)
##   Primary splits:
##       Bl.cromatin  < 3.5 to the left,  improve=4.968944, (0 missing)
##       Cl.thickness < 4.5 to the left,  improve=3.381643, (0 missing)
##       Bare.nuclei  < 1.5 to the left,  improve=2.826087, (0 missing)
##       Mitoses      < 1.5 to the left,  improve=2.522516, (0 missing)
##       Epith.c.size < 2.5 to the left,  improve=1.992754, (0 missing)
##   Surrogate splits:
##       Cl.thickness    < 5.5 to the left,  agree=0.870, adj=0.571, (0 split)
##       Marg.adhesion   < 7   to the left,  agree=0.826, adj=0.429, (0 split)
##       Normal.nucleoli < 2.5 to the left,  agree=0.826, adj=0.429, (0 split)
##       Mitoses         < 1.5 to the left,  agree=0.826, adj=0.429, (0 split)
##       Epith.c.size    < 4   to the left,  agree=0.783, adj=0.286, (0 split)
## 
## Node number 7: 242 observations,    complexity param=0.0125523
##   predicted class=malignant  expected loss=0.08264463  P(node) =0.3543192
##     class counts:    20   222
##    probabilities: 0.083 0.917 
##   left son=14 (68 obs) right son=15 (174 obs)
##   Primary splits:
##       Cell.size     < 4.5 to the left,  improve=5.297663, (0 missing)
##       Bare.nuclei   < 2.5 to the left,  improve=4.093695, (0 missing)
##       Cell.shape    < 4.5 to the left,  improve=2.958548, (0 missing)
##       Bl.cromatin   < 3.5 to the left,  improve=2.805426, (0 missing)
##       Marg.adhesion < 5.5 to the left,  improve=2.754821, (0 missing)
##   Surrogate splits:
##       Cell.shape    < 4.5 to the left,  agree=0.789, adj=0.250, (0 split)
##       Epith.c.size  < 2.5 to the left,  agree=0.777, adj=0.206, (0 split)
##       Marg.adhesion < 1.5 to the left,  agree=0.744, adj=0.088, (0 split)
##       Bl.cromatin   < 2.5 to the left,  agree=0.736, adj=0.059, (0 split)
## 
## Node number 12: 16 observations
##   predicted class=benign     expected loss=0  P(node) =0.02342606
##     class counts:    16     0
##    probabilities: 1.000 0.000 
## 
## Node number 13: 7 observations
##   predicted class=malignant  expected loss=0.2857143  P(node) =0.0102489
##     class counts:     2     5
##    probabilities: 0.286 0.714 
## 
## Node number 14: 68 observations,    complexity param=0.0125523
##   predicted class=malignant  expected loss=0.25  P(node) =0.09956076
##     class counts:    17    51
##    probabilities: 0.250 0.750 
##   left son=28 (14 obs) right son=29 (54 obs)
##   Primary splits:
##       Bare.nuclei     < 2.5 to the left,  improve=7.600529, (0 missing)
##       Cl.thickness    < 6.5 to the left,  improve=3.558824, (0 missing)
##       Marg.adhesion   < 5.5 to the left,  improve=2.615385, (0 missing)
##       Normal.nucleoli < 2.5 to the left,  improve=1.937690, (0 missing)
##       Bl.cromatin     < 3.5 to the left,  improve=1.525641, (0 missing)
## 
## Node number 15: 174 observations
##   predicted class=malignant  expected loss=0.01724138  P(node) =0.2547584
##     class counts:     3   171
##    probabilities: 0.017 0.983 
## 
## Node number 28: 14 observations
##   predicted class=benign     expected loss=0.2857143  P(node) =0.0204978
##     class counts:    10     4
##    probabilities: 0.714 0.286 
## 
## Node number 29: 54 observations
##   predicted class=malignant  expected loss=0.1296296  P(node) =0.07906296
##     class counts:     7    47
##    probabilities: 0.130 0.870

Y se visualiza el árbol de decisiones en la Figura 1.

rpart.plot(class_tree_fit$fit)

Figure 1: Tree with most important variable to predict Class appears to be Cell.size as it forms the first node

Se puede ver ver que la variable más importante para diagnosticar parece ser la magnitud de las células, ya que forma el primer nodo.

La precisión de entrenamiento de este modelo es \(96.5%\).

augment(class_tree_fit, new_data = BCrec4) %>%
  accuracy(truth = Class, estimate = .pred_class)

## # A tibble: 1 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.968

Y su matriz de confusión es:

augment(class_tree_fit, new_data = BCrec4) %>%
  conf_mat(truth = Class, estimate = .pred_class)

##            Truth
## Prediction  benign malignant
##   benign       431         9
##   malignant     13       230

Para corregir sobreajuste, se crea una división de validación y se ajusta el modelo en el conjunto de datos de entrenamiento.

set.seed(1234)
BCrec_split <- initial_split(BCrec4)

BCrec_train <- training(BCrec_split)
BCrec_test <- testing(BCrec_split)

Se ajusta al conjunto de entrenamiento:

class_tree_fit <- fit(class_tree_spec, Class ~ ., data = BCrec_train)

Y se compara a la matriz de confusión para el conjunto de datos de entrenamiento y el conjunto de datos de prueba.

augment(class_tree_fit, new_data = BCrec_train) %>%
 conf_mat(truth = Class, estimate = .pred_class)

##            Truth
## Prediction  benign malignant
##   benign       312         7
##   malignant     10       183

El conjunto de datos de entrenamiento funciona bien como cabría esperar

augment(class_tree_fit, new_data = BCrec_test) %>%
 conf_mat(truth = Class, estimate = .pred_class)

##            Truth
## Prediction  benign malignant
##   benign       119         4
##   malignant      3        45

pero el conjunto de datos de prueba no funciona tan bien y tiene una precisión menor del \(95.9%\).

augment(class_tree_fit, new_data = BCrec_test) %>%
 accuracy(truth = Class, estimate = .pred_class)

## # A tibble: 1 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.959

Para encontrar una complejidad más óptima del árbol de decisión. Se usa el objeto class_tree_spec y se usa la función set_args() para especificar ajuste de cost_complexity. A continuación, se pasa directamente al objeto workflow para evitar la creación de un objeto intermedio.

class_tree_wf <- workflow() %>%
  add_model(class_tree_spec %>% set_args(cost_complexity = tune())) %>%
  add_formula(Class ~ .)

Para poder ajustar la variable se vuelve a muestrear el objeto, usando validación cruzada y una cuadrícula de valores para probar. Se hace con una cuadrícula regular, dado que solo se ajusta \(1\) hiperparámetro.

set.seed(1234)
BCrec_fold <- vfold_cv(BCrec_train)
 
param_grid <- grid_regular(cost_complexity(range = c(-3, -1)), levels = 10)

tune_res <- tune_grid(
  class_tree_wf, 
  resamples = BCrec_fold, 
  grid = param_grid, 
  metrics = metric_set(accuracy)
  )

El uso de autoplot() en la Figura 2 muestra qué valores de cost_complexity parecen producir la mayor precisión.

autoplot(tune_res)

Figure 2: Cost-complexity vs accuracy

Ahora se selecciona el valor de mejor rendimiento con select_best(), se finaliza el workflow actualizando el valor de cost_complexity y ajustando el modelo en el conjunto de datos de entrenamiento completo.

best_complexity <- select_best(tune_res)
 
class_tree_final <- finalize_workflow(class_tree_wf, best_complexity)
 
class_tree_final_fit <- fit(class_tree_final, data = BCrec_train)
class_tree_final_fit

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: decision_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Class ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## n= 512 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 512 190 benign (0.62890625 0.37109375)  
##    2) Cell.size< 2.5 300  10 benign (0.96666667 0.03333333)  
##      4) Bare.nuclei< 5.5 293   4 benign (0.98634812 0.01365188) *
##      5) Bare.nuclei>=5.5 7   1 malignant (0.14285714 0.85714286) *
##    3) Cell.size>=2.5 212  32 malignant (0.15094340 0.84905660)  
##      6) Cell.shape< 2.5 18   3 benign (0.83333333 0.16666667) *
##      7) Cell.shape>=2.5 194  17 malignant (0.08762887 0.91237113)  
##       14) Bare.nuclei< 2.5 25   9 malignant (0.36000000 0.64000000)  
##         28) Cell.size< 3.5 8   0 benign (1.00000000 0.00000000) *
##         29) Cell.size>=3.5 17   1 malignant (0.05882353 0.94117647) *
##       15) Bare.nuclei>=2.5 169   8 malignant (0.04733728 0.95266272) *

Y el modelo se visualiza en la Figura 3 y se ve que el modelo de mejor rendimiento es menos complejo que el modelo original ajustado.

rpart.plot(class_tree_final_fit$fit$fit$fit)

Figure 3: Better-performing model is less complex

2 Árboles de regresión

Para ajustar un árbol de regresión, la principal diferencia aquí es que la respuesta será continua en lugar de categórica. Se reutiliza tree_spec como base para la especificación del árbol de decisión de regresión.

reg_tree_spec <- tree_spec %>%
  set_mode("regression")

Se hace una división de validación:

set.seed(1234)
BCrec_split <- initial_split(BCrec4)
 
BCrec_train <- training(BCrec_split)
BCrec_test <- testing(BCrec_split)

Y se ajusta el modelo al conjunto de datos de entrenamiento:

reg_tree_fit <- fit(reg_tree_spec, Cell.size ~ ., BCrec_train)
reg_tree_fit

## parsnip model object
## 
## Fit time:  18ms 
## n= 512 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 512 4735.21900 3.210938  
##    2) Cell.shape< 3.5 333  412.07810 1.402402  
##      4) Class=benign 307   95.47883 1.201954 *
##      5) Class=malignant 26  158.61540 3.769231  
##       10) Marg.adhesion< 6.5 19   44.52632 2.842105 *
##       11) Marg.adhesion>=6.5 7   53.42857 6.285714 *
##    3) Cell.shape>=3.5 179 1207.73200 6.575419  
##      6) Cell.shape< 6.5 88  358.07950 4.897727  
##       12) Class=benign 12   18.91667 2.916667 *
##       13) Class=malignant 76  284.63160 5.210526 *
##      7) Cell.shape>=6.5 91  362.43960 8.197802  
##       14) Cell.shape< 9.5 51  175.33330 7.333333 *
##       15) Cell.shape>=9.5 40  100.40000 9.300000 *

augment(reg_tree_fit, new_data = BCrec_test) %>%
  rmse(truth = Cell.size, estimate = .pred)

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        1.11

Se pasa la función rpart.plot() en la Figura 4.

rpart.plot(reg_tree_fit$fit)

Figure 4: Cell.shape shows better prediction for Cell.size

El resultado es una variable numérica en lugar de una clase.

Ahora se ajusta de nuevo cost_complexity para encontrar el modelo de mejor rendimiento.

reg_tree_wf <- workflow() %>%
  add_model(reg_tree_spec %>% set_args(cost_complexity = tune())) %>%
  add_formula(Cell.size ~ .)

set.seed(1234)
BCrec_fold <- vfold_cv(BCrec_train)

param_grid <- grid_regular(cost_complexity(range = c(-4, -1)), levels = 10)

tune_res <- tune_grid(
  reg_tree_wf, 
  resamples = BCrec_fold, 
  grid = param_grid
)

Y parece que una mayor complejidad es preferibles de acuerdo con la validación cruzada (Figura 5).

autoplot(tune_res)

Figure 5: Cost-complexity vs mse and rsq

Se selecciona el modelo de mejor rendimiento de acuerdo con rmse y se ajusta el modelo final a todo el conjunto de datos de entrenamiento.

best_complexity <- select_best(tune_res, metric = "rmse")

reg_tree_final <- finalize_workflow(reg_tree_wf, best_complexity)

reg_tree_final_fit <- fit(reg_tree_final, data = BCrec_train)
reg_tree_final_fit

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: decision_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Cell.size ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## n= 512 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##   1) root 512 4735.219000 3.210938  
##     2) Cell.shape< 3.5 333  412.078100 1.402402  
##       4) Class=benign 307   95.478830 1.201954  
##         8) Cell.shape< 1.5 247   31.538460 1.076923 *
##         9) Cell.shape>=1.5 60   44.183330 1.716667  
##          18) Epith.c.size< 2.5 50   32.180000 1.580000  
##            36) Bare.nuclei< 1.5 41   17.951220 1.414634  
##              72) Marg.adhesion< 1.5 33   11.333330 1.333333 *
##              73) Marg.adhesion>=1.5 8    5.500000 1.750000 *
##            37) Bare.nuclei>=1.5 9    8.000000 2.333333 *
##          19) Epith.c.size>=2.5 10    6.400000 2.400000 *
##       5) Class=malignant 26  158.615400 3.769231  
##        10) Marg.adhesion< 6.5 19   44.526320 2.842105 *
##        11) Marg.adhesion>=6.5 7   53.428570 6.285714 *
##     3) Cell.shape>=3.5 179 1207.732000 6.575419  
##       6) Cell.shape< 6.5 88  358.079500 4.897727  
##        12) Class=benign 12   18.916670 2.916667 *
##        13) Class=malignant 76  284.631600 5.210526  
##          26) Bl.cromatin< 7.5 66  214.984800 4.984848  
##            52) Mitoses< 2.5 48  111.666700 4.583333  
##             104) Cell.shape< 5.5 33   44.909090 4.181818  
##               208) Bare.nuclei>=3.5 26   22.653850 3.884615  
##                 416) Marg.adhesion< 2.5 9    3.555556 3.222222 *
##                 417) Marg.adhesion>=2.5 17   13.058820 4.235294 *
##               209) Bare.nuclei< 3.5 7   11.428570 5.285714 *
##             105) Cell.shape>=5.5 15   49.733330 5.466667 *
##            53) Mitoses>=2.5 18   74.944440 6.055556 *
##          27) Bl.cromatin>=7.5 10   44.100000 6.700000 *
##       7) Cell.shape>=6.5 91  362.439600 8.197802  
##        14) Cell.shape< 9.5 51  175.333300 7.333333  
##          28) Bare.nuclei>=5.5 35  120.285700 6.857143  
##            56) Marg.adhesion< 3.5 7   19.428570 5.285714 *
##            57) Marg.adhesion>=3.5 28   79.250000 7.250000  
##             114) Normal.nucleoli>=8.5 8   11.500000 6.250000 *
##             115) Normal.nucleoli< 8.5 20   56.550000 7.650000  
##               230) Cell.shape< 7.5 11   32.000000 7.000000 *
##               231) Cell.shape>=7.5 9   14.222220 8.444444 *
##          29) Bare.nuclei< 5.5 16   29.750000 8.375000 *
##        15) Cell.shape>=9.5 40  100.400000 9.300000  
##          30) Epith.c.size< 5.5 13   59.230770 8.461538 *
##          31) Epith.c.size>=5.5 27   27.629630 9.703704  
##            62) Bl.cromatin< 6.5 7   22.000000 9.000000 *
##            63) Bl.cromatin>=6.5 20    0.950000 9.950000 *

La visualización del modelo (Figura 6) revela un árbol mucho más complejo.

rpart.plot(reg_tree_final_fit$fit$fit$fit)

Figure 6: Regression tree

3 Bagging y random forests

Los bosques aleatorios están diseñados para contrarrestar las deficiencias de los árboles de decisión. Son simplemente conjuntos de árboles de decisión. Cada árbol se entrena con una parte diferente seleccionada al azar de los datos con variables predictoras seleccionadas al azar. El objetivo de introducir la aleatoriedad es reducir la varianza del modelo para que no se sobreajuste, a expensas de un pequeño aumento en el sesgo y cierta pérdida de interpretabilidad. Esta estrategia generalmente mejora el rendimiento del modelo final.

El algoritmo de bosques aleatorios intenta descorrelacionar los árboles para que aprendan cosas diferentes sobre los datos. Lo hace seleccionando un subconjunto aleatorio de variables. Si una o unas pocas variables predictoras son predictores muy fuertes para la variable de respuesta, estas características se seleccionarán en muchos de los árboles, haciendo que se correlacionen. El submuestreo aleatorio de las variables predictoras garantiza que no siempre se seleccionen los mejores predictores en general para cada árbol y que el modelo tenga la oportunidad de aprender otras características de los datos.

Otro método de muestreo introducido al construir bosques aleatorios es el remuestreo previo antes de construir cada árbol. Esto trae la ventaja de la predicción de errores out-of-the-bag (OOB). En este caso, el error de predicción se puede estimar para las muestras de entrenamiento que fueron OOB, lo que significa que no se utilizaron en el entrenamiento, para algún porcentaje de los árboles. El error de predicción para cada muestra se puede estimar a partir de los árboles donde esa muestra fue OOB. Las estimaciones OOB afirmaron ser una buena alternativa a los errores estimados de validación cruzada.

Aquí se usará el paquete randomForest (Liaw and Wiener 2002) como motor. Un modelo de bagging o ensacado es lo mismo que un bosque aleatorio donde mtry es igual al número de predictores, especificando mtry=.cols(), lo que significa que se usa el número de columnas en la matriz del predictor. Esto es útil si desea que la especificación sea más general y utilizable para muchos conjuntos de datos diferentes. .cols () es uno de los muchos descriptores del paquete parsnip. También se establece important = TRUE en set_engine() para decirle al motor que guarde la información relacionada con la importancia de la variable. Esto es necesario para usar el paquete vip.

bagging_spec <- rand_forest(mtry = .cols()) %>%
  set_engine("randomForest", importance = TRUE) %>%
  set_mode("regression")

Se ajusta el modelo

bagging_fit <- fit(bagging_spec, Cell.size ~ ., data = BCrec_train)

Y se visualiza el rendimiento de las pruebas, buscando una mejora con respecto al árbol de decisiones.

augment(bagging_fit, new_data = BCrec_test) %>%
  rmse(truth = Cell.size, estimate = .pred)

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        1.01

En la Figura 7, una gráfica de dispersión rápida entre el valor verdadero y el predicho del diagnóstico.

augment(bagging_fit, new_data = BCrec_test) %>%
  ggplot(aes(Cell.size, .pred)) +
  geom_abline() +
  geom_point(alpha = 0.5)

Figure 7: Cell.size vs pred fit

3.1 Importancia de las variables

Los bosques aleatorios vienen con métricas de importancia variable integradas. Una de las métricas es similar a la “métrica de abandono de variables,” donde las variables predictoras se permutan. En este caso, se utilizan muestras OOB y las variables se permutan una a la vez. Cada vez, las muestras con las variables permutadas se alimentan a la red y se mide la disminución de la precisión. Usando esta cantidad, las variables se pueden clasificar.

Un método menos costoso con un rendimiento similar es utilizar la impureza de gini. Cada vez que se usa una variable en un árbol para hacer una división, la impureza de gini es menor que el nodo madre. Este método suma estas disminuciones de impurezas de gini para cada variable individual en los árboles y las divide por el número de árboles en el bosque. Esta métrica a menudo es consistente con la medida de importancia de la permutación (Breiman 2001).

Para acceder a los valores de importancia y visualizarlos, el resultado se muestra en la Figura 8.

vip(bagging_fit)

Figure 8: Variable importance

3.2 Random forests

Por defecto, se consideran \(p / 3\) variables cuando se construye un bosque aleatorio de árboles de regresión, y seconsideran \(\sqrt(p)\) variables cuando se construye un bosque aleatorio de árboles de clasificación. Aquí, se usa mtry = 3.

rf_spec <- rand_forest(mtry = 3) %>%
set_engine("randomForest", importance = TRUE) %>%
set_mode("regression")

Se ajusta el modelo como de costumbre

rf_fit <- fit(rf_spec, Cell.size ~ ., data = BCrec_train)

este modelo tiene un rendimiento ligeramente mejor que el modelo de bagging

augment(rf_fit, new_data = BCrec_test) %>%
  rmse(truth = Cell.size, estimate = .pred)

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        1.01

Del mismo modo, la grafica del valor real contra el valor predicho se muestra en la Figura 9.

augment(rf_fit, new_data = BCrec_test) %>%
  ggplot(aes(Cell.size, .pred)) +
  geom_abline() +
  geom_point(alpha = 0.5)

Figure 9: Cell.size vs pred
Random forest model

Se ve bien. No hay diferencia apreciable entre esta gráfica y la del modelo de bagging.

La gráfica de importancia de variables también es bastante similar al bagging (Figura 10).

vip(rf_fit)

Figure 10: Variable importance
Random forest model

4 Boosting

La paquetería xgboost (Chen et al. 2021) ofrece una buena implementación de árboles potenciados. Se establece tree = 5000 para que crezcan \(5000\) árboles con una profundidad máxima de \(4\) estableciendo tree_depth = 4.

library(xgboost)
boost_spec <- boost_tree(trees = 5000, tree_depth = 4) %>%
  set_engine("xgboost") %>%
  set_mode("regression")

Ajuste:

boost_fit <- fit(boost_spec, Cell.size ~ ., data = BCrec_train)

y se muestra el rmse:

augment(boost_fit, new_data = BCrec_test) %>%
rmse(truth = Cell.size, estimate = .pred)

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        1.33

Y se observa el diagrama de dispersión en la Figura 11.

augment(boost_fit, new_data = BCrec_test) %>%
    ggplot(aes(Cell.size, .pred)) +
    geom_abline() +
    geom_point(alpha = 0.5)

Figure 11: Cell.size vs pred
Boosting

Referencias

Breiman, Leo. 2001. “Random Forests.” Machine Learning 45 (1): 5–32.

Chen, Tianqi, Tong He, Michael Benesty, Vadim Khotilovich, Yuan Tang, Hyunsu Cho, Kailong Chen, et al. 2021. Xgboost: Extreme Gradient Boosting. https://CRAN.R-project.org/package=xgboost.

Greenwell, Brandon M., and Bradley C. Boehmke. 2020. “Variable Importance Plots—an Introduction to the Vip Package.” The R Journal 12 (1): 343–66. https://doi.org/10.32614/RJ-2020-013.

Kuhn, Max. 2020. Dials: Tools for Creating Tuning Parameter Values. https://CRAN.R-project.org/package=dials.

———. 2021. Tune: Tidy Tuning Tools. https://CRAN.R-project.org/package=tune.

Kuhn, Max, and Davis Vaughan. 2021. Parsnip: A Common API to Modeling and Analysis Functions. https://CRAN.R-project.org/package=parsnip.

Kuhn, Max, and Hadley Wickham. 2020. Tidymodels: A Collection of Packages for Modeling and Machine Learning Using Tidyverse Principles. https://www.tidymodels.org.

———. 2021. Recipes: Preprocessing Tools to Create Design Matrices. https://CRAN.R-project.org/package=recipes.

Leisch, Friedrich, and Evgenia Dimitriadou. 2021. Mlbench: Machine Learning Benchmark Problems.

Liaw, Andy, and Matthew Wiener. 2002. “Classification and Regression by randomForest.” R News 2 (3): 18–22. https://CRAN.R-project.org/doc/Rnews/.

Newman, D. J., S. Hettich, C. L. Blake, and C. J. Merz. 1998. “UCI Repository of Machine Learning Databases.” University of California, Irvine, Dept. of Information; Computer Sciences. http://www.ics.uci.edu/~mlearn/MLRepository.html.

Therneau, Terry, and Beth Atkinson. 2019. Rpart: Recursive Partitioning and Regression Trees. https://CRAN.R-project.org/package=rpart.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.