XGBoost o Extreme Gradient Boosting, es uno de los algoritmos de machine learning de tipo supervisado más usados en la actualidad.

Este algoritmo se caracteriza por obtener buenos resultados de predicción con relativamente poco esfuerzo, en muchos casos equiparables o mejores que los devueltos por modelos más complejos computacionalmente, en particular para problemas con datos heterogéneos.

XGBoost es una herramienta muy útil para un data scientist y cuenta con implementaciones para diferentes lenguajes y entornos de programación.

En este artículo revisaremos la implementación de XGBoost en R. Veremos cómo preparar los datos para usar este algoritmo, sus hiper parámetros básicos y una manera sencilla de evaluar sus resultados.

Pero antes, una breve introducción a XGBoost.

Una introducción informal a XGBoost

XGBoost Extreme Gradient Boosting es un algoritmo predictivo supervisado que utiliza el principio de boosting.

La idea detrás del boosting es generar múltiples modelos de predicción “débiles” secuenciualmente,y que cada uno de estos tome los resultados del modelo anterior, para generar un modelo más “fuerte”, con mejor poder predictivo y mayor estabilidad en sus resultados.

Para conseguir un modelo más fuerte, se emplea un algoritmo de optimización, este caso Gradient Descent (descenso de gradiente).

Durante el entrenamiento, los parámetros de cada modelo débil son ajustados iterativamente tratando de encontrar el mínimo de una función objetivo, que puede ser la proporción de error en la clasificación, el área bajo la curva (AUC), la raíz del error cuadrático medio (RMSE) o alguna otra.

Cada modelo es comparado con el anterior. Si un nuevo modelo tiene mejores resultados, entonces se toma este como base para realizar nuevas modificaciones. Si, por el contrario, tiene peores resultados, se regresa al mejor modelo anterior y se modifica ese de una manera diferente.

Este proceso se repite hasta llegar a un punto en el que la diferencia entre modelos consecutivos es insignificante, lo cual nos indica que hemos encontrado el mejor modelo posible, o cuando se llega al número de iteraciones máximas definido por el usuario.

XGBoost usa como sus modelos débiles árboles de decisión de diferentes tipos, que pueden ser usados para tareas de clasificación y de regresión.

Si quieres conocer más sobre este algoritmo, puedes leer definiciones más formales que incluyen discusión sobre su implementación en los siguientes artículos:

Ahora veamos cómo usar XGBoost en R.

Implementación en R

Instalación

Lo primero que necesitamos es instalar los paquetes que usaremos con install.packages().

xgboost contiene la implementación de este tipo de modelo para R.

tidyverse es un meta paquete que carga varios paquetes a nuestro entorno de trabajo. Para nuestros fines, lo importante es que carga los siguientes:

  • dplyr. Herramientas manipular, transformar y organizar datos.
  • readr. Facilita importar datos desde archivos.
  • purrr. Agrega características de programación funcional.

caret contiene una gran cantidad de utilidades para machine learning, pero en esta ocasión únicamente usaremos una función para generar matrices de confusión.

Cargamos los paquetes a nuestro espacio de trabajo con library().

-- Attaching packages ----------------------------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.2.1     v purrr   0.3.3
v tibble  2.1.3     v dplyr   0.8.3
v tidyr   1.0.0     v stringr 1.4.0
v readr   1.3.1     v forcats 0.4.0
-- Conflicts -------------------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
x dplyr::slice()  masks xgboost::slice()
Loading required package: lattice

Attaching package: 'caret'
The following object is masked from 'package:purrr':

    lift

Datos que usaremos

El conjunto de datos que usaremos es conocido como Agraricus. Contiene características de diferentes hongos y lo que deseamos predecir es si son venenosos o no.

En la práctica es común que tengas que lidiar con datos en formatos no convencionales que requieren preparación antes de ser usables por un algoritmo de machine learning.

Por esta razón, en lugar de usar la versión de estos datos incluida en el paquete xgboost y que ya está lista para usar, trabajaremos con una versión de estos mismos datos que requiere preparación, pues generalmente es la etapa del flujo de trabajo de machine learning que requiere más tiempo.

La versión de los datos que usaremos está disponible en el Machine Learning Repository de UCI.

He copiado los datos a un repositorio de github para asegurar que estés usando la misma versión que aparece en este artículo. Son dos archivos en total, uno con extensión .data que contiene los datos de los hongos, y otro de extensión .names que incluye una descripción de ellos.

