Universidad de Deusto - Facultad de Ingeniería

1 Bank Marketing dataset

Este dataset ha sido descargado del repositorio de ML del Center for Machine Learning and Intelligent Systems de la University of California, Irvine: Bank Marketing Data Set

1.1 Información del dataset

Los datos incluidos en este dataset pertenecen a campañas de marketing directo de una instutución bancaria portuguesa. Estas campañas de marketing se hicieron en base a llamadas telefónicas que, con frecuencia, necesitaron ser más de una para conocer si el producto (un depósito a plazo) fue finalmente suscrito o no.

Objetivo: mediante un modelo de clasificación, predecir si el cliente suscribirá un depósito a plazo o no.
Source: [Moro et al., 2014] S. Moro, P. Cortez and P. Rita. A datos_imputados-Driven Approach to Predict the Success of Bank Telemarketing. Decision Support Systems, Elsevier, 62:22-31, June 2014.

1.2 Atributos/Variables

Cliente
- age [numérica]: edad del cliente
- job [categórica]: tipo de trabajo
  ‘admin.’, ‘blue-collar’, ‘entrepreneur’, ‘housemaid’, ‘management’, ‘retired’, ‘self-employed’, ‘services’, ‘student’, ‘technician’, ‘unemployed’, ‘unknown’
- marital [categórica]: estado civil
  ‘divorced’, ‘married’, ‘single’, ‘unknown’; Note: ‘divorced’ means divorced or widowed
- education [categorical]: nivel de estudios
  ‘basic.4y’, ‘basic.6y’, ‘basic.9y’, ‘high.school’, ‘illiterate’, ‘professional.course’, ‘university.degree’, ‘unknown’
- default [categórica]: ¿tiene crédito pendiente de pago?
  ‘no’, ‘yes’, ‘unknown’
- housing [categórica]: ¿tiene un préstamos hipotecario?
  ‘no’, ‘yes’, ‘unknown’
- loan [categórica]: ¿tiene un préstamo personal? ‘no’, ‘yes’, ‘unknown’
Último contacto de la campaña actual
- contact [binaria]: tipo de comunicaión mediante la que se ha contactado
  ‘cellular’, ‘telephone’
- month [categórica]: mes del año del último contacto
  categorical: ‘jan’, ‘feb’, ‘mar’, …, ‘nov’, ‘dec’
- day_of_week [categórica]: día de la semana del último contacto
  ‘mon’, ‘tue’, ‘wed’, ‘thu’, ‘fri’)
- duration [numérica]: duración de la última comunicación, en segundos
  Importante: este atributo tiene una fuerte relación con la variable de salida (e.g., si duration=0 entonces y=‘no’). La duración no se puede conocer antes de realizar la llamada y, además, tras finalizar ésta se conocerá el valor de y. Por tanto, este atributo debería ser incluído sólo como referencia y eliminado si la intención es obtener un modelo predictivo realista.
Otros atributos
- campaign [numérica]: cuántas veces que se ha contactado con este cliente durante la campaña actual (incluída la última llamada)
- pdays [numérica]: días que han pasado desde la última vez que se contactó al cliente por una campaña anterior
  Importante: 999 significa que el cliente no ha sido contactado con anterioridad (ésta es la primera vez que se contacta al cliente (para una campaña publicitaria)
- previous [numérica]: cuántas veces que se ha contactado con este cliente con anterioridad a esta campaña por cualquier motivo
- poutcome [categórica]: resultado de la anterior campaña de marketing
  ‘failure’, ‘nonexistent’, ‘success’)
Contexto socioeconómico
- emp.var.rate [numérica]: tasa de variación de empleo (indicador trimestral)
- cons.price.idx [numérica]: ??ndice de precios al consumidor (indicador mensual)
- cons.conf.idx [numérica]: ??ndice de confianza del consumidor (indicador mensual)
- euribor3m [numérica]: euribor a 3 meses (indicador diario)
- nr.employed [numérica]: número de empleados (indicador trimestral)
Variable de salida (target deseados)
- y [binaria]: ¿ha suscrito el cliente el depósito a plazo ofertado?
  ‘yes’, ‘no’

2 Exploración del dataset

2.1 Preparar el entorno de trabajo

Borrar los objetos del workspace y establecer el working directory en el que se encuentran los datos. Además cargar los datos pre-procesados (para no hacerlo cada vez que ejecuto el informe):

# Limpiar el global environment
rm(list = ls())

# Configurar el working directory
setwd("./")

Leer el dataset a utilizar, en este caso un CSV con separador “;” y decimales formateados con “.”:

# Leer datos de origen
rawdata <- read.csv2("./data/bank-additional/bank-additional-full.csv",
                     sep = ";", dec = ".")
data <- rawdata # Copia de los datos - "unknown" como clase

# Definir valor semilla para todo el proceso
seed.val <- 1492

2.2 Estructura del dataset

dim(datos_clase): el dataset está compuesto de 41188 observaciones (clientes) caracterizadas por 21 atributos cada una (20 de entrada y 1 de salida).

Para comprobar que la lectura de los atributos ha funcionado correctamente y que coincide con la lista descrita en la sección “Atributos/Variables”, se puede recurrir tanto a la estructura del dataset vía str(datos_imputados) como a ver un resumen del mismo vía summary(datos_imputados). En esta ocasión utilizaré ambas opciones, para poder obervar mejor tanto el tipo de dato de los diferentes atributos como algunas de sus características:

# Estructura del objeto
str(data)

| 'data.frame': 41188 obs. of  21 variables:
|  $ age           : int  56 57 37 40 56 45 59 41 24 25 ...
|  $ job           : Factor w/ 12 levels "admin.","blue-collar",..: 4 8 8 1 8 8 1 2 10 8 ...
|  $ marital       : Factor w/ 4 levels "divorced","married",..: 2 2 2 2 2 2 2 2 3 3 ...
|  $ education     : Factor w/ 8 levels "basic.4y","basic.6y",..: 1 4 4 2 4 3 6 8 6 4 ...
|  $ default       : Factor w/ 3 levels "no","unknown",..: 1 2 1 1 1 2 1 2 1 1 ...
|  $ housing       : Factor w/ 3 levels "no","unknown",..: 1 1 3 1 1 1 1 1 3 3 ...
|  $ loan          : Factor w/ 3 levels "no","unknown",..: 1 1 1 1 3 1 1 1 1 1 ...
|  $ contact       : Factor w/ 2 levels "cellular","telephone": 2 2 2 2 2 2 2 2 2 2 ...
|  $ month         : Factor w/ 10 levels "apr","aug","dec",..: 7 7 7 7 7 7 7 7 7 7 ...
|  $ day_of_week   : Factor w/ 5 levels "fri","mon","thu",..: 2 2 2 2 2 2 2 2 2 2 ...
|  $ duration      : int  261 149 226 151 307 198 139 217 380 50 ...
|  $ campaign      : int  1 1 1 1 1 1 1 1 1 1 ...
|  $ pdays         : int  999 999 999 999 999 999 999 999 999 999 ...
|  $ previous      : int  0 0 0 0 0 0 0 0 0 0 ...
|  $ poutcome      : Factor w/ 3 levels "failure","nonexistent",..: 2 2 2 2 2 2 2 2 2 2 ...
|  $ emp.var.rate  : num  1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 ...
|  $ cons.price.idx: num  94 94 94 94 94 ...
|  $ cons.conf.idx : num  -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 ...
|  $ euribor3m     : num  4.86 4.86 4.86 4.86 4.86 ...
|  $ nr.employed   : num  5191 5191 5191 5191 5191 ...
|  $ y             : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...

