Introducción (simple) a caret

Objetivos

• Introducir algunos conceptos básicos del enfoque del Aprendizaje Automático
• Mostrar el framework caret para automatizar algunas tareas del entrenamiento

El problema

Nuestro problema central es, poder realizar un modelo que logre prededir los ingresos de la ocupación principal (p21) en la EPH del segundo trimestre del 2015.

Puede notarse que se trata de un problema bastante amigable, por así decirlo, al enfoque de Machine Learning:

• tenemos un conjunto de casos en los que desconocemos nuestra variable \(Y\)
• queremos predecirla

Preprocesando los datos

Lo primero que tenemos que hacer es importar las librerías con las que vamos a trabajar:

library(caret)

Loading required package: lattice
Loading required package: ggplot2

library(tidyverse)

-- Attaching packages --------------------------------------- tidyverse 1.2.1 --
v tibble  2.1.3     v purrr   0.2.5
v tidyr   0.8.2     v dplyr   0.8.3
v readr   1.3.1     v stringr 1.3.1
v tibble  2.1.3     v forcats 0.3.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
x purrr::lift()   masks caret::lift()

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'))

Existen en nuestro dataset, datos que no contestaron ingresos. Son datos perdidos y tenemos que resolver qué hacer con ellos. En este ejemplo vamos a eliminarlos. Esta opción está lejos de ser la óptima, pero la seleccionamos para simplificar el problema y la exposición. En caso de que les interese el tema de imputación de missing data aquí y aquí pueden encontrar dos aplicaciones de Machine Learning al problema (y sobre el mismo dataset).

df_imp <- data %>%
        filter(imp_inglab1==1) %>%
        select(-imp_inglab1)
df_train <- data %>%
        filter(imp_inglab1==0) %>%
        select(-imp_inglab1) %>%
        mutate(p21 = case_when(
                        p21==0 ~ 100,
                        TRUE ~ p21))

Algunas cosas a notar

Por un lado, vemos que encadenamos unas cuántas operaciones mediante un operador (%>%) llamado pipe. El pipe es un símbolo que relaciona dos entidades. Dicho en forma más simple, el pipe de R está en familia con otros operadores más convencionales, como +, - o /. Y al igual que los otros operadores, entrega un resultado en base a los operandos que recibe.

Ahora bien… ¿Para qué sirve? En resumidas cuentas, hace que el código necesario para realizar una serie de operaciones de transformación de datos sea mucho más simple de escribir y de interpretar.

Repasemos la primer secuencia

• filtramos los datos con algún perdido (%>% filter(imp_inglab==1))
• eliminamos la columna identificadora de los casos perdidos (select(-imp_inglab))

Estimando el error de generalización

Recordemos: tenemos muchas formas de estimar el error de generalización (train-test split, cross validation, bootstrap). Usaremos una estrategia de validación cruzada. Vamos a generar los índices mediante caret.

Vamos a tener que fijar dos estrategias de estimación del error: la primera para estimar los hiperparámetros de los modelos y la segunda para la estimación final del error de generalización. En ambos casos, utilizaremos validación cruzada, pero sobre dos muestras diferentes.

Primero, fijamos la semilla aleatoria (para asegurarnos la posibilidad de replicabilidad)

set.seed(957)

Podemos usar la función createFolds() para generar los índices. Aquí, pas

cv_index <- createFolds(y = df_train$p21,
                        k=5,
                        list=TRUE,
                        returnTrain=TRUE)

Aquí usamos tres argumentos:

• y = df_train$p21, es el vector de resultados. En nuestro caso, los ingresos de la ocupación principal
• k=5, es la cantidad de grupos para realizar la validación cruzada
• returnTrain=TRUE, le decimos que lo que nos devuelva, sean las posiciones de correspondientes a los datos de entrenamiento en cada posición

Finalmente, especificamos el diseño de remuestreo mediante la función trainControl:

fitControl <- trainControl(
        index=cv_index,
        method="cv",
        number=5)

Entrenando modelos (train())

Tenemos listo nuestro esquema de remuestreo. Podemos pasar a entrenar nuestro primer modelo. Para ello haremos uso extensivo de la función train(). La misma puede usarse para

• evaluar mediante remuestreo el efecto de cada hiperparámetro en la performance
• elegir el modelo “óptimo” (la mejor combinación de parámetros)
• estimar la performance del modelo

Primero, debemos elegir el modelo para entrenar. Actualmente, caret dispone de 238 modelos disponibles. Puede consultarse la seccion correspondiente del sitio para mayores detalles. También, llegado el caso, podrían usarse modelos ad-hoc definidos por el usuario.

Comencemos con un modelo simple, pero efectivo: una regresión lineal. Como podrán ver en el sitio, cada modelo puede ser estimado por diferentes implementaciones en diferentes paquetes. Nosotros usaremos la implementación de r-base lm() por simplicidad.

