Ensamble Learning - Bagging, Random Forest y AdaBoost

Objetivos

• Mostrar la implementación en caret de diferentes algoritmos de Ensamble Learning para problemas de clasificación
• Mostrar los problemas de los datasets desbalanceados y algunos métodos posibles para lidiar con el problema

Carga de librerías y datos

library(caret)
library(tidyverse)
library(rpart)

Luego, cargamos los datos y formateamos un poco algunas etiquetas:

load('../data/EPH_2015_II.RData')

data$pp03i<-factor(data$pp03i, labels=c('1-SI', '2-No', '9-NS'))



data$intensi<-factor(data$intensi, labels=c('1-Sub_dem', '2-SO_no_dem', 
                                            '3-Ocup.pleno', '4-Sobreoc',
                                            '5-No trabajo', '9-NS'))

data$pp07a<-factor(data$pp07a, labels=c('0-NC',
                                        '1-Menos de un mes',
                                        '2-1 a 3 meses',
                                        '3-3 a 6 meses',
                                        '4-6 a 12 meses',
                                        '5-12 a 60 meses',
                                        '6-Más de 60 meses',
                                        '9-NS'))



data <- data %>%
        mutate(imp_inglab1=factor(imp_inglab1, labels=c('non_miss','miss')))

Ahora, nuestra variable a predecir es el indicador imp_inglab. Por ende, elimonamos la \(p21\).

df_train <- data %>%
        select(-p21)

Lo primero que vamos a hacer es crear una partición de datos:

set.seed(1234)
tr_index <- createDataPartition(y=df_train$imp_inglab1,
                                p=0.8,
                                list=FALSE)

Y generamos dos datasets:

train <- df_train[tr_index,]
test <- df_train[-tr_index,]

Random Forest

La idea ahora es poder entrenar un modelo random forest. Pero… lo van a hacer ustedes. Para ello, deberán consultar la documentación de caret y buscar cuál es la mejor función para entrenar el modelo y sus hiperparámetros.

Pista: buscar el método ranger.

Usaremos la partición original entre Train-Test. Recordemos los pasos a seguir

1. Generar el esquema de validación cruzada
2. Generar el grid de hiperparámetros
3. Tunear el modelo
4. Seleccionar y estimar el modelo final
5. Validar y evaluar contra test-set

set.seed(5699)
cv_index_rf <- createFolds(y=train$imp_inglab1,
                        k=5,
                        list=TRUE,
                        returnTrain=TRUE)

fitControlrf <- trainControl(
        index=cv_index_rf,
        method="cv",
        number=5,
        summaryFunction = twoClassSummary,
        classProbs=TRUE
        )

# Generar grid

grid_rf <- expand.grid(mtry=c(5,10,15, 25),
                       min.node.size=c(5,10,15,20),
                       splitrule='gini'
        )

# Tunear
t0 <- proc.time()
rf_fit_orig <-  train(imp_inglab1 ~ . , 
                 data = train, 
                 method = "ranger", 
                 trControl = fitControlrf,
                 tuneGrid = grid_rf,
                 metric='ROC'
)
proc.time() -  t0

saveRDS(rf_fit_orig, '../models/rf_fit_orig.RDS')

rf_fit_orig

Random Forest 

19519 samples
   25 predictor
    2 classes: 'non_miss', 'miss' 

No pre-processing
Resampling: Cross-Validated (5 fold) 
Summary of sample sizes: 15616, 15615, 15615, 15615, 15615 
Resampling results across tuning parameters:

  mtry  min.node.size  ROC        Sens       Spec      
   5     5             0.7872581  0.9958132  0.04006167
   5    10             0.7867366  0.9963285  0.03605447
   5    15             0.7872137  0.9966506  0.03730697
   5    20             0.7862695  0.9969082  0.03129758
  10     5             0.7885261  0.9781643  0.14671880
  10    10             0.7902176  0.9789372  0.13645503
  10    15             0.7902504  0.9799678  0.13144689
  10    20             0.7907318  0.9820290  0.12518875
  15     5             0.7872475  0.9691465  0.18202139
  15    10             0.7884391  0.9712721  0.17325855
  15    15             0.7901788  0.9734622  0.16775167
  15    20             0.7899749  0.9751369  0.15998883
  25     5             0.7849605  0.9598068  0.20831051
  25    10             0.7876549  0.9638647  0.20155050
  25    15             0.7880576  0.9661192  0.19253766
  25    20             0.7890911  0.9691465  0.18427640