Descargamos ambos archivos a nuestra carpeta de trabajo usando la download.files() con el argumento mode = "wb" para asegurar que los archivos se guarden correctamente.

Exploración de los datos

Empezamos con la exploración del contenido de estos archivos.

Podemos dar una mirada a su contenido con algún procesador de texto externo, como notepad++ o gedit, pero también podemos hacer esto directamente en R con las funciones read_lines() y head().

Veamos los primeros renglones del archivo “agaricus-lepiota.data”.

[1] "p,x,s,n,t,p,f,c,n,k,e,e,s,s,w,w,p,w,o,p,k,s,u"
[2] "e,x,s,y,t,a,f,c,b,k,e,c,s,s,w,w,p,w,o,p,n,n,g"
[3] "e,b,s,w,t,l,f,c,b,n,e,c,s,s,w,w,p,w,o,p,n,n,m"
[4] "p,x,y,w,t,p,f,c,n,n,e,e,s,s,w,w,p,w,o,p,k,s,u"
[5] "e,x,s,g,f,n,f,w,b,k,t,e,s,s,w,w,p,w,o,e,n,a,g"
[6] "e,x,y,y,t,a,f,c,b,n,e,c,s,s,w,w,p,w,o,p,k,n,g"

Los datos se encuentran en una estructura tabular, con columnas separada por comas. Para fines prácticos, es equivalente a un archivo de extensión .csv pero con una extensión diferente. Eso son buenas noticias.

Para evitar errores en la lectura, importaremos su contenido con la función read.table Esta función no trata de convertir los datos a un tipo en particular, de modo que todo será importado como de tipo caracter, lo cual previene errores más adelante.

Llamamos a esta función especificando que el delimitador de columnas es una coma y convertimos su resultado a un tibble para mejorar la compatibilidad con los paquetes del tidyverse y mejorar su presentación en pantalla.

Veamos el resultado.

# A tibble: 8,124 x 23
   V1    V2    V3    V4    V5    V6    V7    V8    V9    V10   V11   V12   V13  
   <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct>
 1 p     x     s     n     t     p     f     c     n     k     e     e     s    
 2 e     x     s     y     t     a     f     c     b     k     e     c     s    
 3 e     b     s     w     t     l     f     c     b     n     e     c     s    
 4 p     x     y     w     t     p     f     c     n     n     e     e     s    
 5 e     x     s     g     f     n     f     w     b     k     t     e     s    
 6 e     x     y     y     t     a     f     c     b     n     e     c     s    
 7 e     b     s     w     t     a     f     c     b     g     e     c     s    
 8 e     b     y     w     t     l     f     c     b     n     e     c     s    
 9 p     x     y     w     t     p     f     c     n     p     e     e     s    
10 e     b     s     y     t     a     f     c     b     g     e     c     s    
# ... with 8,114 more rows, and 10 more variables: V14 <fct>, V15 <fct>,
#   V16 <fct>, V17 <fct>, V18 <fct>, V19 <fct>, V20 <fct>, V21 <fct>,
#   V22 <fct>, V23 <fct>

Continuamos con la preparación de los datos.

Preparación de los datos

Asignación de nombres a las variables

Aunque este paso no es estrictamente necesario, agregaremos el nombre de cada columna, es decir, de las variables o features.

Es frecuente que en la práctica tengas que trabajar con datos a los que se les ha ocultado intencionalmente el nombre de los features por seguridad o confidencialidad, entre otras razones. Sin embargo, cuando cuenta con los nombres de las variables pueden encontrar insights sobre nuestros datos que pueden ser útiles para realizar un buen análisis.

Además, en nuestro ejemplo es esencial conocer cuál de las columnas es la variable objetivo, aquella en la que se encuentra etiquetado un hongo como venenoso o no.

Los nombres de las variables han sido obtenidos de la información en el archivo agaricus-lepiota.names. La variable objetivo es la primera columna con el nombre target.

Haremos la asignación de nombres con la función names().

Nuestro resultado es el siguiente.

# A tibble: 6 x 23
  target cap_shape cap_surface cap_color bruises odor  gill_attachment
  <fct>  <fct>     <fct>       <fct>     <fct>   <fct> <fct>          