Entrenemos una regresión lineal con caret: comencemos con un modelo simple, sexo y edad.

lm_p21 <- train(p21 ~ ch04 + ch06, data = df_train, 
                 method = "lm", 
                 trControl = fitControl)
lm_p21

Linear Regression 

19406 samples
    2 predictor

No pre-processing
Resampling: Cross-Validated (5 fold) 
Summary of sample sizes: 15523, 15525, 15526, 15525, 15525 
Resampling results:

  RMSE      Rsquared    MAE     
  5600.527  0.05340693  3836.812

Tuning parameter 'intercept' was held constant at a
 value of TRUE

Veamos los coeficientes…

lm_p21$finalModel

Call:
lm(formula = .outcome ~ ., data = dat)

Coefficients:
(Intercept)    ch04Mujer         ch06  
    5022.91     -2019.82        70.68  

¿Qué se puede ver?

Veamos, ahora, un modelo más complejo:

lm_p21_b <- train(p21 ~ ., data = df_train, 
                 method = "lm", 
                 trControl = fitControl)

lm_p21_b

Linear Regression 

19406 samples
   25 predictor

No pre-processing
Resampling: Cross-Validated (5 fold) 
Summary of sample sizes: 15523, 15525, 15526, 15525, 15525 
Resampling results:

  RMSE      Rsquared   MAE     
  4019.713  0.5127815  2476.719

Tuning parameter 'intercept' was held constant at a
 value of TRUE

Los modelos de machine learning tienen ciertos parámetros que deben ser seleccionados antes de estimar el modelo, propiamente dicho: se llaman hiperparámetros. Si bien la regresión lineal no es estrictamente hablando un modelo de machine learning (aunque muches lo consideran como tal) sí tiene algo que se le parece bastante a un hiperparámetro… la existencia de un intercepto. En efecto, nosotros estimamos un modelo de la siguiente forma:

\(y_{i} = \beta_{0} + \sum_{p=1}^P \beta_{p} X_{i}\)

Pero podríamos haber estimado

\(y_{i} = \sum_{p=1}^P \beta_{p} X_{i}\)

Más allá de la discusión sobre si la regresión es ML o no, lo interesante es ver que la decisión sobre el entrenamiento de un modelo lineal con intercepto o no, es una decisión que se toma antes de entrenar el modelo propiamente dicho.

Ahora bien, vamos a buscar otro modelo con mejores hiperparámetros para tunear: un árbol de decisión. Si bien, lo vamos a ver en detalle la próxima clase, vamos a revisar su implementación en caret.

Tuneando hiperparámetros…

Podemos entonces, comparar la performance de un modelo con y sin hiperparámetros. Para ello, primero tenemos que construir la grilla de hiperparámetros.

grid <- expand.grid(maxdepth=c(1, 2, 4, 8, 16))

Y volvemos a entrenar el modelo:

cart_p21 <- train(p21 ~ . , 
                 data = df_train, 
                 method = "rpart2", 
                 trControl = fitControl,
                 tuneGrid =grid)
cart_p21

CART 

19406 samples
   25 predictor

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

  maxdepth  RMSE      Rsquared   MAE     
   1        5254.093  0.1672134  3380.275
   2        5119.021  0.2100596  3288.701
   4        4874.119  0.2837009  3146.744
   8        4672.406  0.3415589  2923.876
  16        4647.562  0.3485327  2901.520

RMSE was used to select the optimal model using
 the smallest value.
The final value used for the model was maxdepth = 16.

En este caso, hemos realizado una búsqueda exhaustiva, es decir, hemos recorrido la totalidad de la grilla de hiperparámetros y hemos seleccionado el mejor modelo. Como puede verse, esto ha llevado un tiempo de cómputo nada despreciable.

Es por ello que existe una segunda opción…

Random search

En este caso, en lugar de realizar una búsqueda exhaustiva, podemos reducir notablemente el tiempo de cómputo buscando en una muestra aleatoria de la grilla de hiperparámetros. Para esto, solamente debemos agregar un parámetro en el objeto fitControl:

fitControl_rand <- trainControl(
        index=cv_index, 
        method="cv",
        number=5,
        search = 'random')

Y volvemos a entrenar el modelo:

cart_p21_rand <- train(p21 ~ ., data = df_train, 
                 method = "rpart2", 
                 trControl = fitControl_rand,
                 tuneLength = 2)
cart_p21_rand

CART 

19406 samples
   25 predictor

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

  maxdepth  RMSE      Rsquared   MAE     
  5         4791.211  0.3078338  3019.845
  8         4647.562  0.3485327  2901.520

RMSE was used to select the optimal model using
 the smallest value.
The final value used for the model was maxdepth = 8.

Seleccionando el mejor modelo