# Resumen del objeto
summary(data)

|       age                 job            marital     
|  Min.   :17.00   admin.     :10422   divorced: 4612  
|  1st Qu.:32.00   blue-collar: 9254   married :24928  
|  Median :38.00   technician : 6743   single  :11568  
|  Mean   :40.02   services   : 3969   unknown :   80  
|  3rd Qu.:47.00   management : 2924                   
|  Max.   :98.00   retired    : 1720                   
|                  (Other)    : 6156                   
|                education        default         housing     
|  university.degree  :12168   no     :32588   no     :18622  
|  high.school        : 9515   unknown: 8597   unknown:  990  
|  basic.9y           : 6045   yes    :    3   yes    :21576  
|  professional.course: 5243                                  
|  basic.4y           : 4176                                  
|  basic.6y           : 2292                                  
|  (Other)            : 1749                                  
|       loan            contact          month       day_of_week
|  no     :33950   cellular :26144   may    :13769   fri:7827   
|  unknown:  990   telephone:15044   jul    : 7174   mon:8514   
|  yes    : 6248                     aug    : 6178   thu:8623   
|                                    jun    : 5318   tue:8090   
|                                    nov    : 4101   wed:8134   
|                                    apr    : 2632              
|                                    (Other): 2016              
|     duration         campaign          pdays          previous    
|  Min.   :   0.0   Min.   : 1.000   Min.   :  0.0   Min.   :0.000  
|  1st Qu.: 102.0   1st Qu.: 1.000   1st Qu.:999.0   1st Qu.:0.000  
|  Median : 180.0   Median : 2.000   Median :999.0   Median :0.000  
|  Mean   : 258.3   Mean   : 2.568   Mean   :962.5   Mean   :0.173  
|  3rd Qu.: 319.0   3rd Qu.: 3.000   3rd Qu.:999.0   3rd Qu.:0.000  
|  Max.   :4918.0   Max.   :56.000   Max.   :999.0   Max.   :7.000  
|                                                                   
|         poutcome      emp.var.rate      cons.price.idx  cons.conf.idx  
|  failure    : 4252   Min.   :-3.40000   Min.   :92.20   Min.   :-50.8  
|  nonexistent:35563   1st Qu.:-1.80000   1st Qu.:93.08   1st Qu.:-42.7  
|  success    : 1373   Median : 1.10000   Median :93.75   Median :-41.8  
|                      Mean   : 0.08189   Mean   :93.58   Mean   :-40.5  
|                      3rd Qu.: 1.40000   3rd Qu.:93.99   3rd Qu.:-36.4  
|                      Max.   : 1.40000   Max.   :94.77   Max.   :-26.9  
|                                                                        
|    euribor3m      nr.employed     y        
|  Min.   :0.634   Min.   :4964   no :36548  
|  1st Qu.:1.344   1st Qu.:5099   yes: 4640  
|  Median :4.857   Median :5191              
|  Mean   :3.621   Mean   :5167              
|  3rd Qu.:4.961   3rd Qu.:5228              
|  Max.   :5.045   Max.   :5228              
|

Cabe destacar lo siguiente:

Imbalanced data con un ratio ≈8:1 en la variable de salida y. Abordaré este tema en el apartado Addressing Imbalanced data.
month tiene 10 niveles (‘mar’, …‘dec’) en lugar de 12 como cabría esperar.
nr.employed es de tipo numeric. Parece ser el Portugal labour force expresado en miles.
pdays es de tipo numeric, pero el valor 999 es como una categoría, puesto que significa que no se ha contactado nunca con el cliente en cuestión en ninguna campaña anterior.
En el dataset no hay NAs (missing values), sino que estos fueron etiquetados como “unknown” en algunas variables categóricas. Ante una situación como esta, ¿qué se puede hacer?:
1. eliminar la categoría.
2. tratar la categoría como una más.
3. imputar los missing values (sustituir valores no informados en una observación -NA- por otros en función de la distribución resto de categorías disponibles para un atributo concreto).

2.3 Preprocesado

2.3.1 Factorizar y eliminar variables

# Antes del cambio
table(data$pdays)

| 
|     0     1     2     3     4     5     6     7     8     9    10    11 
|    15    26    61   439   118    46   412    60    18    64    52    28 
|    12    13    14    15    16    17    18    19    20    21    22    25 
|    58    36    20    24    11     8     7     3     1     2     3     1 
|    26    27   999 
|     1     1 39673

# Renombrar "pdays" como "prev.campaign" y factorizarla
names(data)[13] <- "prev.campaign"
data$prev.campaign <- cut(x = data$prev.campaign,
                          breaks = c(seq(0,7*4+1,7),999+1),
                          right = FALSE,
                          labels = c("Week 1","Week 2",
                                     "Week 3","Week 4","Never"))
# Después del cambio
table(data$prev.campaign)

| 
| Week 1 Week 2 Week 3 Week 4  Never 
|   1117    316     74      8  39673

# Eliminar la variable "duration" (en este caso por exceso de información)
data$duration <- NULL

2.3.2 Tratamiento de NAs

En este dataset todas las observaciones están completas (sum(complete.cases(datos_orig))), por lo que no hay NAs como tal. Sin embargo, hay algunas variables categóricas con un nivel denominado “unknown” que hace las veces de missing value (NA).

¿Cuántos NA hay por en cada varible?:

# Recuento de NAs por variable
sapply(X = data, FUN = function(x) sum(length(which(is.na(x)))))

|            age            job        marital      education        default 
|              0              0              0              0              0 
|        housing           loan        contact          month    day_of_week 
|              0              0              0              0              0 
|       campaign  prev.campaign       previous       poutcome   emp.var.rate 
|              0              0              0              0              0 
| cons.price.idx  cons.conf.idx      euribor3m    nr.employed              y 
|              0              0              0              0              0

# Recuento de "unknown" por variable
sapply(X = data, FUN = function(x) sum(x == "unknown"))

|            age            job        marital      education        default 
|              0            330             80           1731           8597 
|        housing           loan        contact          month    day_of_week 
|            990            990              0              0              0 
|       campaign  prev.campaign       previous       poutcome   emp.var.rate 
|              0              0              0              0              0 
| cons.price.idx  cons.conf.idx      euribor3m    nr.employed              y 
|              0              0              0              0              0

La variable default parece tener demasiadas observaciones del nivel “unknown” (observaciones no informadas), a diferencia de job, marital, education, housing y loan. Veamos qué porcentaje del total de observaciones supone:

# Recuento de "unknown" en la variable "default"
round(prop.table(table(data$default=="unknown")), digits = 3)

| 
| FALSE  TRUE 
| 0.791 0.209

# Demasiados "unknown"; mejor descartar esta variable
data$default <- NULL

La proporción es tal (20.9 %), que decido no utilizar la variable default para la clasificación.

Para el resto de atributos con valores “unknown”, decido hacer una imputación utilizando el paquete MICE (Multivariate Imputation by Chained Equations):

if (!require("mice")) {
  install.packages("mice",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(mice)
data.na <- data

# Convertir los valores "unknown" en NAs
data.na[data.na == "unknown"] = NA

El patrón que siguen los NAs se puede observar en forma de tabla con el comando md.pattern(datos_na). Esta tabla no es muy intuitiva y analizarla requiere tiempo, por lo que voy a intentar recoger la misma información de forma visual:

if (!require("VIM")) {
  install.packages("VIM",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(VIM)

# Distribución de missing values
distro.na <- aggr(data.na, col = c("lightblue", "darkorange"),
     combined = TRUE,
     numbers = TRUE,
     prop = FALSE,
     sortVars = TRUE,
     labels = names(data.na),
     cex.axis = .7,
     cex.numbers = .7,
     only.miss = TRUE,
     ylab = "Patrón de NAs")

| 
|  Variables sorted by number of missings: 
|        Variable Count
|       education  1731
|         housing   990
|            loan   990
|             job   330
|         marital    80
|             age     0
|         contact     0
|           month     0
|     day_of_week     0
|        campaign     0
|   prev.campaign     0
|        previous     0
|        poutcome     0
|    emp.var.rate     0
|  cons.price.idx     0
|   cons.conf.idx     0
|       euribor3m     0
|     nr.employed     0
|               y     0

Una vez hecho esto, vamos a imputar los missing values de cada atributo dándoles como valor alguno de los posibles niveles de tal categoría. Para ello, utilizaré el método “pmm” (Predictive Mean Matching). Otros métodos que funcionan con variables categóricas son “rf” (Random Forest) y “cart” (Classification and Regression Trees) (ver methods(mice)):

# Imputar los datos; 5 posibles imputaciones por NA y 5 iteraciones
data.temp <- mice(data.na, m = 5, maxit = 5,
                  method = "pmm", seed = seed.val)

| 
|  iter imp variable
|   1   1  job  marital  education  housing  loan
|   1   2  job  marital  education  housing  loan
|   1   3  job  marital  education  housing  loan
|   1   4  job  marital  education  housing  loan
|   1   5  job  marital  education  housing  loan
|   2   1  job  marital  education  housing  loan
|   2   2  job  marital  education  housing  loan
|   2   3  job  marital  education  housing  loan
|   2   4  job  marital  education  housing  loan
|   2   5  job  marital  education  housing  loan
|   3   1  job  marital  education  housing  loan
|   3   2  job  marital  education  housing  loan
|   3   3  job  marital  education  housing  loan
|   3   4  job  marital  education  housing  loan
|   3   5  job  marital  education  housing  loan
|   4   1  job  marital  education  housing  loan
|   4   2  job  marital  education  housing  loan
|   4   3  job  marital  education  housing  loan
|   4   4  job  marital  education  housing  loan
|   4   5  job  marital  education  housing  loan
|   5   1  job  marital  education  housing  loan
|   5   2  job  marital  education  housing  loan
|   5   3  job  marital  education  housing  loan
|   5   4  job  marital  education  housing  loan
|   5   5  job  marital  education  housing  loan

# Comprobar las primeras observaciones con datos imputados
# (por ejemplo para la variable "education")
head(data.temp$imp$education)

# Reemplazar NAs con los valores imputados en el 1º de los 5 dataset (por ejemplo)
data.imp <- complete(data.temp, 1)

# Eliminar las categorías "unknown" (de las que ahora hay 0 observaciones)
if (!require("gdata")) {
  install.packages("gdata",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(gdata)
data.imp <- drop.levels(data.imp)

# Borrar cálculos intermedios
rm(data, data.na, data.temp, distro.na)

2.3.3 Normalización y binarización

El dataset con las variables numéricas normalizadas y las variables categóricas binarizadas se crea como se indica a continuación:

# Normalizar variables numéricas
if (!require("clusterSim")) {
  install.packages("clusterSim",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(clusterSim)
data.norm <- data.Normalization(data.imp[sapply(data.imp, is.numeric)],
                                type = "n4", # (x-min)/range
                                normalization = "column")

# Crear matriz de las variables categóricas
form <- as.formula(
    paste0('~', paste(names(data.imp[sapply(data.imp,is.factor)][!names(
        data.imp[sapply(data.imp, is.factor)]) %in% 'y']),collapse = ' + ')))
matrix.cat <- model.matrix(form, data = data.imp)

# Crear datasets finales
y <- data.imp$y
data.norm <- cbind(data.norm, matrix.cat, y)

# Borrar cálculos intermedios
rm(form, y, matrix.cat)

De esta forma, los datos con los que trabajar quedan as?? (se muestran únicamente los primeros 15 atributos):

# Estructura del objeto (mostrar sólo los 15 primeros atributos)
str(data.norm, list.len=15)

| 'data.frame': 41188 obs. of  50 variables:
|  $ age                         : num  0.481 0.494 0.247 0.284 0.481 ...
|  $ campaign                    : num  0 0 0 0 0 0 0 0 0 0 ...
|  $ previous                    : num  0 0 0 0 0 0 0 0 0 0 ...
|  $ emp.var.rate                : num  0.938 0.938 0.938 0.938 0.938 ...
|  $ cons.price.idx              : num  0.699 0.699 0.699 0.699 0.699 ...
|  $ cons.conf.idx               : num  0.603 0.603 0.603 0.603 0.603 ...
|  $ euribor3m                   : num  0.957 0.957 0.957 0.957 0.957 ...
|  $ nr.employed                 : num  0.86 0.86 0.86 0.86 0.86 ...
|  $ (Intercept)                 : num  1 1 1 1 1 1 1 1 1 1 ...
|  $ jobblue-collar              : num  0 0 0 0 0 0 0 1 0 0 ...
|  $ jobentrepreneur             : num  0 0 0 0 0 0 0 0 0 0 ...
|  $ jobhousemaid                : num  1 0 0 0 0 0 0 0 0 0 ...
|  $ jobmanagement               : num  0 0 0 0 0 0 0 0 0 0 ...
|  $ jobretired                  : num  0 0 0 0 0 0 0 0 0 0 ...
|  $ jobself-employed            : num  0 0 0 0 0 0 0 0 0 0 ...
|   [list output truncated]

2.3.4 Train/Test datasets

Una vez finalizado el pre-procesado de los datos, se generan los dataset de entrenamiento y de test:

if (!require("caret")) {
  install.packages("caret",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(caret)

# Dividir los datos
split = 0.80
set.seed(seed.val)
ind = createDataPartition(data.norm[,ncol(data.norm)],
                          p = split, tim = 1, list = TRUE)

# Train y Test datasets
train.norm = data.norm[ind$Resample1, ]
test.norm = data.norm[-ind$Resample1, ]

# Train y Test labels
trainlabels = train.norm$y
testlabels = test.norm$y

# Eliminar datos intermedios
rm(ind, split)

Generados los dataset de entrenamiento y de test, se comprueba que las proporciones del target son semejantes en ambos dataset:

# Comprobar si el reparto de la target variable ha sido equitativo
train.prop <-round(prop.table(table(train.norm$y)), digits = 6)
test.prop <- round(prop.table(table(test.norm$y)), digits = 6)
rbind(train.prop, test.prop)

|                  no      yes
| train.prop 0.887348 0.112652
| test.prop  0.887338 0.112662

# Borrar cálculos intermedios
rm(train.prop, test.prop)

2.4 Estudio de variables categóricas

Vamos a observar la distribución de edades de los clientes contactados as?? como la distribución de los empleos de los mismos:

# ¿Qué edades tienen los clientes contactados?
if (!require("ggplot2")) {
    install.packages("ggplot2",
                     repos = "https://cloud.r-project.org/",
                     dependencies = TRUE)};
library(ggplot2)
plot <- ggplot(data.imp, aes(age)) +
    geom_histogram(aes(y=..density..),binwidth = 1) +
    stat_function(fun = dnorm,
                  color = "red",
                  args = list(mean=mean(data.imp$age),
                              sd=sd(data.imp$age))) +
    ggtitle("Distribución de edades de los clientes contactados") +
    theme_bw()

# Representación utilizando Plot.ly
if (!require("plotly")) {
  install.packages("plotly",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(plotly)
ggplotly(p = plot)

# ¿En qué trabajan los clientes contactados?
ggplot(data.imp, aes(x=job)) +
    geom_bar(aes(y= (..count..)/sum(..count..) )) +
    coord_flip() +
    ggtitle("¿En qué trabajan los clientes contactados?") +
    ylab("density") +
    theme_bw()

Podemos observar también en qué meses se contratan, proporcionalmente, más créditos. El mayor ??ndice de eficiencia se da en los meses de marzo, septiembre, octubre y diciembre, aunque la mayor cantidad de contratos se dieron entre los meses de mayo y agosto:

# Re-factorizar la variable "month"
data.imp$month.all <- factor(x = data.imp$month,
                     levels = c("jan","feb","mar","apr","may","jun",
                                "jul","aug","sep","oct","nov","dec"),
                     labels = c("enero","febrero","marzo",
                                "abril","mayo","junio",
                                "julio","agosto","septiembre",
                                "octubre","noviembre","diciembre"))

# Ratio de contratación mensual
plot1 <- ggplot(data.imp, aes(month.all, fill = y)) +
    geom_bar(position = "fill") +
    ggtitle("Ratio de contratación mensual") +
    theme_classic()

# Distribución por meses del éxito (o no) de la contratación
plot2 <- ggplot(data.imp, aes(month.all,y)) +
    geom_count(aes(size = ..prop.., group = y)) +
    scale_size_area(max_size = 20) +
    ggtitle("Distribución por meses de la contratación") +
    theme_bw()

if (!require("gridExtra")) {
  install.packages("gridExtra",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(gridExtra)
grid.arrange(plot1, plot2, nrow=2, ncol=1)

# Borrar cálculos intermedios
rm(plot, plot1, plot2)

2.5 Estudio de variables numéricas

A continuación se calcula la matriz de correlación de las variables numéricas y se representa ordenadamente en función de la importancia de cada atributo:

# Calcular la matriz de correlación de las variables numéricas
matrix.corr <- cor(data.imp[sapply(data.imp, is.numeric)],
                         use = "complete.obs")

# Visualizar la matriz de correlación
if (!require("corrplot")) {
  install.packages("corrplot",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(corrplot)
corrplot(matrix.corr, method="ellipse", tl.col = "black", tl.srt = 45,
         order = "FPC", type = "lower", mar = c(0,0,2,0),
         main = substitute("Correlación ordenada por First Principal Component"))

# Borrar cálculos intermedios
rm(matrix.corr)

Se observa que el a mayor fuerza de trabajo (“nr.employed”) y por tanto mayor tasa de variación de empleo (“emp.var.rate”), crece el IPC (“cons.price.idx”) y el euribor (“euribor3m”), principal tipo variable de referencia del mercado hipotecario.

Insight: estos indicadores son inversamente proporcionales al número de contactos telefónicos que hace este banco (“previous”), es decir, te contactan cuando por la situación económica mundial y/o del país necesitan vender productos (en este caso créditos a plazo).

2.6 Abordar un dataset no balanceado

En el dataset utilizado en este ejercicio, la variable a predecir (y) se compone de 2 clases: “yes”/“no”. Puesto que la diferencia de observaciones que pertenecen a cada clase es tan grande (≈90% -> y=“no”; ≈10% -> y=“yes”) estamos ante un dataset con el target no balanceado.

Un muy buen post acerca de cómo abordar este problema se puede consultar aqu??: 8 Tactics To Combat Imbalanced Data. A modo de resumen, las tácticas serían las siguientes:

¿Puedes recoger más datos?
Cambiar la métrica de rendimiendo/desempeño
Muestrear el dataset de nuevo: Over-sampling / Under-sampling
Generar muestras sintéticas
Probar diferentes algoritmos
Usar modelos de penalización
Adoptar una perspectiva diferente
¡Ser creativo!

Los árboles de clasificación generalmente funcionan bien con datasets no balanceados, por lo que voy a utilizarlos directamente sin realizar ningún procesado adicional (las reglas de división que se utilizan pueden forzar a abordar ambas clases). Algunos de los algoritmos de árboles de decisión más populares son C5.0, CART y Random Forest; todos ellos utilizados en el apartado Clasificación.

Mediré el rendimiento de los diferentes métodos de clasificación utilizando las métricas Precision, Recall y F-Score, además de medir del tiempo que lleva ajustar cada modelo probado.

3 Clasificación

3.1 C5.0

El algoritmo C4.5, desarrollado por Ross Quinlan, construye árboles de decisión desde un grupo de datos de entrenamiento de la misma forma en que lo hace ID3 (también desarrollado por Quinlan), usando el concepto de entropía de información. Por tanto, C4.5 está casi siempre referido como un clasificador estadístico.

En C4.5 se hicieron un número de mejoras a ID3. Algunas de ellas son:

Manejo de ambos atributos continuos y discretos.
A fin de manejar atributos continuos, C4.5 crea un umbral y luego se divide la lista en aquellos cuyo valor de atributo es superior al umbral y los que son menores o iguales a él.
Manejo de los datos de formación con valores de atributos faltantes.
C4.5 permite valores de los atributos para ser marcado como? para faltantes. Los valores faltantes de los atributos simplemente no se usa en los cálculos de la ganancia y la entropía.
Manejo de atributos con costos diferentes.
Podando árboles después de la creación.
C4.5 se remonta a través del árbol una vez que ha sido creado e intenta eliminar las ramas que no ayudan, reemplazándolos con los nodos de hoja.

Quinlan continuó con la creación del C5.0 y el See5 (C5.0 para Unix / Linux, See5 para Windows) con fines comerciales. C5.0 ofrece una serie de mejoras en el C4.5. Algunas de estas son:

Velocidad -> C5.0 es significativamente más rápido que el C4.5 (varios órdenes de magnitud).
El uso de memoria -> C5.0 es más eficiente que el C4.5.
Árboles de decisión más pequeños -> C5.0 obtiene resultados similares a C4.5 con árboles de decisión mucho más pequeños.
Soporte para boosting -> Boosting mejora los árboles y les da una mayor precisión.
Ponderación -> C5.0 le permite ponderar los distintos casos y tipos de errores de clasificación.
Winnowing -> una opción automática de C5.0 winnows los atributos para eliminar aquellos que pueden ser de poca ayuda.

[…]

Leer más vía C4.5 Algorithm - Wikipedia (EN)

3.1.1 C5.0 por defecto

if (!require("C50")) {
  install.packages("C50",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(C50)
# Ajustar el modelo
train_label_index = grep(pattern = "^y$", colnames(train.norm))
set.seed(seed.val)
tic <- proc.time()
fit.C50 <- C5.0(x = train.norm[, -train_label_index],
                y = train.norm$y)
time <- proc.time() - tic

# Predecir sobre las muestras de test
test_label_index = grep(pattern = "^y$", colnames(test.norm))
C50.pred <- predict(fit.C50, test.norm[, -test_label_index], type = "class")

# Confusion Matrix
if (!require("gmodels")) {
  install.packages("gmodels",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(gmodels)
CrossTable(x = C50.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      7098 |       666 |      7764 | 
|              |      0.91 |      0.09 |      0.94 | 
|              |      0.97 |      0.72 |           | 
|              |      0.86 |      0.08 |           | 
| -------------|-----------|-----------|-----------|
|          yes |       211 |       262 |       473 | 
|              |      0.45 |      0.55 |      0.06 | 
|              |      0.03 |      0.28 |           | 
|              |      0.03 |      0.03 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(C50.pred == testlabels)/length(testlabels)
prec <- posPredValue(C50.pred, testlabels, positive = "yes")
rec <- sensitivity(C50.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)

results <- data.frame(as.numeric(prec), as.numeric(rec), as.numeric(f),
                      as.numeric(time[3]),
                      "C5.0", "fit.C50",
                      stringsAsFactors = FALSE)
names(results) <- c("precision", "recall", "f-score", "time",
                    "algoritmo", "descripción")

Observaciones:

accuracy: el 89% de acierto general (TP+TN)
precision: el 55% de los predichos como “yes” lo son de verdad
recall: el 28% de los “yes” originales son detectados

Estableciendo sólo un ≈6% de los contactos orignales (473/8237), se consiguen el 28% de los contratos.

Para mejorar este punto, voy a añadir al algoritmo ya utilizado (C5.0) el método ensemble adaptive boosting…

3.1.2 C5.0 con adaptive boosting

Utilizando el mismo algoritmo de árboles de decisión que en el apartado anterior (C5.0), se configura para realizar 10 iteraciones con el fin de mejorar el rendimiento del mismo:

# Ajustar el modelo
set.seed(seed.val)
tic <- proc.time()
fit.C50.boost10 <- C5.0(x = train.norm[, -train_label_index],
                        y = train.norm$y,
                        trials = 10)
time <- proc.time() - tic

# Predecir sobre las muestras de test
C50.boost10.pred <- predict(fit.C50.boost10, test.norm[, -test_label_index],
                            type = "class")

# Confusion Matrix
CrossTable(x = C50.boost10.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      7089 |       657 |      7746 | 
|              |      0.92 |      0.08 |      0.94 | 
|              |      0.97 |      0.71 |           | 
|              |      0.86 |      0.08 |           | 
| -------------|-----------|-----------|-----------|
|          yes |       220 |       271 |       491 | 
|              |      0.45 |      0.55 |      0.06 | 
|              |      0.03 |      0.29 |           | 
|              |      0.03 |      0.03 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(C50.boost10.pred == testlabels)/length(testlabels)
prec <- posPredValue(C50.boost10.pred, testlabels, positive = "yes")
rec <- sensitivity(C50.boost10.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)
results <- rbind(results, c(as.numeric(prec), as.numeric(rec), as.numeric(f),
                            as.numeric(time[3]),
                            "C5.0", "fit.C50.boost10"))

Observaciones:

accuracy: el 89% de acierto general (TP+TN)
precision: el 55% de los predichos como “yes” lo son de verdad
recall: el 29% de los “yes” originales son detectados

Estableciendo sólo un ≈6% de los contactos orignales (491/8237), se consiguen el 29% de los contratos, siendo por tanto más eficiente que el modelo anterior pero quedándose aún corto en cuanto a la cantidad de contratos conseguidos.

3.1.3 C5.0 con penalización

En esta nueva aproximación se penalizan los falsos negativos (FN) 5 veces más que que los falsos positivos (FP). ¿Por qué? Porque en este caso es preferible contactar con un potencial comprador que finalmente no lo sea, a que compradores reales no sean contactados.

if (!require("gmodels")) {
  install.packages("gmodels",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(gmodels)
# Montar la matriz de penalizaciones
matrix.dim <- list(c("no","yes"), c("no","yes"))
names(matrix.dim) <- c("PREDICCION", "REAL")
error.cost <- matrix(c(0,1,5,0), nrow = 2,
                     dimnames = matrix.dim)

# Ajustar el modelo
set.seed(seed.val)
tic <- proc.time()
fit.C50.cost <- C5.0(x = train.norm[, -train_label_index],
                        y = train.norm$y,
                        costs = error.cost)
time <- proc.time() - tic

# Predecir sobre las muestras de test
C50.cost.pred <- predict(fit.C50.cost, test.norm[, -test_label_index],
                         type = "class")

# Confusion Matrix
CrossTable(x = C50.cost.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      6520 |       381 |      6901 | 
|              |      0.94 |      0.06 |      0.84 | 
|              |      0.89 |      0.41 |           | 
|              |      0.79 |      0.05 |           | 
| -------------|-----------|-----------|-----------|
|          yes |       789 |       547 |      1336 | 
|              |      0.59 |      0.41 |      0.16 | 
|              |      0.11 |      0.59 |           | 
|              |      0.10 |      0.07 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(C50.cost.pred == testlabels)/length(testlabels)
prec <- posPredValue(C50.cost.pred, testlabels, positive = "yes")
rec <- sensitivity(C50.cost.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)
results <- rbind(results, c(as.numeric(prec), as.numeric(rec), as.numeric(f),
                            as.numeric(time[3]),
                            "C5.0", "fit.C50.cost"))

Observaciones:

accuracy: el 86% de acierto general (TP+TN)
precision: el 41% de los predichos como “yes” lo son de verdad
recall: el 59% de los “yes” originales son detectados

Estableciendo sólo un ≈16% de los contactos originales (1336/8237), se consiguen el 59% de los contratos.

3.1.4 CART

CART (Classification and Regression Trees), es el término introducido por Leo Breiman para referirse a los algoritmos de árboles de decisión usados tanto en problemas predictivos de clasificación como de regresión. […]

Leer más vía Decision Tree Learning - Wikipedia (EN)

if (!require("ipred")) {
  install.packages("ipred",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(ipred)
# Ajustar el modelo
set.seed(seed.val)
tic <- proc.time()
fit.bagging10 <- bagging(y ~.,
                         data = train.norm,
                         nbagg = 10)
time <- proc.time() - tic

# Predecir sobre las muestras de test
bagging10.pred <- predict(fit.bagging10, test.norm[, -test_label_index],
                          type = "class")

# Confusion Matrix
CrossTable(x = bagging10.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      7010 |       658 |      7668 | 
|              |      0.91 |      0.09 |      0.93 | 
|              |      0.96 |      0.71 |           | 
|              |      0.85 |      0.08 |           | 
| -------------|-----------|-----------|-----------|
|          yes |       299 |       270 |       569 | 
|              |      0.53 |      0.47 |      0.07 | 
|              |      0.04 |      0.29 |           | 
|              |      0.04 |      0.03 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(bagging10.pred == testlabels)/length(testlabels)
prec <- posPredValue(bagging10.pred, testlabels, positive = "yes")
rec <- sensitivity(bagging10.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)
results <- rbind(results, c(as.numeric(prec), as.numeric(rec), as.numeric(f),
                            as.numeric(time[3]),
                            "CART", "fit.bagging10"))

Observaciones:

accuracy: el 88% de acierto general (TP+TN)
precision: el 47% de los predichos como “yes” lo son de verdad
recall: el 29% de los “yes” originales son detectados

Estableciendo sólo un ≈7% de los contactos orignales (569/8237), se consiguen el 29% de los contratos.

3.2 Random Forest

El algoritmo Random Forest es una combinación de árboles predictores tal que cada árbol depende de los valores de un vector aleatorio probado independientemente y con la misma distribución para cada uno de estos. Es una modificación sustancial de bagging que construye una larga colección de árboles no correlacionados y luego los promedia.

Cada árbol es construido usando el siguiente algoritmo:

Sea N el número de casos de prueba, M es el número de variables en el clasificador.
Sea m el número de variables de entrada a ser usado para determinar la decisión en un nodo dado; m debe ser mucho menor que M.
Elegir un conjunto de entrenamiento para este árbol y usar el resto de los casos de prueba para estimar el error.
Para cada nodo del árbol, elegir aleatoriamente m variables en las cuales basar la decisión. Calcular la mejor partición a partir de las m variables del conjunto de entrenamiento.

Para la predicción un nuevo caso es empujado hacia abajo por el árbol. Luego se le asigna la etiqueta del nodo terminal donde termina. Este proceso es iterado por todos los árboles en el ensamblado, y la etiqueta que obtenga la mayor cantidad de incidencias es reportada como la predicción.
[…]

Leer más vía Random Forest - Wikipedia (EN)

if (!require("randomForest")) {
  install.packages("randomForest",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(randomForest)
# Entrenar el modelo
set.seed(seed.val)
tic <- proc.time()
fit.forest100 <- randomForest(x = train.norm[, -train_label_index],
                              y = train.norm$y,
                              ntree = 100)
time <- proc.time() - tic

# Predecir sobre las muestras de test
forest100.pred = predict(fit.forest100, test.norm[, -test_label_index],
                         type = "class")

# Confusion Matrix
CrossTable(x = forest100.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      7118 |       668 |      7786 | 
|              |      0.91 |      0.09 |      0.95 | 
|              |      0.97 |      0.72 |           | 
|              |      0.86 |      0.08 |           | 
| -------------|-----------|-----------|-----------|
|          yes |       191 |       260 |       451 | 
|              |      0.42 |      0.58 |      0.05 | 
|              |      0.03 |      0.28 |           | 
|              |      0.02 |      0.03 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(forest100.pred == testlabels)/length(testlabels)
prec <- posPredValue(forest100.pred, testlabels, positive = "yes")
rec <- sensitivity(forest100.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)
results <- rbind(results, c(as.numeric(prec), as.numeric(rec), as.numeric(f),
                            as.numeric(time[3]),
                            "Random Forest", "fit.forest100"))

Observaciones:

accuracy: el 89% de acierto general (TP+TN)
precision: el 58% de los predichos como “yes” lo son de verdad
recall: el 28% de los “yes” originales son detectados

Estableciendo sólo un ≈5% de los contactos orignales (451/8237), se consiguen el 28% de los contratos.

3.3 k-NN

El método k-NN (k-Nearest Neighbors) es un método de clasificación supervisada que sirve para estimar la función de densidad de las predictoras x por cada clase C.

Es un método de clasificación no paramétrico que estima el valor de la función de densidad de probabilidad, o directamente la probabilidad a posteriori, de que un elemento x pertenezca a la clase C a partir de la información proporcionada por el conjunto de prototipos. En el proceso de aprendizaje no se hace ninguna suposición acerca de la distribución de las variables predictoras.

En el reconocimiento de patrones, el algoritmo k-NN es usado como método de clasificación de objetos (elementos) basado en un entrenamiento mediante ejemplos cercanos en el espacio de los elementos. k-NN es un tipo de “Lazy Learning”, donde la función se aproxima sólo localmente y todo el cómputo es diferido a la clasificación.
[…]

Leer más vía k-Nearest Neighbors Algorithm - Wikipedia (EN)

if (!require("class")) {
  install.packages("class",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(class)
# k = 3; Entrenar el modelo y predecir sobre las muestras de test
set.seed(seed.val)
tic <- proc.time()
knn3.pred <- knn(train = train.norm[,-train_label_index],
                 test = test.norm[,-test_label_index],
                 cl = train.norm$y,
                 k = 3)
time <- proc.time() - tic

# Confusion Matrix
CrossTable(x = knn3.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      6987 |       699 |      7686 | 
|              |      0.91 |      0.09 |      0.93 | 
|              |      0.96 |      0.75 |           | 
|              |      0.85 |      0.08 |           | 
| -------------|-----------|-----------|-----------|
|          yes |       322 |       229 |       551 | 
|              |      0.58 |      0.42 |      0.07 | 
|              |      0.04 |      0.25 |           | 
|              |      0.04 |      0.03 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(knn3.pred == testlabels)/length(testlabels)
prec <- posPredValue(knn3.pred, testlabels, positive = "yes")
rec <- sensitivity(knn3.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)
results <- rbind(results, c(as.numeric(prec), as.numeric(rec), as.numeric(f),
                            as.numeric(time[3]),
                            "k-NN", "knn3.pred"))

Observaciones:

accuracy: el 88% de acierto general (TP+TN)
precision: el 42% de los predichos como “yes” lo son de verdad
recall: el 25% de los “yes” originales son detectados

Estableciendo sólo un ≈7% de los contactos orignales (551/8237), se consiguen el 25% de los contratos.

3.4 SVM

Las máquinas de vector soporte (Support Vector Machines, SVMs) son un conjunto de algoritmos de aprendizaje supervisado desarrollados por Vladimir Vapnik y su equipo en los laboratorios AT&T.

Dado un conjunto de puntos, subconjunto de un conjunto mayor (espacio), en el que cada uno de ellos pertenece a una de dos posibles categorías, un algoritmo basado en SVM construye un modelo capaz de predecir si un punto nuevo (cuya categoría desconocemos) pertenece a una categoría o a la otra. Como en la mayoría de los métodos de clasificación supervisada, los datos de entrada (los puntos) son vistos como un vector p-dimensional (una lista ordenada de “p” números).

SVM busca un hiperplano que separe de forma óptima a los puntos de una clase de la de otra, que eventualmente han podido ser previamente proyectados a un espacio de dimensionalidad superior.

En ese concepto de “separación óptima” es donde reside la característica fundamental de las SVM: este tipo de algoritmos buscan el hiperplano que tenga la máxima distancia (margen) con los puntos que estén más cerca de él mismo. Por eso también a veces se les conoce a las SVM como clasificadores de margen máximo. De esta forma, los puntos del vector que son etiquetados con una categoría estarán a un lado del hiperplano y los casos que se encuentren en la otra categoría estarán al otro lado.

Los algoritmos SVM pertenecen a la familia de los clasificadores lineales. También pueden ser considerados un caso especial de la regularización de Tikhonov. […]

Leer más vía Support Vector Machine - Wikipedia (EN)

Nota: para el entrenamiento de este modelo se ha utilizado el dataset “data.svm.norm”. Este dataset, proporcionado en la misma página web que el original, es un conjunto con menos observaciones que el utilizado en el resto de casos, al que se le ha aplicado el mismo preprocesado. El motivo de su utilización es la gran carga computacional que conlleva la ejecucuón de este tipo de algoritmo.

if (!require("kernlab")) {
  install.packages("kernlab",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(kernlab)
# Cargar los datos pre-procesados
load(file = "./environment/svm_data.RData")
# Ajustar el modelo
set.seed(seed.val)
tic <- proc.time()
fit.svm <- ksvm(y ~ ., data = data.svm.norm,
                      kernel = "vanilladot")

|  Setting default kernel parameters

time <- proc.time() - tic

# Predecir sobre las muestras de test
svm.pred <- predict(fit.svm, test.norm[, -test_label_index],
                    type = "response")

# Confusion Matrix
CrossTable(x = svm.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      7215 |       747 |      7962 | 
|              |      0.91 |      0.09 |      0.97 | 
|              |      0.99 |      0.80 |           | 
|              |      0.88 |      0.09 |           | 
| -------------|-----------|-----------|-----------|
|          yes |        94 |       181 |       275 | 
|              |      0.34 |      0.66 |      0.03 | 
|              |      0.01 |      0.20 |           | 
|              |      0.01 |      0.02 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(svm.pred == testlabels)/length(testlabels)
prec <- posPredValue(svm.pred, testlabels, positive = "yes")
rec <- sensitivity(svm.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)
results <- rbind(results, c(as.numeric(prec), as.numeric(rec), as.numeric(f),
                            as.numeric(time[3]),
                            "SVM", "fit.svm"))

Observaciones:

accuracy: el 90% de acierto general (TP+TN)
precision: el 66% de los predichos como “yes” lo son de verdad
recall: el 20% de los “yes” originales son detectados

Estableciendo sólo un ≈3% de los contactos orignales (275/8237), se consiguen el 20% de los contratos.

Utilizando el mismo algoritmo, se crea a continuación un modelo aplicando diferentes pesos a las clases del target; en este caso, se da una importancia a la clase “yes” 9 veces la de la variables “no”.

# Ajustar el modelo con pesos
set.seed(seed.val)
tic <- proc.time()
fit.svm.cost <- ksvm(y ~ ., data = data.svm.norm,
                class.weights = c("yes"=9, "no"=1),
                kernel = "vanilladot")

|  Setting default kernel parameters

time <- proc.time() - tic

# Predecir sobre las muestras de test
svm.cost.pred <- predict(fit.svm.cost, test.norm[, -test_label_index],
                         type = "response")

# Confusion Matrix
CrossTable(x = svm.cost.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      6288 |       381 |      6669 | 
|              |      0.94 |      0.06 |      0.81 | 
|              |      0.86 |      0.41 |           | 
|              |      0.76 |      0.05 |           | 
| -------------|-----------|-----------|-----------|
|          yes |      1021 |       547 |      1568 | 
|              |      0.65 |      0.35 |      0.19 | 
|              |      0.14 |      0.59 |           | 
|              |      0.12 |      0.07 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(svm.cost.pred == testlabels)/length(testlabels)
prec <- posPredValue(svm.cost.pred, testlabels, positive = "yes")
rec <- sensitivity(svm.cost.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)
results <- rbind(results, c(as.numeric(prec), as.numeric(rec), as.numeric(f),
                            as.numeric(time[3]),
                            "SVM", "fit.svm.cost"))

Observaciones:

accuracy: el 83% de acierto general (TP+TN)
precision: el 35% de los predichos como “yes” lo son de verdad
recall: el 59% de los “yes” originales son detectados

Estableciendo sólo un ≈19% de los contactos orignales (1568/8237), se consiguen el 59% de los contratos. Como sucedía en el caso del árbol C5.0 con pesos, el recall del modelo aumenta considerablemente.

3.5 Neural Network

Las redes de neuronas artificiales (en inglés como: “ANN”) son un paradigma de aprendizaje y procesamiento automático inspirado en la forma en que funciona el sistema nervioso biológico. Se trata de un sistema de interconexión de neuronas que colaboran entre s?? para producir un estímulo de salida. En inteligencia artificial (IA) es frecuente referirse a ellas como redes de neuronas o redes neuronales.

Una red neuronal se compone de unidades llamadas neuronas. Cada neurona recibe una serie de entradas a través de interconexiones y emite una salida. Esta salida viene dada por tres funciones:

Una función de propagación (también conocida como función de excitación), que por lo general consiste en el sumatorio de cada entrada multiplicada por el peso de su interconexión (valor neto). Si el peso es positivo, la conexión se denomina excitatoria; si es negativo, se denomina inhibitoria.
Una función de activación, que modifica a la anterior. Puede no existir, siendo en este caso la salida la misma función de propagación.
Una función de transferencia, que se aplica al valor devuelto por la función de activación. Se utiliza para acotar la salida de la neurona y generalmente viene dada por la interpretación que queramos darle a dichas salidas. Algunas de las más utilizadas son la función sigmoidea (para obtener valores en el intervalo (0,1)) y la tangente hiperbólica (para obtener valores en el intervalo (-1,1)).

[…]

Leer más vía Artificial Neural Network - Wikipedia (EN)

# Feed-forward Neural Network con UNA ÚNICA hidden layer de 25 nodos
if (!require("nnet")) {
  install.packages("nnet",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(nnet)

# Ajustar el modelo
set.seed(seed.val)
tic <- proc.time()
fit.nn1.25 <- nnet(train.norm[,-50], class.ind(train.norm$y),
                   size = 25, softmax = TRUE,
                   maxit = 250, MaxNWts = 1500) # reltol = 0.000001
time <- proc.time() - tic

# Predecir sobre las muestras de test y factorizar la salida
nn1.25.pred <- predict(fit.nn1.25, test.norm[,-50], type = "class")
nn1.25.pred <- factor(nn1.25.pred)

# Confusion Matrix
CrossTable(x = nn1.25.pred, y = testlabels,
           prop.chisq = FALSE, prop.r = TRUE, prop.c = TRUE, digits = 2,
           dnn = c("PREDICCION", "REAL"))

| 
|  
|    Cell Contents
| |-------------------------|
| |                       N |
| |           N / Row Total |
| |           N / Col Total |
| |         N / Table Total |
| |-------------------------|
| 
|  
| Total Observations in Table:  8237 
| 
|  
|              | REAL 
|   PREDICCION |        no |       yes | Row Total | 
| -------------|-----------|-----------|-----------|
|           no |      7005 |       668 |      7673 | 
|              |      0.91 |      0.09 |      0.93 | 
|              |      0.96 |      0.72 |           | 
|              |      0.85 |      0.08 |           | 
| -------------|-----------|-----------|-----------|
|          yes |       304 |       260 |       564 | 
|              |      0.54 |      0.46 |      0.07 | 
|              |      0.04 |      0.28 |           | 
|              |      0.04 |      0.03 |           | 
| -------------|-----------|-----------|-----------|
| Column Total |      7309 |       928 |      8237 | 
|              |      0.89 |      0.11 |           | 
| -------------|-----------|-----------|-----------|
| 
|

# Performance metrics...
acc <- sum(nn1.25.pred == testlabels)/length(testlabels)
prec <- posPredValue(nn1.25.pred, testlabels, positive = "yes")
rec <- sensitivity(nn1.25.pred, testlabels, positive = "yes")
f <- (2 * prec * rec) / (prec + rec) # = (2 * TP) / (2 * TP + FP + FN)
results <- rbind(results, c(as.numeric(prec), as.numeric(rec), as.numeric(f),
                            as.numeric(time[3]),
                            "Neural Network", "fit.nn1.25"))

Observaciones:

accuracy: el 88% de acierto general (TP+TN)
precision: el 46% de los predichos como “yes” lo son de verdad
recall: el 28% de los “yes” originales son detectados

Estableciendo sólo un ≈7% de los contactos orignales (564/8237), se consiguen el 28% de los contratos.

Representación de la red neuronal. Los pesos positivos se representan en color verde y los negativos en rojo.

if (!require("NeuralNetTools")) {
  install.packages("NeuralNetTools",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};

| package 'RSNNS' successfully unpacked and MD5 sums checked
| package 'NeuralNetTools' successfully unpacked and MD5 sums checked
| 
| The downloaded binary packages are in
|   C:\Users\ivazquez\AppData\Local\Temp\RtmpaMr0sj\downloaded_packages

library(NeuralNetTools)

plotnet(mod_in = fit.nn1.25, alpha = 0.2, cex_val = 0.8,
        pos_col = "lightgreen", neg_col = "red",
        bord_col = "lightgrey", circle_col = "lightgrey")

3.6 Comparativa de modelos

Para comparar el rendimiento de los diferentes modelos creados se pueden utilizar diferentes métricas. Puesto que las clases del target están imbalanceadas y el foco está puesto en detectar a los potenciales compradores del producto (y=“yes”), descarto la métrica accuracy (que no sólo tiene en cuenta los TP sino también los TN).

Decido por tanto utilizar las métricas precision, recall y especialmente F-Score, métrica que deriva de las 2 anteriores. Más abajo se representan estos resultados de diferentes formas:

# Convertir las 4 primeras variables de 'results' en numéricas
results[,1] <- as.numeric(results[,1])
results[,2] <- as.numeric(results[,2])
results[,3] <- as.numeric(results[,3])
results[,4] <- as.numeric(results[,4])

# Precision vs. Recall
plot1 <- ggplot(data = results,
                aes(precision, recall, label = results[6])) +
    geom_label(aes(fill = factor(algoritmo)),
               colour = "white", fontface = "bold") +
    labs(title = "Precision vs. Recall",
         x = "Precision = TP/(TP+FP)", y = "Recall = TP/(TP+FN)") +
    theme(legend.position = "none") +
    coord_cartesian(xlim = c(0.25,0.85), ylim = c(0.1,0.6))

# Zoom
plot2 <- ggplot(data = results,
                aes(precision, recall, label = results[6])) +
    geom_label(aes(fill = factor(algoritmo)),
               colour = "white", fontface = "bold") +
    labs(title = "Precision vs. Recall (zoom)",
         x = "Precision = TP/(TP+FP)", y = "Recall = TP/(TP+FN)") +
    theme(legend.position = "none") +
    coord_cartesian(xlim = c(0.42,0.65), ylim = c(0.27,0.31))

grid.arrange(plot1, plot2, nrow=1, ncol=2)

En esta visualización, un modelo ideal se situaría en la esquina superior derecha, significando una precisión y un recall del 100%. Los modelos con un mejor comportamiento parecen ser los 2 en los que se han aplicado diferentes pesos/penalizaciones (fit.C50.cost y fit.svm.pesos). La diferencia entre ambos es menor; mientras uno es ligeramente superior en términos de precisión, el otro lo es en recall.

Para decidir entre ambos candidatos, se visualiza a continuación la métrica F-Score en función del tiempo del entramiento que ha supuesto cada modelo:

# F-Score vs. Tiempo de entrenamiento
ggplot(data = results,
       aes(time, `f-score`, label = results[6])) +
    geom_label(aes(fill = factor(algoritmo)),
               colour = "white", fontface = "bold") +
    labs(title = "F-Score vs. Tiempo de Entrenamiento",
         x = "Tiempo de entrenamiento [s]", y = "F-Score") +
    theme(legend.position = "none") +
    coord_cartesian(xlim = c(0,70))

La métrica F-Score lleva a decantarse finalmente por el algoritmo de árbol de decisión C5.0 con penalizaciones en las decisiones. Además, todos los modelos C5.0 han resultado ser los más rápidos, más aún teniendo en cuenta que para el entranamiento de los SVM se ha utilizado un dataset especial reducido (precisamente por su alta exigencia computacional).

3.7 Optimizar un modelo

Cada algoritmo de clasificación tiene sus propios parámetros de configuración que, afinados adecuadamente, puede resultar en un mayor rendimiento del modelo para la aplicación en la que se esté usando.

A modo de ejemplo, a continuación he “optimizado” un modelo basado en el algoritmo C5.0 con adaptive boosting:

if (!require("caret")) {
  install.packages("caret",
                   repos = "https://cloud.r-project.org/",
                   dependencies = TRUE)};
library(caret)
# Definir el proceso de entrenamiento - Ver ?trainControl
fitControl <- trainControl(method = "repeatedcv",
                           number = 10,
                           repeats = 3,
                           selectionFunction = "best")

# Data Frame con todas las combinaciones de configuraciones
grid <- expand.grid(.model = "tree",
                    .winnow = c("FALSE","TRUE"),
                    .trials = c(3, 5))

# Entrenar los direntes modelos con el algoritmo C5.0
# Performance metric: Kappa
set.seed(seed.val)
modelo <- train(y ~ ., data = data.norm,
                method="C5.0",
                tuneGrid = grid,
                trControl = fitControl,
                metric = "Kappa",
                verbose = FALSE)

El mejor modelo de los 4 probados es el C5.0 con adaptive boosting de 5 iteraciones y con feature selection.

# Mejor configuración de las propuestas
modelo$bestTune

# Resultados CON y SIN feature selection (winnow/nowinnow)
winnow <- modelo$results[c(2,4),]
nowinnow <- modelo$results[c(1,3),]

# Trials vs. Kappa
ggplot(data = winnow, aes(trials, Kappa)) +
    geom_line(colour = "blue") +
    geom_point(colour = "red", size = 2) +
    geom_line(data = nowinnow, colour = "green") +
    geom_point(data = nowinnow, colour = "red", size = 2) +
    labs(title = "Kappa vs. Trials", x = "Trials", y = "Kappa (Repeated CV)") +
    coord_cartesian(xlim = c(3,5)) +
    scale_x_continuous(breaks=c(3,5))

4 Conclusiones

El banco parece lanzar esta campaña de marketing cuando la situación económica del país es débil.
¿Merece la pena aplicar el modelo? Calcular CAC (Cost of Customer Acquisition) en cada caso:
1. Rentabilidad de contactar con 41.188 clientes colocando 4.640 créditos (11% de los contactos).
2. Rentabilidad de contactar con 6.590 clientes colocando 2.701 créditos (41% de los contactos, que además es el 59% del total posible).
¿Hasta qué punto es rentable aumentar el ratio de FP a cambio de mejorar el ratio de TP? -> Curva ROC
Otros enfoques:
1. ¿Se puede ahorrar en el gasto telefónico? ¿Qué plan se tiene (RTC vs. VoIP)?
2. Fórmulas de marketing digital (RR.SS., AdWords, …)
3. Adaptar condiciones del producto al cliente (según edad, situación laboral, …)
Muchas más que saldrían con más tiempo y sobre todo… ¡trabajando en equipo! :-)

M3.1 Business Analytics

Arkaitz Etxezarreta

16/12/2016

1 Bank Marketing dataset

1.1 Información del dataset

1.2 Atributos/Variables

2 Exploración del dataset

2.1 Preparar el entorno de trabajo

2.2 Estructura del dataset

2.3 Preprocesado

2.3.1 Factorizar y eliminar variables

2.3.2 Tratamiento de NAs

2.3.3 Normalización y binarización

2.3.4 Train/Test datasets

2.4 Estudio de variables categóricas

2.5 Estudio de variables numéricas

2.6 Abordar un dataset no balanceado

3 Clasificación

3.1 C5.0

3.1.1 C5.0 por defecto

3.1.2 C5.0 con adaptive boosting

3.1.3 C5.0 con penalización

3.1.4 CART

3.2 Random Forest

3.3 k-NN

3.4 SVM

3.5 Neural Network

3.6 Comparativa de modelos

3.7 Optimizar un modelo

4 Conclusiones