1 p      x         s           n         t       p     f              
2 e      x         s           y         t       a     f              
3 e      b         s           w         t       l     f              
4 p      x         y           w         t       p     f              
5 e      x         s           g         f       n     f              
6 e      x         y           y         t       a     f              
# ... with 16 more variables: gill_spacing <fct>, gill_size <fct>,
#   gill_color <fct>, stalk_shape <fct>, stalk_root <fct>,
#   stalk_surface_above_ring <fct>, stalk_surface_below_ring <fct>,
#   stalk_color_above_ring <fct>, stalk_color_below_ring <fct>,
#   veil_type <fct>, veil_color <fct>, ring_number <fct>, ring_type <fct>,
#   spore_print_color <fct>, population <fct>, habitat <fct>

Conversión a numérico

Hasta aquí, todo luce bien, sin embargo, xgboost requiere matrices numéricas para funcionar correctamente. Tenemos que convertir nuestra columnas de datos de tipo caracter a tipo numérico.

Utilizaremos las funciones map_df() de purrr para realizar esta conversión de manera eficiente.

Cada columna de nuestros datos será convertida a factor con as.factor(), y después a número con as.numeric(). En R si intentamos convertir un dato de tipo caracter a numérico obtenemos únicamente NAs, por eso necesitamos un paso intermedio.

Finalmente, restamos 1 al resultado de la conversión porque XGBoost espera valores de 0 y 1 para la variable objetivo. En nuestro ejemplo, el 1 representa un hongo venenoso y 0 un hongo que no lo es.

Si no realizamos esta conversión, es probable que XGBoost devuelva resultados incorrectos.

Nuestro resultado es el siguiente.

# A tibble: 6 x 23
  target cap_shape cap_surface cap_color bruises  odor gill_attachment
   <dbl>     <dbl>       <dbl>     <dbl>   <dbl> <dbl>           <dbl>
1      1         5           2         4       1     6               1
2      0         5           2         9       1     0               1
3      0         0           2         8       1     3               1
4      1         5           3         8       1     6               1
5      0         5           2         3       0     5               1
6      0         5           3         9       1     0               1
# ... with 16 more variables: gill_spacing <dbl>, gill_size <dbl>,
#   gill_color <dbl>, stalk_shape <dbl>, stalk_root <dbl>,
#   stalk_surface_above_ring <dbl>, stalk_surface_below_ring <dbl>,
#   stalk_color_above_ring <dbl>, stalk_color_below_ring <dbl>,
#   veil_type <dbl>, veil_color <dbl>, ring_number <dbl>, ring_type <dbl>,
#   spore_print_color <dbl>, population <dbl>, habitat <dbl>

Creación de una lista

Este paso tampoco es estrictamente necesario, pero a mí me resulta más fácil guardar todos los objetos relacionados con un mismo proceso en una lista.

De este modo, mi espacio de trabajo es más fácil organizar. Mientras sea consistente con la estructura de esta lista es posible reproducir o hacer ajustes a los análisis e identificar en qué paso del proceso se han presentado problemas.

Creamos una lista llamada hongo y lo primero que haremos es guardar nuestros datos originales en ella.

Sets de entrenamiento y prueba

Como es el caso para todos los algoritmos de predicción supervisados, necesitamos dividir nuestros datos en un conjunto de entrenamiento, que aprenderá las características de los datos y generará un modelo de predicción; y un conjunto de prueba, con el que validamos el modelo generado.

Generamos nuestro set de entrenamiento con la función sample_frac() de dplyr. Extraemos una muestra aleatoria del 70% de nuestros datos originales. Utilizamos set.seed() para asegurar que estos resultados son replicables.

El resto de los datos, 30% de ellos, será el conjunto de prueba (test).

Usamos setdiff() para seleccionarlos.

En ambos casos, hemos guardado los resultados en nuestra lista hongo.

Veamos el tamaño de estos conjuntos de datos con la función dim.

[1] 5687   23
[1] 2437   23

Convertir a DMatrix

Como ya lo mencionamos, la implementación XGBoost de R requiere que los datos que usemos sean matrices, específicamente de tipo DMatrix, así que necesitamos convertir nuestros sets de entrenamiento y prueba a este tipo de estructura.

Usaremos la función xgb.DMatrix() de xgboost para la conversión.

Esta función espera una matriz numérica como primer argumento y también se pueden especificar algunos atributos adicionales al objeto que devolverá. Nosotros definiremos el atributo label para identificar la variable objetivo en nuestros datos.

Al usar esta función es muy importante que tu data no incluya la columna con la variable objetivo, de lo contrario, obtendrás una precisión perfecta en tus predicciones, que será inútil con datos nuevos.