Una vez finalizado el proceso de tunning de los hiperparámetros, podemos proceder a elegir cuál es el mejor modelo.

cart_p21

CART 

19406 samples
   25 predictor

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

  maxdepth  RMSE      Rsquared   MAE     
   1        5254.093  0.1672134  3380.275
   2        5119.021  0.2100596  3288.701
   4        4874.119  0.2837009  3146.744
   8        4672.406  0.3415589  2923.876
  16        4647.562  0.3485327  2901.520

RMSE was used to select the optimal model using
 the smallest value.
The final value used for the model was maxdepth = 16.

Podemos persistir el modelo en disco (si quisiéramos):

saveRDS(cart_p21, '../models/p21_cart.rds')

Podemos realizar un plot del efecto de los hiperparámetros:

ggplot(cart_p21)

Existen diferentes métricas de selección, las cuales deben ser definidas en la función train, usando el argumento selectionFunction que puede tomar tres valores:

• "best": se selecciona el mejor modelo con la menor métrica de error (la que usaremos aquí)
• "oneSE": utiliza la regla de “un desvío estándar” de Breiman et al (1986)
• "tolerance; que busca seleccionar el modelo menos complejo dentro de un margen de tolerancia respecto al mejor modelo

También podrían definirse métodos ad-hoc para esta selección.

cart_p21$bestTune

¿Cuál es el mejor modelo (en términos absolutos)?

Realizando la evaluación final

Una vez que hemos seleccionado el mejor modelo, podemos pasar a la evaluación final y persistimos el modelo para reutilizarlo en otras aplicaciones.

Primero, tenemos que volver a generar un esquema de validación cruzada:

set.seed(7412)
cv_index_final <- createFolds(y = df_train$p21,
                        k=5,
                        list=TRUE,
                        returnTrain=TRUE)
fitControl_final <- trainControl(
        indexOut=cv_index_final, 
        method="cv",
        number=5)

Y entrenamos una vez más:

cart_final<-train(p21 ~ ., data = df_train,
                method = "rpart2", 
                trControl = fitControl_final, 
                tuneGrid = cart_p21$bestTune,
                metric='RMSE')
#saveRDS(rf_final, '../models/rf_final.RDS')
cart_final

CART 

19406 samples
   25 predictor

No pre-processing
Resampling: Cross-Validated (5 fold) 
Summary of sample sizes: 15524, 15526, 15525, 15525, 15524 
Resampling results:

  RMSE      Rsquared   MAE     
  4608.328  0.3594745  2888.506

Tuning parameter 'maxdepth' was held constant at a
 value of 16

Vemos entonces que el modelo seleccionado performa con un \(R^2=0.36\) y un \(RMSE=4603\). Solamente nos queda entrenar el modelos sobre la totalidad del dataset de entrenamiento:

cart_final_f<-train(p21~., data=df_train,
                  method = "rpart2",
                  tuneGrid = cart_p21$bestTune)
cart_final_f

CART 

19406 samples
   25 predictor

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 19406, 19406, 19406, 19406, 19406, 19406, ... 
Resampling results:

  RMSE      Rsquared  MAE     
  4679.353  0.344868  2915.385

Tuning parameter 'maxdepth' was held constant at a
 value of 16

Obteniendo las predicciones finales

El último paso es obtener las predicciones finales (es decir, nuestras imputaciones). De forma interesante, podemos utilizar lso datos perdidos como datos “nuevos” y desconocidos.

Es decir que, finalmente, habremos realizado una imputación de datos perdidos. Para ello, llamamos a la función predict() que toma como primer argumento al objeto que contiene al modelo final y como segundo argumento el data.frame con los datos a imputar:

y_preds_cart <- predict(cart_final_f, df_imp)

Comparemos, ahora, las distribuciones de datos imputados por el INDEC (mediante el método Hot Deck) y los que hemos imputado con rpart2. Para ello, organizamos todo en un data frame que, luego, llevamos al formato tidy.

preds <- cbind(y_preds_cart,
               df_imp$p21
)
colnames(preds) <- c('CART', 'Hot_Deck')
preds <- preds %>% as.data.frame() %>% gather(model, value)

Finalmente, ploteamos un gráfico de densidad para comparar las distribuciones de los casos imputados con ambos métodos.

ggplot(preds) +
        geom_density(aes(x=value, fill=model), alpha=0.5)

ggplot(preds) +
        geom_histogram(aes(x=value, fill=model), alpha=0.5,
                       bins=50)

Práctica independiente: entrenando un árbol para predecir la no respuesta en ingresos

La idea ahora es que ustedes entrenen otro modelo. Vamos a entrenar y evaluar otro modelo en otro problema. Tratemos de predecir la probabilidad de que una persona no conteste ingresos. Usemos para ello un arbol de decisión.

rr ###