Tuning parameter 'splitrule' was held constant at a value of gini
ROC was used to select the optimal model using the largest value.
The final values used for the model were mtry = 10, splitrule =
 gini and min.node.size = 20.

confusionMatrix(predict(rf_fit_orig, test), test$imp_inglab1)

Confusion Matrix and Statistics

          Reference
Prediction non_miss miss
  non_miss     3787  871
  miss           94  127
                                          
               Accuracy : 0.8022          
                 95% CI : (0.7908, 0.8133)
    No Information Rate : 0.7954          
    P-Value [Acc > NIR] : 0.1241          
                                          
                  Kappa : 0.1449          
                                          
 Mcnemar's Test P-Value : <2e-16          
                                          
            Sensitivity : 0.9758          
            Specificity : 0.1273          
         Pos Pred Value : 0.8130          
         Neg Pred Value : 0.5747          
             Prevalence : 0.7954          
         Detection Rate : 0.7762          
   Detection Prevalence : 0.9547          
      Balanced Accuracy : 0.5515          
                                          
       'Positive' Class : non_miss        
                                          

• ¿Qué pueden decir de este modelo?
• ¿Funciona mejor que un árbol individual?
• ¿Qué sucede con su performance entre los miss

Datasets desbalanceados… ¿qué hacer?

Efectivamente, parte del problema de este dataset es que tenemos una variable dependiente sumamente desbalanceada:

ggplot(train) + 
        geom_bar(aes(x=imp_inglab1)) +
        theme_minimal()

Matriz de confusión

Ya estamos familiarizades con la matriz de confusión. Básicamente, es una forma rápida de comparar las predicciones del modelo con los valores reales.

Así, la métrica más simple es la siguiente:

\(Accuracy = \frac{TP + TN}{N} = \frac{TP + TN}{TP + TN + FP + FN}\)

Simplemente, calcula la proporción de casos bien clasificados (\(TP + TN\)) sobre el total. Pero esta métrica tiene algunos problemas. El más importante es que tiende a ser una medida engañosa de la performance de un modelo en un dataset desbalanceado.

Un primer punto importante es poder seleccionar una métrica adecuada para evaluar nuestro modelo. Existen dos grandes familias de medidas:

• Dependiente del umbral: Esto incluye métricas como precisión, recall y puntaje F1, que requieren una matriz de confusión para ser calculada usando un límite duro en las probabilidades pronosticadas. Estas métricas suelen ser bastante pobres en el caso de clases desequilibradas, ya que el software estadístico utiliza de manera inapropiada un umbral predeterminado de 0,50, lo que hace que el modelo prediga que todas las observaciones pertenecen a la clase mayoritaria.

\(Precision = \frac{TP}{TP+FP}\) Ratio entre verdaderos positivos y los que el modelo clasifica como positivos. ¿Qué proporción de los que clasifcamos como positivos son efectivamente positivos?

\(Recall = \frac{TP}{TP+FN}\) También llamada, \(Sensitividad\), Ratio entre los positivos bien clasificados y el total de los que son positivos en realidad. ¿Qué proporción de los que efectivamente positivos hemos clasificado bien?

\(F1-Score = 2 \times \frac{Precision \times Recall}{Recall + Precision}\)

• Independiente del umbral: Esto incluye métricas como área bajo la curva ROC (AUC), que cuantifica la tasa positiva verdadera en función de la tasa de falsos positivos para una variedad de umbrales de clasificación. Otra forma de interpretar esta métrica es la probabilidad de que una instancia positiva aleatoria tenga una probabilidad estimada más alta que una instancia negativa aleatoria.