Entonces, quitamos la variable objetivo de los datos usando la función select() de dplyr, convertiremos nuestros datos a matriz con as.matrix(), convertimos esta matriz con xgb.Dmatrix() y guardamos el resultado en nuestra lista hongo.

Nuestro resultado es el siguiente.

xgb.DMatrix  dim: 5687 x 22  info: label  colnames: yes

Realizamos el mismo procedimiento con nuestro set de prueba.

Obtenemos el siguiente resultado.

xgb.DMatrix  dim: 2437 x 22  info: label  colnames: yes

¡Listo! Hemos concluido la parte más laboriosa de todo el proceso. Podemos comenzar con el entrenamiento del modelo predictivo.

Entrenamiento del modelo predictivo

Para entrenar un modelo usamos la función xgboost().

Tenemos a nuestra disposición una amplia cantidad de hiper parámetros para ajustar, pero para este ejemplo introductorio haremos ajustes solo a los siguientes:

Como los datos de nuestro ejemplo son sencillos, definiremos valores muy conservadores para todos estos hiper parámetros.

Entrenamos el modelo y lo guardamos nuestro modelo en hongo$modelo_01. Para datos más complejos, este proceso puede ser tardado.

Por defecto, nos es mostrado el resultado de la función objetivo de cada iteración. De este modo podemos analizar el desempeño del modelo que hemos especificado.

[1] train-error:0.089678 
[2] train-error:0.073853 
[3] train-error:0.045894 
[4] train-error:0.036047 
[5] train-error:0.029189 
[6] train-error:0.026200 
[7] train-error:0.024618 
[8] train-error:0.024618 
[9] train-error:0.015298 
[10]    train-error:0.015298 

Nuestro resultado es el siguiente. Notarás que la salida incluye información de los hiper parámetros y los datos que usamos, así como un resumen de la ejecución del modelo.

##### xgb.Booster
raw: 4.2 Kb 
call:
  xgb.train(params = params, data = dtrain, nrounds = nrounds, 
    watchlist = watchlist, verbose = verbose, print_every_n = print_every_n, 
    early_stopping_rounds = early_stopping_rounds, maximize = maximize, 
    save_period = save_period, save_name = save_name, xgb_model = xgb_model, 
    callbacks = callbacks, objective = "binary:logistic", max.depth = 2, 
    eta = 0.3, nthread = 2)
params (as set within xgb.train):
  objective = "binary:logistic", max_depth = "2", eta = "0.3", nthread = "2", silent = "1"
xgb.attributes:
  niter
callbacks:
  cb.print.evaluation(period = print_every_n)
  cb.evaluation.log()
# of features: 22 
niter: 10
nfeatures : 22 
evaluation_log:
    iter train_error
       1    0.089678
       2    0.073853
---                 
       9    0.015298
      10    0.015298

Generación de predicciones

El siguiente paso es utilizar la función predict() con el set de prueba hongo$test_mat para generar las predicciones de nuestro modelo hongo$modelo_01.

Esta función espera un modelo de predicción y datos nuevos con la misma estructura que los usados para entrenar al modelo. Esto es muy importante: si tus datos de prueba tienen una estructura diferente a los de entrenamiento, no podrás obtener predicciones.

Guardamos las predicciones en la lista hongo.

Nuestro resultado es un vector de valores numéricos, cada uno representando la probabilidad de que un caso en particular pertenezca al valor 1 de nuestra variable objetivo. Es decir, la probabilidad de que ese hongo sea venenoso.

[1] 0.09554099 0.11355446 0.04708923 0.86634964 0.04708923 0.86634964

Para este ejemplo, tomaremos las probabilidades mayores a 0.5 como una predicción de pertenencia al valor 1 de nuestra variable objetivo.

[1] FALSE FALSE FALSE  TRUE FALSE  TRUE

Evaluación del modelo

Para evaluar nuestro modelo comparamos nuestras predicciones con las categorías reales de nuestro set de prueba.

Para esto, usaremos la función confusionMatrix() de caret para generar un matriz de confusión.

La manera más sencilla de utilizar esta función es darle como argumento un objeto de tipo table.

Unimos nuestras predicciones en hongo$predict_01 y los valores reales en hongo$test_df$target con cbind() para generar un data frame con data.frame(), y con ella un objeto table con table().

Confusion Matrix and Statistics

   X2
