Introducción a Caret
Introducción
En este notebook de R haremos una breve introducción al paquete caret (Classification And REgression Training). Toda la información expuesta aquí está basada, principalmente, en el libro Applied Predictive Modeling, de Max Kuhn, el creador de caret, y de la referencia web del paquete.
Además de conocer las funciones básicas de caret, utilizaremos algunos paquetes adicionales. Ejecutando el código siguiente (quitando el comentario) los instalaremos:
#install.packages(c("caret", "AppliedPredictiveModeling", "readr", "e1071", "RANN", "corrplot", "ggplot2", "plyr", "dplyr","ranger", "party", "mboost", "doMC", "DMwR", "ROSE", "caTools", "kernlab"))
library(readr)
library(caret)
library(e1071)
library(ggplot2)
library(corrplot)
library(plyr)
library(dplyr)
library(ranger)
library(party)
library(mboost)
library(DMwR)
library(ROSE)
library(caTools)
library(kernlab)caret es un paquete orientado a simplificar el proceso de creación de modelos predictivos. Permite separar datos, preprocesarlos, seleccionar variables, afinar modelos y estimar la importancia de las variables, entre otras muchas funciones. Entre todas las aportaciones que realiza, para mí la más importante de todas es proporcionar una interfaz única a una gran cantidad de paquetes adicionales de modelado.
La estructura de la charla será la siguiente:
Funcionalidades básicas
- Preprocesado
- Separación de datos
- Entrenamiento, evaluación y ajuste de modelos
Funcionalidades avanzadas
- Procesamiento en paralelo
- Submuestreo para compensar desequilibrios entre clases
- Importancia de las variables
- Selección de variables mediante algoritmo genético
1. Preprocesado de datos
La gran mayoría del preprocesado es accesible mediante la función preProcess.
1.1 Transformaciones sobre variables individuales
- Centrado: la sustracción de la media.
- Normalización: cálculo de \(z = \frac{x-\mu}{\sigma}\), donde \(\mu\) y \(\sigma\) son la media y desviación estándar de la variable, respectivamente.
- deskew (desinclinación): corrección de la curtosis en una variable. La mejor manera de hacerla mediante
caretes aplicando una transformación Box Cox: \[ x^*=\cases{\frac{x^\lambda-1}{\lambda}\text{if }\lambda\neq0\\\log(x)\text{ if }\lambda=0} \] El paquete encuentra el valor óptimo del parámetro y lo devuelve, simplificando en gran medida esta parte.
#Utilizaremos un dataset sintético de registros de sesiones de compras de clientes en un portal web
dat <- read_csv(file = "sesiones.csv")## Parsed with column specification:
## cols(
## clicks = col_integer(),
## meanTime = col_double(),
## sessionTime = col_double(),
## articlesSeen = col_integer(),
## articlesBought = col_integer(),
## categories = col_integer(),
## money = col_double(),
## userAgent = col_character()
## )summary(dat)## clicks meanTime sessionTime articlesSeen
## Min. : 0.0 Min. : 0.00193 Min. : 0.042 Min. : 1.000
## 1st Qu.: 16.0 1st Qu.: 16.92835 1st Qu.: 136.084 1st Qu.: 4.000
## Median : 34.0 Median : 67.02895 Median : 287.745 Median : 7.000
## Mean : 87.2 Mean : 63.30748 Mean : 557.231 Mean : 8.563
## 3rd Qu.: 84.0 3rd Qu.: 97.27688 3rd Qu.: 660.216 3rd Qu.:12.000
## Max. :811.0 Max. :227.25200 Max. :6311.810 Max. :35.000
## articlesBought categories money userAgent
## Min. : 0.000 Min. : 1.000 Min. : 0.00 Length:62122
## 1st Qu.: 0.000 1st Qu.: 3.000 1st Qu.: 0.00 Class :character
## Median : 2.000 Median : 6.000 Median : 33.85 Mode :character
## Mean : 4.431 Mean : 7.401 Mean : 64.65
## 3rd Qu.: 6.000 3rd Qu.:11.000 3rd Qu.: 92.85
## Max. :30.000 Max. :32.000 Max. :853.28dat$userAgent <- factor(dat$userAgent)
# Centrado ----
pre <- preProcess(x = dat, method = c("center"))
pre## Created from 62122 samples and 8 variables
##
## Pre-processing:
## - centered (7)
## - ignored (1)summary(predict(object = pre, newdata = dat))## clicks meanTime sessionTime articlesSeen
## Min. :-87.195 Min. :-63.306 Min. :-557.2 Min. :-7.563
## 1st Qu.:-71.195 1st Qu.:-46.379 1st Qu.:-421.1 1st Qu.:-4.563
## Median :-53.195 Median : 3.721 Median :-269.5 Median :-1.563
## Mean : 0.000 Mean : 0.000 Mean : 0.0 Mean : 0.000
## 3rd Qu.: -3.195 3rd Qu.: 33.969 3rd Qu.: 103.0 3rd Qu.: 3.437
## Max. :723.805 Max. :163.945 Max. :5754.6 Max. :26.437
## articlesBought categories money
## Min. :-4.431 Min. :-6.401 Min. :-64.65
## 1st Qu.:-4.431 1st Qu.:-4.401 1st Qu.:-64.65
## Median :-2.431 Median :-1.401 Median :-30.81
## Mean : 0.000 Mean : 0.000 Mean : 0.00
## 3rd Qu.: 1.569 3rd Qu.: 3.599 3rd Qu.: 28.20
## Max. :25.569 Max. :24.599 Max. :788.63
## userAgent
## Chrome Generic MacOSX: 7619
## Chrome Generic Win10 :24125
## Chrome Generic Win7 : 7617
## Firefox Generic Linux: 7581
## Firefox Generic Win10:15180
## # Centrado y escalado ----
pre <- preProcess(x = dat, method = c("center", "scale"))
pre## Created from 62122 samples and 8 variables
##
## Pre-processing:
## - centered (7)
## - ignored (1)
## - scaled (7)summary(predict(object = pre, newdata = dat))## clicks meanTime sessionTime articlesSeen
## Min. :-0.71199 Min. :-1.54333 Min. :-0.7761 Min. :-1.2369
## 1st Qu.:-0.58134 1st Qu.:-1.13068 1st Qu.:-0.5866 1st Qu.:-0.7462
## Median :-0.43437 Median : 0.09073 Median :-0.3753 Median :-0.2556
## Mean : 0.00000 Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.:-0.02609 3rd Qu.: 0.82814 3rd Qu.: 0.1434 3rd Qu.: 0.5621
## Max. : 5.91021 Max. : 3.99681 Max. : 8.0150 Max. : 4.3235
## articlesBought categories money
## Min. :-0.7911 Min. :-1.1281 Min. :-0.7282
## 1st Qu.:-0.7911 1st Qu.:-0.7756 1st Qu.:-0.7282
## Median :-0.4341 Median :-0.2469 Median :-0.3470
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.2801 3rd Qu.: 0.6343 3rd Qu.: 0.3176
## Max. : 4.5651 Max. : 4.3352 Max. : 8.8821
## userAgent
## Chrome Generic MacOSX: 7619
## Chrome Generic Win10 :24125
## Chrome Generic Win7 : 7617
## Firefox Generic Linux: 7581
## Firefox Generic Win10:15180
## # Transformación BoxCox ----
pre <- preProcess(x = as.data.frame(dat), method = c("center", "scale", "BoxCox"))
pre## Created from 62122 samples and 8 variables
##
## Pre-processing:
## - Box-Cox transformation (4)
## - centered (7)
## - ignored (1)
## - scaled (7)
##
## Lambda estimates for Box-Cox transformation:
## 0.6, 0.1, 0.3, 0.2boxDat <- predict(object = pre, newdata = as.data.frame(dat))
apply(X = dat[,1:7], MARGIN = 2, FUN = skewness)## clicks meanTime sessionTime articlesSeen articlesBought
## 2.02035681 0.02649531 2.71775159 0.85244464 1.48697854
## categories money
## 0.90652971 2.40388211apply(X = boxDat[,1:7], MARGIN = 2, FUN = skewness)## clicks meanTime sessionTime articlesSeen articlesBought
## 2.02035681 -0.33015893 -0.32230048 -0.02762570 1.48697854
## categories money
## -0.09347684 2.40388211qplot(x = sessionTime, geom = "density", data = dat)
qplot(x = sessionTime, geom = "density", data = boxDat)
1.2 Transformaciones sobre múltiples variables
Transformaciones para tratar con outliers: en caso de querer utilizar algún modelo que sea sensible a outliers, se puede utilizar la transformación spatial sign, que proyecta todos los puntos sobre una esfera multidimensional, haciendo que estén a la misma distancia del centro.
Reducción de dimensionalidad:
caretpermite aplicar PCA (Análisis de Componentes Principales) de manera sencilla usandopreProcess. Es recomendable realizar un análisis de la varianza explicada por los componentes antes de preprocesar los datos, aunque si no lo hacemos,caretescoge un umbral del 95% de varianza por defecto.
dat <- as.data.frame(dat)
#Transformación spatial sign ----
pre <- preProcess(x = dat, method = "spatialSign")
pre## Created from 62122 samples and 8 variables
##
## Pre-processing:
## - centered (7)
## - ignored (1)
## - scaled (7)
## - spatial sign transformation (7)spaDat <- predict(object = pre, newdata = dat)
qplot(data = dat, x = clicks, y =categories)
qplot(data = spaDat, x = clicks, y = categories)
#PCA ----
pre <- preProcess(x = dat, method = c("BoxCox", "center", "scale"))
plot(prcomp(x = predict(object = pre,newdata = dat)[,-8]), type = "l")
pre <- preProcess(x = dat, method = c("BoxCox", "pca"), pcaComp = 5)
pre## Created from 62122 samples and 8 variables
##
## Pre-processing:
## - Box-Cox transformation (4)
## - centered (7)
## - ignored (1)
## - principal component signal extraction (7)
## - scaled (7)
##
## Lambda estimates for Box-Cox transformation:
## 0.6, 0.1, 0.3, 0.2
## PCA used 5 components as specifiedpcaDat <- predict(object = pre, newdata = dat)
summary(pcaDat)## userAgent PC1 PC2
## Chrome Generic MacOSX: 7619 Min. :-2.8879 Min. :-3.76398
## Chrome Generic Win10 :24125 1st Qu.:-1.4054 1st Qu.:-0.93745
## Chrome Generic Win7 : 7617 Median :-0.6443 Median : 0.04163
## Firefox Generic Linux: 7581 Mean : 0.0000 Mean : 0.00000
## Firefox Generic Win10:15180 3rd Qu.: 0.9411 3rd Qu.: 0.84063
## Max. : 7.8295 Max. : 7.09301
## PC3 PC4 PC5
## Min. :-6.40904 Min. :-4.3525 Min. :-5.10667
## 1st Qu.:-0.62434 1st Qu.:-0.3090 1st Qu.:-0.20500
## Median :-0.09428 Median : 0.1155 Median : 0.01477
## Mean : 0.00000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.57278 3rd Qu.: 0.4140 3rd Qu.: 0.17686
## Max. : 2.93416 Max. : 3.3992 Max. : 1.654871.3 Tratar con datos incompletos
Entre las distintas opciones disponibles cuando no tenemos todos los datos en un conjunto, preProcess incluye de manera nativa la posibilidad de interpolar los datos faltantes mediante el método k-Nearest Neighbours.
#Creamos valores faltantes en el dataset
missingDat <- dat
samples <- 100
missingDat[sample(x = 1:dim(dat)[1], size = samples, replace = FALSE), 1] <- NaN
#Comparamos la distribución de la variable clicks antes y después de crear nuevos valores
summary(missingDat$clicks)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 16.00 34.00 87.19 84.00 811.00 100pre1 <- preProcess(x = missingDat, method = "knnImpute")
pre1 <- predict(object = pre1, newdata = missingDat)
summary(pre1$clicks)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.71200 -0.58140 -0.43440 0.00006 -0.02605 5.91100pre2 <- preProcess(x = missingDat, method = c("scale", "center"))
pre2 <- predict(object = pre2, newdata = missingDat)
summary(pre2$clicks)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.71200 -0.58140 -0.43440 0.00000 -0.02605 5.91100 1001.4 Eliminar variables
En ocasiones es conveniente o incluso necesario eliminar variables de nuestros datasets. Situaciones como correlaciones altas o variables con muy pocos distintos valores pueden alterar o incluso perjudicar el modelo. caret propone la función nearZeroVar para hallar las variables que presentan problemas de baja varianza.
1.4.1 Encontrar variables correlacionadas entre sí
En el libro de referencia se incluye una discusión muy interesante sobre el análisis de correlaciones entre variables usando corrplot. Si sólo se desea utilizar caret, como es el caso del taller, se puede utilizar directamente la función findCorrelation, que encuentra en una matriz de correlaciones las variables correladas entre sí, aplica un algoritmo de decisión (ver en el libro, página 47) y devuelve un vector numérico de variables a eliminar.
# Variables con poca varianza
nearZeroVar(x = dat)## integer(0)# Variables correladas
findCorrelation(x = cor(x = dat[,-8]), names = TRUE)## [1] "articlesBought" "categories"corrplot(cor(x = dat[,-8]))
2. Separación de datos
Tras realizar el preprocesado deseado, lo usual es separar los datos siguiendo algún tipo de técnica como bootstrapping, validación cruzada, etc. caret proporciona una manera sencilla de realizar tanto bootstrapping (muestreos aleatorios con reemplazamiento del dataset) como particiones en entrenamiento/test o en k-folds (ambos estratificando en función de la variable objetivo).
#Bootstrapping
resamples <- createResample(y = dat$userAgent, times = 5)
summary(resamples)## Length Class Mode
## Resample1 62122 -none- numeric
## Resample2 62122 -none- numeric
## Resample3 62122 -none- numeric
## Resample4 62122 -none- numeric
## Resample5 62122 -none- numeric #Comprobamos cómo se distribuye la variable objetivo
userAgentDist <- t(sapply(X = resamples, FUN = function(X){
res <- dat[X,]
return(as.numeric(table(res$userAgent)*100/dim(res)[1]))
}))
userAgentDist <- rbind(as.numeric(table(dat$userAgent)*100)/dim(dat)[1], userAgentDist)
userAgentDist## [,1] [,2] [,3] [,4] [,5]
## 12.26458 38.83487 12.26136 12.20341 24.43579
## Resample1 12.27906 38.90409 12.29355 12.24365 24.27964
## Resample2 12.32575 38.72863 12.12453 12.38370 24.43740
## Resample3 12.48189 39.02965 11.99253 12.03599 24.45993
## Resample4 12.29355 38.87834 12.30965 12.12292 24.39554
## Resample5 12.16960 38.73829 12.36921 12.40623 24.31667#Entrenamiento/Test
partitions <- createDataPartition(y = dat$userAgent, times = 1, p = 0.8)
summary(partitions)## Length Class Mode
## Resample1 49699 -none- numerichead(partitions$Resample1)## [1] 1 2 4 5 6 7 #El resultado es una lista con tantos vectores numéricos como particiones deseadas.
#Para separar los datos, entonces, solamente tenemos que hacer:
train <- dat[partitions$Resample1,]
test <- dat[-partitions$Resample1,]
table(train$userAgent)*100/dim(train)[1]##
## Chrome Generic MacOSX Chrome Generic Win10 Chrome Generic Win7
## 12.26584 38.83378 12.26182
## Firefox Generic Linux Firefox Generic Win10
## 12.20346 24.43510table(test$userAgent)*100/dim(test)[1]##
## Chrome Generic MacOSX Chrome Generic Win10 Chrome Generic Win7
## 12.25952 38.83925 12.25952
## Firefox Generic Linux Firefox Generic Win10
## 12.20317 24.43854 #Si deseamos hacer más particiones, cambiamos el argumento times al número deseado.
#K-folds
folds <- createFolds(y = dat$userAgent, k = 10, list = TRUE, returnTrain = TRUE)
summary(folds)## Length Class Mode
## Fold01 55911 -none- numeric
## Fold02 55909 -none- numeric
## Fold03 55911 -none- numeric
## Fold04 55910 -none- numeric
## Fold05 55910 -none- numeric
## Fold06 55909 -none- numeric
## Fold07 55909 -none- numeric
## Fold08 55910 -none- numeric
## Fold09 55909 -none- numeric
## Fold10 55910 -none- numerichead(folds$Fold01)## [1] 1 2 3 4 5 6fold1 <- dat[folds$Fold01,]Máxima distancia
Además de en base a una variable, caret también permite muestrear un dataset para añadirlo a otro asegurando que las nuevas muestras son lo más distintas posibles de las antiguas, mediante maxDissim. Vamos a verlo intentando separar a la gente que usa Chrome y la gente que usa Firefox en nuestro dataset (¡ojo! Este método utiliza por defecto la distancia euclidea).
##Tarda en ejecutar, así que dejo el código pero no lo ejecutaré
#Separamos a los usuarios en función del navegador:
#chrome <- filter(dat, userAgent %in% levels(dat$userAgent)[1:3])
#firefox <- filter(dat, userAgent %in% levels(dat$userAgent)[4:5])
#Añadimos 1000 usuarios de firefox a los de chrome
#usuarios <- maxDissim(a = chrome, b = firefox, n = 10)3. Entrenamiento, evaluación y ajuste de modelos
Llegamos a la parte “estrella” de caret. Se ha hecho tan conocido por añadir el equivalente a una API común para una cantidad enorme de modelos. De manera equivalente a preProcess para el preprocesado de datos, en esta fase la función básica es train. Un ejemplo básico, utilizando un modelo lineal para predecir el dinero que se gastan los clientes, sería:
modelo1 <- train(x = select(dat, -c(articlesBought,money)), y = dat$money, method = "lm")
#equivalente a lm(formula = money ~ ., data = select(dat, -articlesBought))... pero no del todo
modelo1## Linear Regression
##
## 62122 samples
## 6 predictor
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 62122, 62122, 62122, 62122, 62122, 62122, ...
## Resampling results:
##
## RMSE Rsquared
## 76.94432 0.2477551
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
## Train control
Como podemos ver, train ha hecho más cosas que entrenar un único modelo: ha entrenado 25 modelos mediante bootstrap y calculado el error cuadrático medio y \(R^2\) de todos. Ahora bien, ¿Y si queremos utilizar otro tipo de esquema de validación? Entra en juego el parámetro trainControl:
#10 folds de validación cruzada
trControl <- trainControl(method = "cv", number = 5)
#trControl <- trainControl(method = "repeatedcv", number = 5, repeats = 5)
modelo2 <- train(x = select(dat, -c(articlesBought,money)), y = dat$money, method = "lm", trControl = trControl)
modelo2## Linear Regression
##
## 62122 samples
## 6 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 49697, 49697, 49698, 49698, 49698
## Resampling results:
##
## RMSE Rsquared
## 76.97616 0.2483774
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
## modelo2$finalModel##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Coefficients:
## (Intercept) clicks
## 9.2387 0.1048
## meanTime sessionTime
## 0.6457 -0.0249
## articlesSeen categories
## -11.2014 16.6261
## userAgentChrome Generic Win10 userAgentChrome Generic Win7
## -21.0283 0.3597
## userAgentFirefox Generic Linux userAgentFirefox Generic Win10
## 1.1083 0.5293Preprocesado
A la hora de interpretar estos coeficientes, nos vendría bien que todas las variables estuvieran en las mismas unidades. Afortunadamente, ¡Sabemos hacer eso con preProcess! En realidad, ni siquiera necesitamos llamar a la función. Es posible incluir todo el preprocesado dentro de la misma llamada a train, lo que nos aligera más aún el código:
#Normalizando las variables y aplicando PCA
trControl <- trainControl(method = "cv", number = 5)
modelo3 <- train(x = select(dat, -c(articlesBought,money)), y = dat$money, method = "lm", trControl = trControl, preProcess = c("center", "scale"))
modelo3## Linear Regression
##
## 62122 samples
## 6 predictor
##
## Pre-processing: centered (5), scaled (5), ignore (1)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 49698, 49698, 49697, 49698, 49697
## Resampling results:
##
## RMSE Rsquared
## 76.97357 0.2484465
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
## modelo3$finalModel##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Coefficients:
## (Intercept) clicks
## 72.5097 12.8356
## meanTime sessionTime
## 26.4841 -17.8771
## articlesSeen categories
## -68.4938 94.3395
## userAgentChrome Generic Win10 userAgentChrome Generic Win7
## -21.0283 0.3597
## userAgentFirefox Generic Linux userAgentFirefox Generic Win10
## 1.1083 0.5293Ajuste de parámetros
Veamos qué sucede si, en vez de regresión lineal, utilizamos un random forest:
trimDat <- dat[createDataPartition(y = dat$money, p = 0.1, list = FALSE),]
trControl <- trainControl(method = "cv", number = 5, verboseIter = TRUE)
modelo4 <- train(x = select(trimDat, -c(articlesBought,money)), y = trimDat$money, method = "ranger", trControl = trControl)## + Fold1: mtry=2
## - Fold1: mtry=2
## + Fold1: mtry=4
## - Fold1: mtry=4
## + Fold1: mtry=6
## - Fold1: mtry=6
## + Fold2: mtry=2
## - Fold2: mtry=2
## + Fold2: mtry=4
## - Fold2: mtry=4
## + Fold2: mtry=6
## - Fold2: mtry=6
## + Fold3: mtry=2
## - Fold3: mtry=2
## + Fold3: mtry=4
## - Fold3: mtry=4
## + Fold3: mtry=6
## - Fold3: mtry=6
## + Fold4: mtry=2
## - Fold4: mtry=2
## + Fold4: mtry=4
## - Fold4: mtry=4
## + Fold4: mtry=6
## - Fold4: mtry=6
## + Fold5: mtry=2
## - Fold5: mtry=2
## + Fold5: mtry=4
## - Fold5: mtry=4
## + Fold5: mtry=6
## - Fold5: mtry=6
## Aggregating results
## Selecting tuning parameters
## Fitting mtry = 2 on full training setmodelo4## Random Forest
##
## 6214 samples
## 6 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 4971, 4972, 4971, 4970, 4972
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared
## 2 46.26980 0.7381076
## 4 46.63874 0.7350057
## 6 47.03515 0.7304645
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 2.Con determinados modelos, train realiza también una evaluación de distintos valores de sus parámetros. Podemos ver qué parámetros tiene cada modelo llamando a modelLookup:
modelLookup("blackboost")## model parameter label forReg forClass probModel
## 1 blackboost mstop #Trees TRUE TRUE TRUE
## 2 blackboost maxdepth Max Tree Depth TRUE TRUE TRUECon esta información, podemos personalizar la búsqueda que se realiza de los parámetros. Vamos a probar con un boosted tree, que tiene 2 parámetros, el número de árboles y la profundidad máxima de cada uno.
busquedaPars <- expand.grid(mstop = c(10, 50, 100, 200))
#modelo5 <- train(x = select(trimDat, -c(articlesBought,money)), y = factor(trimDat$money > 100), method = #"blackboost", trControl = trControl, preProcess = c("center", "scale"), tuneGrid = busquedaPars)Es necesario darle todos los valores. Buscando en la ayuda de la función podemos utilizar los valores por defecto como valores orientativos inciales.
busquedaPars <- expand.grid(mstop = c(10, 50, 100, 200),
maxdepth = c(5,10,20,50))
modelo6 <- train(x = select(trimDat, -c(articlesBought,money,userAgent)), y = trimDat$money, method = "blackboost", trControl = trControl, tuneGrid = busquedaPars)## + Fold1: maxdepth= 5, mstop=200
## - Fold1: maxdepth= 5, mstop=200
## + Fold1: maxdepth=10, mstop=200
## - Fold1: maxdepth=10, mstop=200
## + Fold1: maxdepth=20, mstop=200
## - Fold1: maxdepth=20, mstop=200
## + Fold1: maxdepth=50, mstop=200
## - Fold1: maxdepth=50, mstop=200
## + Fold2: maxdepth= 5, mstop=200
## - Fold2: maxdepth= 5, mstop=200
## + Fold2: maxdepth=10, mstop=200
## - Fold2: maxdepth=10, mstop=200
## + Fold2: maxdepth=20, mstop=200
## - Fold2: maxdepth=20, mstop=200
## + Fold2: maxdepth=50, mstop=200
## - Fold2: maxdepth=50, mstop=200
## + Fold3: maxdepth= 5, mstop=200
## - Fold3: maxdepth= 5, mstop=200
## + Fold3: maxdepth=10, mstop=200
## - Fold3: maxdepth=10, mstop=200
## + Fold3: maxdepth=20, mstop=200
## - Fold3: maxdepth=20, mstop=200
## + Fold3: maxdepth=50, mstop=200
## - Fold3: maxdepth=50, mstop=200
## + Fold4: maxdepth= 5, mstop=200
## - Fold4: maxdepth= 5, mstop=200
## + Fold4: maxdepth=10, mstop=200
## - Fold4: maxdepth=10, mstop=200
## + Fold4: maxdepth=20, mstop=200
## - Fold4: maxdepth=20, mstop=200
## + Fold4: maxdepth=50, mstop=200
## - Fold4: maxdepth=50, mstop=200
## + Fold5: maxdepth= 5, mstop=200
## - Fold5: maxdepth= 5, mstop=200
## + Fold5: maxdepth=10, mstop=200
## - Fold5: maxdepth=10, mstop=200
## + Fold5: maxdepth=20, mstop=200
## - Fold5: maxdepth=20, mstop=200
## + Fold5: maxdepth=50, mstop=200
## - Fold5: maxdepth=50, mstop=200
## Aggregating results
## Selecting tuning parameters
## Fitting mstop = 50, maxdepth = 20 on full training setmodelo6## Boosted Tree
##
## 6214 samples
## 5 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 4971, 4970, 4971, 4972, 4972
## Resampling results across tuning parameters:
##
## maxdepth mstop RMSE Rsquared
## 5 10 55.62235 0.7120240
## 5 50 47.61873 0.7292896
## 5 100 47.61873 0.7292896
## 5 200 47.61873 0.7292896
## 10 10 53.65429 0.7406727
## 10 50 46.37944 0.7436290
## 10 100 46.37944 0.7436290
## 10 200 46.37944 0.7436290
## 20 10 53.65395 0.7406761
## 20 50 46.36326 0.7438808
## 20 100 46.36326 0.7438808
## 20 200 46.36326 0.7438808
## 50 10 53.65395 0.7406761
## 50 50 46.36326 0.7438808
## 50 100 46.36326 0.7438808
## 50 200 46.36326 0.7438808
##
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mstop = 50 and maxdepth = 20.plot(modelo6)
plot(modelo6, metric = "Rsquared")
ggplot(modelo6)## Warning: Ignoring unknown aesthetics: shape
Métricas de evaluación alternativas
Como hemos visto, caret basa sus cálculos en el RMSE para regresión (y en la precisión (accuracy) para clasificación). Por supuesto, estas métricas no serán adecuadas en todos los casos. Una vez más, caret nos permite personalizar qué métrica se utiliza. Vamos a verlo en un ejemplo de clasificación utilizando el área bajo la curva ROC.
trControl <- trainControl(method = "cv", number = 5, summaryFunction = twoClassSummary, classProbs = TRUE, verboseIter = TRUE)
busquedaPars <- expand.grid(mtry = c(2,4,6))
dat$vip <- factor(dat$money > 100)
levels(dat$vip) <- c("no", "yes")
vipDatTrain <- createDataPartition(y = dat$vip, times = 1, p = 0.1)$Resample1
vipDat <- dat[vipDatTrain,]
modelo7 <- train(x = select(vipDat, -c(articlesBought,money, vip)), y = vipDat$vip, method = "ranger", trControl = trControl, metric = "ROC", tuneGrid = busquedaPars)## + Fold1: mtry=2
## - Fold1: mtry=2
## + Fold1: mtry=4
## - Fold1: mtry=4
## + Fold1: mtry=6
## - Fold1: mtry=6
## + Fold2: mtry=2
## - Fold2: mtry=2
## + Fold2: mtry=4
## - Fold2: mtry=4
## + Fold2: mtry=6
## - Fold2: mtry=6
## + Fold3: mtry=2
## - Fold3: mtry=2
## + Fold3: mtry=4
## - Fold3: mtry=4
## + Fold3: mtry=6
## - Fold3: mtry=6
## + Fold4: mtry=2
## - Fold4: mtry=2
## + Fold4: mtry=4
## - Fold4: mtry=4
## + Fold4: mtry=6
## - Fold4: mtry=6
## + Fold5: mtry=2
## - Fold5: mtry=2
## + Fold5: mtry=4
## - Fold5: mtry=4
## + Fold5: mtry=6
## - Fold5: mtry=6
## Aggregating results
## Selecting tuning parameters
## Fitting mtry = 2 on full training setmodelo7## Random Forest
##
## 6213 samples
## 6 predictor
## 2 classes: 'no', 'yes'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 4970, 4970, 4970, 4971, 4971
## Resampling results across tuning parameters:
##
## mtry ROC Sens Spec
## 2 0.9448882 0.9187287 0.7289832
## 4 0.9437988 0.9172736 0.7289705
## 6 0.9428121 0.9183140 0.7339705
##
## ROC was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.plot(modelo7)
Predecir nuevos datos
Tras crear el modelo, es usual que comprobemos su rendimiento con datos que no hayan sido vistos durante el entrenamiento. caret sigue la sintaxis básica de R: la función predict.
vipTest <- select(dat[-vipDatTrain,], -c(articlesBought, money))
predict(object = modelo7, newdata = head(select(vipTest, -vip)), type = "raw")## [1] no no no no no no
## Levels: no yespredict(object = modelo7, newdata = head(select(vipTest, -vip)), type = "prob")## no yes
## 1 1.0000 0.0000
## 2 0.9976 0.0024
## 3 1.0000 0.0000
## 4 1.0000 0.0000
## 5 1.0000 0.0000
## 6 0.9970 0.0030confusionMatrix(data = table(predict(object = modelo7, newdata = select(vipTest, -vip), type = "raw"),vipTest$vip))## Confusion Matrix and Statistics
##
##
## no yes
## no 39605 3211
## yes 3692 9401
##
## Accuracy : 0.8765
## 95% CI : (0.8738, 0.8792)
## No Information Rate : 0.7744
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6513
## Mcnemar's Test P-Value : 7.592e-09
##
## Sensitivity : 0.9147
## Specificity : 0.7454
## Pos Pred Value : 0.9250
## Neg Pred Value : 0.7180
## Prevalence : 0.7744
## Detection Rate : 0.7084
## Detection Prevalence : 0.7658
## Balanced Accuracy : 0.8301
##
## 'Positive' Class : no
## Comparar modelos
Vista la facilidad para entrenar con una gran cantidad de modelos usando train, es normal que queramos comparar rápidamente el rendimiento de cada uno. Lo haremos con la función resamples, siempre que hayan sido entrenando siguiendo el mismo esquema de validación:
modelos <- resamples(list(lm = modelo2, lm2 = modelo3, rf = modelo4))
summary(modelos)##
## Call:
## summary.resamples(object = modelos)
##
## Models: lm, lm2, rf
## Number of resamples: 5
##
## RMSE
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## lm 75.69 76.52 77.27 76.98 77.58 77.82 0
## lm2 74.95 76.16 77.59 76.97 77.79 78.38 0
## rf 43.76 44.09 45.12 46.27 47.34 51.04 0
##
## Rsquared
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## lm 0.2429 0.2454 0.2474 0.2484 0.2486 0.2576 0
## lm2 0.2401 0.2455 0.2488 0.2484 0.2535 0.2543 0
## rf 0.7104 0.7240 0.7494 0.7381 0.7497 0.7570 0bwplot(modelos, metric = "RMSE")
Segunda parte
Pasamos ahora a la segunda parte del taller. Hemos visto cómo caret nos ayuda a preprocesar, separar, entrenar, evaluar, ajustar y comparar modelos. Ese es el núcleo del paquete y, con diferencia, aquello que mejor hace y por lo que es conocido e importante. En esta segunda parte estudiaremos algunas pequeñas funcionalidades extra que el paquete incluye que también son interesantes. Por supuesto, lo que podamos ver es sólo una pequeña parte de todo lo que ofrece. Tanto en el libro como en la web referenciados al principio se pueden encontrar estos ejemplos y mucho más material.
4. Procesamiento en paralelo
Uno de los problemas de R es el trabajo con grandes cantidades de datos. La capacidad de cómputo que permite la distribución de tareas en herramientas como Hadoop o Spark empequeñece al código de R más eficiente. Hay paquetes en R que permiten aprovechar los núcleos de nuestros procesadores para paralelizar ciertas tareas, como doMC en Mac/Linux o doParallel en Windows. De manera análoga a los modelos de Machine Learning, caret nos ofrece una interfaz simple a la paralelización:
trimDat <- dat[createDataPartition(y = dat$money, p = 0.02, list = FALSE),]
trControl <- trainControl(method = "cv", number = 5, verboseIter = TRUE)
system.time({modeloSimple <- train(x = select(trimDat, -c(articlesBought,money)), y = trimDat$money, method = "rf", trControl = trControl)})## Loading required package: randomForest## randomForest 4.6-12## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ranger':
##
## importance## The following object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:ggplot2':
##
## margin## + Fold1: mtry=2
## - Fold1: mtry=2
## + Fold1: mtry=4
## - Fold1: mtry=4
## + Fold1: mtry=7
## - Fold1: mtry=7
## + Fold2: mtry=2
## - Fold2: mtry=2
## + Fold2: mtry=4
## - Fold2: mtry=4
## + Fold2: mtry=7
## - Fold2: mtry=7
## + Fold3: mtry=2
## - Fold3: mtry=2
## + Fold3: mtry=4
## - Fold3: mtry=4
## + Fold3: mtry=7
## - Fold3: mtry=7
## + Fold4: mtry=2
## - Fold4: mtry=2
## + Fold4: mtry=4
## - Fold4: mtry=4
## + Fold4: mtry=7
## - Fold4: mtry=7
## + Fold5: mtry=2
## - Fold5: mtry=2
## + Fold5: mtry=4
## - Fold5: mtry=4
## + Fold5: mtry=7
## - Fold5: mtry=7
## Aggregating results
## Selecting tuning parameters
## Fitting mtry = 2 on full training set## user system elapsed
## 25.590 0.775 26.613library(doMC)## Loading required package: foreach## Loading required package: iteratorsregisterDoMC(cores = 4)
system.time({modeloParalelo <- train(x = select(trimDat, -c(articlesBought,money)), y = trimDat$money, method = "rf", trControl = trControl)})## Aggregating results
## Selecting tuning parameters
## Fitting mtry = 2 on full training set## user system elapsed
## 23.741 1.229 11.629En caso de utilizar Windows, este enlace da una solución similar a la planteada aquí.
5. Submuestreo para compensar desequilibrios entre clases
Cuando realizamos clasificación binaria, es muy común encontrar datasets muy desequilibrados. Es decir, aquellos en los que hay muchas más muestras de una clase que de otra (normalmente, muchas más de la clase que no nos interesa predecir). Esto puede afectar al rendimiento de algunos modelos al entrenar. Una manera sencilla de atacar este problema es alterando la proporción entre clases, ya sea disminuyendo la clase más abundante, aumentando las muestras de la clase más escasa o mediante técnicas de muestreo híbrido. caret permite utilizar todas de la misma manera. Hay dos opciones. La primera es hacerlo antes de entrenar:
dat$strange <- factor(dat$clicks > 300 & dat$sessionTime < 60)
levels(dat$strange) <- c("notStrange", "strange")
summary(dat$strange)## notStrange strange
## 61742 380summary(downSample(x = dat, y = dat$strange)$strange)## notStrange strange
## 380 380summary(upSample(x = dat, y = dat$strange)$strange)## notStrange strange
## 61742 61742summary(SMOTE(form = strange ~ ., data = dat, perc.over = 500, k = 10)$strange)## notStrange strange
## 3800 2280summary(ROSE(formula = strange ~., data = dat)$data$strange)## notStrange strange
## 30875 31247La segunda opción (y, a veces, más conveniente) es incluir el muestreo dentro del bucle de validación cruzada. Lo único que tenemos que cambiar es la opción sampling dentro de nuestro objeto trainControl.
trControl <- trainControl(method = "cv", number = 5, verboseIter = TRUE)
TrainStrange <- function(){
return(train(x = select(dat, -c(clicks,sessionTime, strange)), y = dat$strange, method = "ranger", trControl = trControl, metric = "Kappa"))
}
#normal <- TrainStrange()
#trControl$sampling <- "down"
#downSampled <- TrainStrange()
#trControl$sampling <- "smote"
#smoted <- PredictStrange()
#trControl$sampling <- "up"
#upSampled <- PredictStrange()
#trControl$sampling <- "rose"
#rosed <- PredictStrange()
load("samplingModels.Rdata")
comparison <- resamples(x = samplingModels)
summary(comparison, metric = "Kappa")##
## Call:
## summary.resamples(object = comparison, metric = "Kappa")
##
## Models: normal, downsampling, upsampling, smote, rose
## Number of resamples: 5
##
## Kappa
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## normal 0.7553 0.7701 0.7850 0.7839 0.7988 0.8101 0
## downsampling 0.3995 0.4598 0.4698 0.4656 0.4792 0.5198 0
## upsampling 0.7668 0.7778 0.7834 0.7876 0.8012 0.8089 0
## smote 0.5487 0.5558 0.5879 0.5978 0.6325 0.6640 0
## rose 0.4062 0.4190 0.4301 0.4327 0.4485 0.4598 0bwplot(comparison, metric = "Kappa")
6. Importancia de las variables
Además de los métodos propios de cada modelo, caret añade la posibilidad de calcular una métrica de importancia de las variables del modelo. Para clasificación, construye el área bajo la curva ROC de cada variable y para regresión, ajusta un modelo (a elegir entre un modelo lineal o un loess smoother) y devuelve la \(R^2\).
#Importancia en clasificación
trimDat <- dat[createDataPartition(y = dat$money, p = 0.1, list = FALSE),]
trControl <- trainControl(method = "cv", number = 5, verboseIter = TRUE)
modeloRF <- train(x = select(trimDat, -c(articlesBought,money)), y = factor(trimDat$money > 0), method = "rf", trControl = trControl, importance = TRUE)## Aggregating results
## Selecting tuning parameters
## Fitting mtry = 5 on full training set #Importancia del modelo
varImpPlot(x =modeloRF$finalModel, scale = FALSE)
#Importancia de caret
varImp(object = modeloRF, useModel = FALSE) ## ROC curve variable importance
##
## Importance
## clicks 100.00
## vip 73.93
## strange 49.65
## userAgent 46.80
## meanTime 43.36
## categories 38.59
## articlesSeen 28.99
## sessionTime 0.00 varImp(object = modeloRF, useModel = FALSE, scale = FALSE)## ROC curve variable importance
##
## Importance
## clicks 0.8340
## vip 0.6635
## strange 0.5047
## userAgent 0.4861
## meanTime 0.4636
## categories 0.4324
## articlesSeen 0.3696
## sessionTime 0.1800 plot(varImp(object = modeloRF, useModel = FALSE, scale = FALSE))
#Importancia en regresión
modeloLM <- train(x = select(dat, -c(articlesBought,money)), y = dat$money, method = "lm", trControl = trControl)## Aggregating results
## Fitting final model on full training set #Importancia del modelo
varImp(object = modeloLM, useModel = TRUE, scale = FALSE)## lm variable importance
##
## Overall
## vipyes 260.178
## categories 40.442
## articlesSeen 38.245
## sessionTime 36.175
## meanTime 31.273
## userAgentChrome Generic Win10 9.941
## strangestrange 3.498
## userAgentChrome Generic Win7 2.588
## userAgentFirefox Generic Linux 2.103
## userAgentFirefox Generic Win10 1.256
## clicks 0.215 summary(modeloLM$finalModel)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -127.38 -24.61 -7.84 18.08 655.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.825e+01 1.056e+00 17.272 < 2e-16 ***
## clicks -6.761e-04 3.144e-03 -0.215 0.82977
## meanTime 2.524e-01 8.072e-03 31.273 < 2e-16 ***
## sessionTime -1.349e-02 3.730e-04 -36.175 < 2e-16 ***
## articlesSeen -5.757e+00 1.505e-01 -38.245 < 2e-16 ***
## categories 7.438e+00 1.839e-01 40.442 < 2e-16 ***
## userAgentChrome Generic Win10 -7.438e+00 7.483e-01 -9.941 < 2e-16 ***
## userAgentChrome Generic Win7 2.232e+00 8.625e-01 2.588 0.00965 **
## userAgentFirefox Generic Linux 1.816e+00 8.635e-01 2.103 0.03549 *
## userAgentFirefox Generic Win10 9.384e-01 7.474e-01 1.256 0.20927
## vipyes 1.528e+02 5.874e-01 260.178 < 2e-16 ***
## strangestrange -9.765e+00 2.792e+00 -3.498 0.00047 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53.23 on 62110 degrees of freedom
## Multiple R-squared: 0.6407, Adjusted R-squared: 0.6406
## F-statistic: 1.007e+04 on 11 and 62110 DF, p-value: < 2.2e-16 #Importancia de caret (tarda :( )
#varImp(object = modeloLM, useModel = FALSE, scale = FALSE)7. Selección de variables mediante algoritmo genético
Además de los métodos más conocidos (stepwise, forward, backward) para seleccionar variables, caret ha incorporado varios métodos más elaborados de cara a esta tarea. Aquí vamos a ver, como pequeño ejemplo, el uso de un algoritmo genético para esta tarea. Lo utilizaremos con el modelo lineal que intenta predecir el dinero gastado:
#Desgraciadamente, tarda demasiado.
#ctrl <- gafsControl(functions = caretGA)
#modeloLMG <- gafs(x = select(dat, -c(articlesBought,money)),
#y = dat$money,iters = 100, gafsControl = ctrl, method = "lm")8. Optimización de hiperparámetros
Por´último, vamos a ver dos alternativas a la búsqueda de rejilla a la hora de optimizar hiperparámetros de un modelo: la búsqueda aleatoria y lo que caret llama “adaptive resampling”. La búsqueda aleatoria es la más sencilla conceptualmente: únicamente tenemos que definir cuántas combinaciones (aleatorias) de parámetros queremos que se prueben y lanzarlo.
#trimDat <- dat[createDataPartition(y = dat$money, p = 0.7, list = FALSE),]
#trControl <- trainControl(method = "cv", number = 5, verboseIter = TRUE, search = "random")
#system.time({modeloRandom <- train(x = select(trimDat, -c(articlesBought,money, userAgent)), y = trimDat$money, method = "svmRadial", trControl = trControl, tuneLength = 15)})
#modeloRandomLa segunda opción, “adaptive resampling”, adapta los valores que va probando en la dirección de los valores óptimos (siguiendo lo dicho en este paper). Para hacerlo en R:
#trControl <- trainControl(method = "cv", number = 5, verboseIter = TRUE,
# adaptive = list(min = 5, alpha = 0.05, method = "gls", complete = TRUE))
#system.time({modeloAdaptive <- train(x = select(trimDat, -c(articlesBought,money, userAgent)), y = trimDat$money, method = "svmRadial", trControl = trControl, tuneLength = 15)})
#modeloAdaptive