En líneas generales, existen cuatro formas generales de lidiar con el problema de las clases desbalanceadas:

• Usar modelos ponderados
• Down-sampling es decir, submuestrear las instancias sobrerrepresentadas
• Up-sampling, es decir, sobremuestrar las instancias subrrepresentadas
• Synthetic minority sampling technique (SMOTE): hace down-sampling de la clase mayoritaria y “sintetiza” nuevas instancias minoritarias mediante interpolación

Utilizar pesos en el modelo

Para esto podemos utilizar el argumento weights en la función train (esto supone que el modelo seleccionado puede manejar pesos).

model_weights <- ifelse(train$imp_inglab1 == "non_miss",
                        (1/table(train$imp_inglab1)[1]) * 0.5,
                        (1/table(train$imp_inglab1)[2]) * 0.5)

fitControlrf$seeds <- rf_fit_orig$control$seeds

Y entrenamos el modelo:

t0 <- proc.time()
rf_fit_wei <-  train(imp_inglab1 ~ . , 
                 data = train, 
                 method = "ranger", 
                 trControl = fitControlrf,
                 tuneGrid = grid_rf,
                 weights = model_weights,
                 metric='ROC'
)
proc.time() -  t0

    user   system  elapsed 
2863.763   13.410  481.441 

saveRDS(rf_fit_wei, '../models/rf_fit_wei.RDS')

Métodos de remuestreo

Probemos, ahora, métodos de remuestreo: up sampling y down sampling.

# Nuevo fitControl con upsample
fitControlrf_imb <- trainControl(
        index=cv_index_rf,
        method="cv",
        number=5,
        summaryFunction = twoClassSummary,
        classProbs=TRUE,
        sampling = 'up'
        )

# Tunear upsample

fitControlrf_imb$seeds <- rf_fit_orig$control$seeds

t0 <- proc.time()
rf_fit_up <-  train(imp_inglab1 ~ . , 
                 data = train, 
                 method = "ranger", 
                 trControl = fitControlrf_imb,
                 tuneGrid = grid_rf,
                 metric='ROC'
)
proc.time() -  t0

saveRDS(rf_fit_up, '../models/rf_fit_up.RDS')

# Tunear downsample

fitControlrf_imb$sampling <- "down"
fitControlrf_imb$seeds <- rf_fit_orig$control$seeds


t0 <- proc.time()
rf_fit_down <-  train(imp_inglab1 ~ . , 
                 data = train, 
                 method = "ranger", 
                 trControl = fitControlrf_imb,
                 tuneGrid = grid_rf,
                 metric='ROC'
)
proc.time() -  t0

saveRDS(rf_fit_down, '../models/rf_fit_down.RDS')

Comparemos todos los modelos entrenados hasta ahora. Pero ya que estamos, incorporemos el modelo de árbol simple que habíamos entrenado la semana pasada:

cart_final <- train(imp_inglab1 ~ . , 
                 data = train, 
                 method = "rpart2", 
                 tuneGrid = data.frame(maxdepth=6),
                 control = rpart.control(cp=0.000001)
)

Colocamos todos los modelos en una lista:

model_list <- list(cart = cart_final,
                   rf_original = rf_fit_orig,
                   rf_weighted = rf_fit_wei,
                   rf_down = rf_fit_down,
                   rf_up = rf_fit_up)

Extraigamos algunas métricas resumen de las matrices de confusión:

extract_conf_metrics <- function(model, data, obs){
        preds <- predict(model, data)        
        c<-confusionMatrix(preds, obs)
        results <- c(c$overall[1], c$byClass)
        return(results)
}

model_metrics <- model_list %>%
        map(extract_conf_metrics, data=test, obs = test$imp_inglab1) %>%
        do.call(rbind,.) %>%
        as.data.frame()

Puede verse, entonces, cómo la performace de los diferentes modelos cambia en los diferentes tunnings…

model_metrics %>%
        select('Accuracy', 'Precision', 'Recall') %>%
        t()

               cart rf_original rf_weighted   rf_down     rf_up
Accuracy  0.7913507   0.8022136   0.7427752 0.6632507 0.7675753
Precision 0.8076510   0.8130099   0.8819081 0.8951271 0.8642270
Recall    0.9683071   0.9757794   0.7812419 0.6531822 0.8397320

Observando específicamente las métricas de \(Precision\), \(Recall\) y \(Accuracy\) puede verse que el modelo de random forest original parece mejorar las diferentes métricas.

AdaBoost

Por último, vamos a entrenar un modelo de Adaboost.

fitControlrf$seeds <- rf_fit_orig$control$seeds
#grid_ada <- expand.grid(mfinal=c(100,150), maxdepth=c(10,20), coeflearn='Breiman')
grid_ada <- expand.grid(nIter=c(100,150), method='Adaboost.M1')

t0 <- proc.time()
adaboost_fit_orig <-  train(imp_inglab1 ~ . , 
                 data = train, 
                 method = "adaboost", 
                 trControl = fitControlrf,
                 tuneGrid = grid_ada,
                 metric='ROC'
)

adaboost_fit_orig

AdaBoost Classification Trees 

19519 samples
   25 predictor
    2 classes: 'non_miss', 'miss' 

No pre-processing
Resampling: Cross-Validated (5 fold) 
Summary of sample sizes: 15616, 15615, 15615, 15615, 15615 
Resampling results across tuning parameters:

  nIter  ROC        Sens       Spec     
  100    0.7661744  0.9148470  0.3064605
  150    0.7683175  0.9162641  0.3037039

Tuning parameter 'method' was held constant at a value
 of Adaboost.M1
ROC was used to select the optimal model using the
 largest value.
The final values used for the model were nIter = 150 and
 method = Adaboost.M1.

Comparando los modelos generados

Volvamos a armar la lista de modelos:

model_list <- list(cart = cart_final,
                   rf_original = rf_fit_orig,
                   rf_weighted = rf_fit_wei,
                   rf_down = rf_fit_down,
                   rf_up = rf_fit_up,
                   adaboost = adaboost_fit_orig)

Y, finalmente, generemos una tabla que contenga todo:

model_metrics <- model_list %>%
        map(extract_conf_metrics, data=test, obs = test$imp_inglab1) %>%
        do.call(rbind,.) %>%
        as.data.frame() %>%
        t()

model_metrics

                          cart rf_original rf_weighted
Accuracy             0.7913507   0.8022136   0.7427752
Sensitivity          0.9683071   0.9757794   0.7812419
Specificity          0.1032064   0.1272545   0.5931864
Pos Pred Value       0.8076510   0.8130099   0.8819081
Neg Pred Value       0.4557522   0.5746606   0.4108258
Precision            0.8076510   0.8130099   0.8819081
Recall               0.9683071   0.9757794   0.7812419
F1                   0.8807124   0.8869891   0.8285285
Prevalence           0.7954499   0.7954499   0.7954499
Detection Rate       0.7702398   0.7761836   0.6214388
Detection Prevalence 0.9536790   0.9547038   0.7046526
Balanced Accuracy    0.5357568   0.5515170   0.6872142
                       rf_down     rf_up  adaboost
Accuracy             0.6632507 0.7675753 0.7893011
Sensitivity          0.6531822 0.8397320 0.9149704
Specificity          0.7024048 0.4869739 0.3006012
Pos Pred Value       0.8951271 0.8642270 0.8357261
Neg Pred Value       0.3424524 0.4386282 0.4761905
Precision            0.8951271 0.8642270 0.8357261
Recall               0.6531822 0.8397320 0.9149704
F1                   0.7552510 0.8518035 0.8735547
Prevalence           0.7954499 0.7954499 0.7954499
Detection Rate       0.5195737 0.6679647 0.7278131
Detection Prevalence 0.5804468 0.7729043 0.8708752
Balanced Accuracy    0.6777935 0.6633530 0.6077858

Puede verse que cart, rf_original y adaboost presentan valores similares en accuracy. A su vez, al observar la relación entre precision y recall puede verse que el modelo que mejor balancea los dos es rf_up. Cuál elegir dependerá del objetivo del análisis.