X1     0    1
  0 1246   45
  1    0 1146
                                          
               Accuracy : 0.9815          
                 95% CI : (0.9754, 0.9865)
    No Information Rate : 0.5113          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.963           
                                          
 Mcnemar's Test P-Value : 5.412e-11       
                                          
            Sensitivity : 1.0000          
            Specificity : 0.9622          
         Pos Pred Value : 0.9651          
         Neg Pred Value : 1.0000          
             Prevalence : 0.5113          
         Detection Rate : 0.5113          
   Detection Prevalence : 0.5297          
      Balanced Accuracy : 0.9811          
                                          
       'Positive' Class : 0               
                                          

Nada mal. Tuvimos una precisión del 98% y sobresale que tuvimos un valor predictivo par la clase negativa, un hongo no venenoso, del 100%.

Como podrás imaginar, con datos reales rara vez obtenemos resultados tan buenos con tan poco esfuerzo, pero si comparas estos resultados contra los de árboles de decisión convencionales, notarás una gran diferencia en desempeño a favor de XGBoost.

Después de preparar nuestros datos, la tarea que más tiempo consume al usar este modelo es encontrar los mejores hiper parámetros para alcanzar la mayor precisión posible de un modelo.

Veamos qué pasa si ajustamos nuestros hiper parámetros con un segundo modelo.

Segundo modelo.

Este segundo modelo tiene un número de iteraciones mayor que el anterior, de 100 en lugar de 10, y una mayor profundidad en los árboles generados, de 2 a 4.

Además, hemos ajustado el hiper parámetro early_stopping_rounds = 10, para que el entrenamiento se detenga si después de diez iteraciones consecutivas no se mejora el modelo. Este hiper parámetro es sumamente útil para acortar el tiempo de entrenamiento de un modelo, pues evita que el proceso continúe si ya no se obtienen mejores resultados de predicción.

En este mismo bloque de código generamos una matriz de confusión para evaluar nuestro segundo modelo.

[1] train-error:0.020573 
Will train until train_error hasn't improved in 10 rounds.

[2] train-error:0.020573 
[3] train-error:0.020573 
[4] train-error:0.011078 
[5] train-error:0.009320 
[6] train-error:0.009320 
[7] train-error:0.008440 
[8] train-error:0.008440 
[9] train-error:0.001231 
[10]    train-error:0.003693 
[11]    train-error:0.001231 
[12]    train-error:0.001231 
[13]    train-error:0.001231 
[14]    train-error:0.001231 
[15]    train-error:0.001231 
[16]    train-error:0.000703 
[17]    train-error:0.000000 
[18]    train-error:0.000000 
[19]    train-error:0.000000 
[20]    train-error:0.000000 
[21]    train-error:0.000000 
[22]    train-error:0.000000 
[23]    train-error:0.000000 
[24]    train-error:0.000000 
[25]    train-error:0.000000 
[26]    train-error:0.000000 
[27]    train-error:0.000000 
Stopping. Best iteration:
[17]    train-error:0.000000
Confusion Matrix and Statistics

   X2
X1     0    1
  0 1246    0
  1    0 1191
                                     
               Accuracy : 1          
                 95% CI : (0.9985, 1)
    No Information Rate : 0.5113     
    P-Value [Acc > NIR] : < 2.2e-16  
                                     
                  Kappa : 1          
                                     
 Mcnemar's Test P-Value : NA         
                                     
            Sensitivity : 1.0000     
            Specificity : 1.0000     
         Pos Pred Value : 1.0000     
         Neg Pred Value : 1.0000     
             Prevalence : 0.5113     
         Detection Rate : 0.5113     
   Detection Prevalence : 0.5113     
      Balanced Accuracy : 1.0000     
                                     
       'Positive' Class : 0          
                                     

El entrenamiento se ha detenido después de 17 iteraciones y ha producido un modelo con una precisión del 100%. Nada mal, aunque vale la pena mencionar que podría haber sobre ajustado.

Para concluir

En este artículo hemos revisado, de manera general, la implementación para R del algoritmo XGBoost. En el proceso, también dimos un vistazo al proceso para preparar datos con formatos no convencionales para ser usados en este algoritmo.

Esta revisión no ha sido exhaustiva, hay algunos temas que es importante estudiar para obtener mejores resultados al usar XGBoost:

Si quieres conocer más sobre estos temas un buen punto de partida es la documentación de XGBoost.


Consultas, dudas, comentarios y correcciones son bienvenidas:

El código y los datos usados en este documento se encuentran en Github: