A-INTRODUCCION Y OBJETIVOS.

Este modulo de datos de muestra rastrea el abandono de clientes de una empresa de telecomunicaciones ficticia en funcion de varios factores.

La columna de abandono indica si el cliente se fue en el ultimo mes. Otras columnas incluyen genero, dependientes, cargos mensuales y muchas con informacion sobre los tipos de servicios que tiene cada cliente.

Se predice el abandono o no del cliente de la empresa.

A posteriori se realizara una evaluacion del negocio en base a la realizacion de una campaña de fidelizacion de clientes(clientes retenidos), en la cual se tomaran datos ficticios sobre la campaña y se mostraran varias opciones para determinar su mejor rendimiento posible y finalmente el ROI de la misma.

B-CONTEXTO Y DEFINICION DE LAS VARIABLES

customerID: Recoge el ID de cada cliente

VARIABLES PERSONALES:

gender: Indica si el cliente es hombre o mujer

SeniorCitizen: Indica si el cliente tiene 65 años o mas: Si, No

Partner: Indica si el cliente tiene pareja o no

Dependents: Indica si el cliente vive con algunn dependiente: Si, No.  Los dependientes pueden ser hijos, padres, abuelos, etc.

VARIABLE QUE INDICA ANTIGUEDAD:

tenure: Indica el nº de meses que el cliente ha permanecido en la cia, o que permanece

VARIABLES ASOCIADAS AL SERVICIO SUSCRITO:

PhoneService: Indica si el cliente tiene servicio telefonico o no

MultipleLines: Indica si el cliente tiene multiples lineas o no

InternetService: Indica el proveedor de servicios de internet del cliente

OnlineSecurity: Indica si el cliente tiene seguridad en linea o no

OnlineBackup: Indica si el cliente tiene respaldo en linea o no

DeviceProtection: Indica si el cliente tiene proteccion del dispositivo o no

TechSupport: Indica si el cliente tiene soporte tecnico o no

StreamingTV: Indica si el cliente tiene servicio de TV en streaming o no

StreamingMovies: Indica si el cliente dispone del servicio de peliculas en streaming o no.

VARIABLES ASOCIADAS AL TIPO DE CONTRATO:

Contract: Indica el plazo del contrato

PaperlessBilling: Indica si el cliente dispone de facturacion electronica o no

PaymentMethod: Indica el metodo de pago del cliente

VARIABLES ASOCIADAS AL IMPORTE Y TOTAL CARGADO:

MonthlyCharges:Indica el importe cobrado mensualmente al cliente

TotalCharges:Indica la cantidad total cargada al cliente

VARIABLE OBJETIVO:

Churn: Indica si el cliente abandona o no, sera la variable a predecir

0. CARGA DE LIBRERIAS, OPCIONES Y DATOS

#Carga de librerias:

library(doParallel)
library(data.table)
library(tidyverse)
library(ROCR)
library(rpart)
library(rpart.plot)
library(smbinning)
library(randomForest)
library(corrplot)
library(ggpubr)
library(cowplot)
library(tidymodels)
library(rattle)
library(RColorBrewer)
library(DataExplorer)
library(inspectdf)
library(kableExtra)

options(scipen=999)#Desactiva la notacion cientifica

#Carga de datos:

dt <- fread('Telco-Customer-Churn.csv')

1. ANALISIS EXPLORATORIO PRELIMINAR

as.data.frame(sort(names(dt))) # Visualizacion nombre de variables datos
##     sort(names(dt))
## 1             Churn
## 2          Contract
## 3        customerID
## 4        Dependents
## 5  DeviceProtection
## 6            gender
## 7   InternetService
## 8    MonthlyCharges
## 9     MultipleLines
## 10     OnlineBackup
## 11   OnlineSecurity
## 12 PaperlessBilling
## 13          Partner
## 14    PaymentMethod
## 15     PhoneService
## 16    SeniorCitizen
## 17  StreamingMovies
## 18      StreamingTV
## 19      TechSupport
## 20           tenure
## 21     TotalCharges
str(dt) # Visualizacion estructura datos 1
## Classes 'data.table' and 'data.frame':   7043 obs. of  21 variables:
##  $ customerID      : chr  "7590-VHVEG" "5575-GNVDE" "3668-QPYBK" "7795-CFOCW" ...
##  $ gender          : chr  "Female" "Male" "Male" "Male" ...
##  $ SeniorCitizen   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Partner         : chr  "Yes" "No" "No" "No" ...
##  $ Dependents      : chr  "No" "No" "No" "No" ...
##  $ tenure          : int  1 34 2 45 2 8 22 10 28 62 ...
##  $ PhoneService    : chr  "No" "Yes" "Yes" "No" ...
##  $ MultipleLines   : chr  "No phone service" "No" "No" "No phone service" ...
##  $ InternetService : chr  "DSL" "DSL" "DSL" "DSL" ...
##  $ OnlineSecurity  : chr  "No" "Yes" "Yes" "Yes" ...
##  $ OnlineBackup    : chr  "Yes" "No" "Yes" "No" ...
##  $ DeviceProtection: chr  "No" "Yes" "No" "Yes" ...
##  $ TechSupport     : chr  "No" "No" "No" "Yes" ...
##  $ StreamingTV     : chr  "No" "No" "No" "No" ...
##  $ StreamingMovies : chr  "No" "No" "No" "No" ...
##  $ Contract        : chr  "Month-to-month" "One year" "Month-to-month" "One year" ...
##  $ PaperlessBilling: chr  "Yes" "No" "Yes" "No" ...
##  $ PaymentMethod   : chr  "Electronic check" "Mailed check" "Mailed check" "Bank transfer (automatic)" ...
##  $ MonthlyCharges  : num  29.9 57 53.9 42.3 70.7 ...
##  $ TotalCharges    : num  29.9 1889.5 108.2 1840.8 151.7 ...
##  $ Churn           : chr  "No" "No" "Yes" "No" ...
##  - attr(*, ".internal.selfref")=<externalptr>
glimpse(dt) # Visualizacion estructura datos 2
## Rows: 7,043
## Columns: 21
## $ customerID       <chr> "7590-VHVEG", "5575-GNVDE", "3668-QPYBK", "7795-CFOCW~
## $ gender           <chr> "Female", "Male", "Male", "Male", "Female", "Female",~
## $ SeniorCitizen    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## $ Partner          <chr> "Yes", "No", "No", "No", "No", "No", "No", "No", "Yes~
## $ Dependents       <chr> "No", "No", "No", "No", "No", "No", "Yes", "No", "No"~
## $ tenure           <int> 1, 34, 2, 45, 2, 8, 22, 10, 28, 62, 13, 16, 58, 49, 2~
## $ PhoneService     <chr> "No", "Yes", "Yes", "No", "Yes", "Yes", "Yes", "No", ~
## $ MultipleLines    <chr> "No phone service", "No", "No", "No phone service", "~
## $ InternetService  <chr> "DSL", "DSL", "DSL", "DSL", "Fiber optic", "Fiber opt~
## $ OnlineSecurity   <chr> "No", "Yes", "Yes", "Yes", "No", "No", "No", "Yes", "~
## $ OnlineBackup     <chr> "Yes", "No", "Yes", "No", "No", "No", "Yes", "No", "N~
## $ DeviceProtection <chr> "No", "Yes", "No", "Yes", "No", "Yes", "No", "No", "Y~
## $ TechSupport      <chr> "No", "No", "No", "Yes", "No", "No", "No", "No", "Yes~
## $ StreamingTV      <chr> "No", "No", "No", "No", "No", "Yes", "Yes", "No", "Ye~
## $ StreamingMovies  <chr> "No", "No", "No", "No", "No", "Yes", "No", "No", "Yes~
## $ Contract         <chr> "Month-to-month", "One year", "Month-to-month", "One ~
## $ PaperlessBilling <chr> "Yes", "No", "Yes", "No", "Yes", "Yes", "Yes", "No", ~
## $ PaymentMethod    <chr> "Electronic check", "Mailed check", "Mailed check", "~
## $ MonthlyCharges   <dbl> 29.85, 56.95, 53.85, 42.30, 70.70, 99.65, 89.10, 29.7~
## $ TotalCharges     <dbl> 29.85, 1889.50, 108.15, 1840.75, 151.65, 820.50, 1949~
## $ Churn            <chr> "No", "No", "Yes", "No", "Yes", "Yes", "No", "No", "Y~
# Visualicemos la estructura de datos y sus dimensiones de forma mas grafica.


plot_intro(dt)

# Inspecionamos visualmente  variables categoricas(niveles):


dt %>%
  inspect_cat() %>% 
  show_plot()

Reservamos variables independientes y variable objetivo susceptibles a pasar a factores.

dt_factor <- c('gender','SeniorCitizen', 'Partner', 'Dependents', 'PhoneService', 'MultipleLines', 'InternetService', 'OnlineSecurity','OnlineBackup', 'DeviceProtection','TechSupport', 'StreamingTV', 'StreamingMovies','Contract', 
'PaperlessBilling','PaymentMethod','Churn')

2. CALIDAD DE DATOS

2.1 Estadisticos Basicos

lapply(dt,summary)
## $customerID
##    Length     Class      Mode 
##      7043 character character 
## 
## $gender
##    Length     Class      Mode 
##      7043 character character 
## 
## $SeniorCitizen
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1621  0.0000  1.0000 
## 
## $Partner
##    Length     Class      Mode 
##      7043 character character 
## 
## $Dependents
##    Length     Class      Mode 
##      7043 character character 
## 
## $tenure
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   29.00   32.37   55.00   72.00 
## 
## $PhoneService
##    Length     Class      Mode 
##      7043 character character 
## 
## $MultipleLines
##    Length     Class      Mode 
##      7043 character character 
## 
## $InternetService
##    Length     Class      Mode 
##      7043 character character 
## 
## $OnlineSecurity
##    Length     Class      Mode 
##      7043 character character 
## 
## $OnlineBackup
##    Length     Class      Mode 
##      7043 character character 
## 
## $DeviceProtection
##    Length     Class      Mode 
##      7043 character character 
## 
## $TechSupport
##    Length     Class      Mode 
##      7043 character character 
## 
## $StreamingTV
##    Length     Class      Mode 
##      7043 character character 
## 
## $StreamingMovies
##    Length     Class      Mode 
##      7043 character character 
## 
## $Contract
##    Length     Class      Mode 
##      7043 character character 
## 
## $PaperlessBilling
##    Length     Class      Mode 
##      7043 character character 
## 
## $PaymentMethod
##    Length     Class      Mode 
##      7043 character character 
## 
## $MonthlyCharges
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.25   35.50   70.35   64.76   89.85  118.75 
## 
## $TotalCharges
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    18.8   401.4  1397.5  2283.3  3794.7  8684.8      11 
## 
## $Churn
##    Length     Class      Mode 
##      7043 character character

CONCLUSION: A priori podemos comprobar que la variable TotalCharges contiene 11 NA’s.

2.2 Analisis de nulos(NA’s)

-Visualizacion grafica de posibles NA’s.

plot_missing(dt)

data.frame(colSums(is.na(dt)))
##                  colSums.is.na.dt..
## customerID                        0
## gender                            0
## SeniorCitizen                     0
## Partner                           0
## Dependents                        0
## tenure                            0
## PhoneService                      0
## MultipleLines                     0
## InternetService                   0
## OnlineSecurity                    0
## OnlineBackup                      0
## DeviceProtection                  0
## TechSupport                       0
## StreamingTV                       0
## StreamingMovies                   0
## Contract                          0
## PaperlessBilling                  0
## PaymentMethod                     0
## MonthlyCharges                    0
## TotalCharges                     11
## Churn                             0

CONCLUSION: Corroboramos que hay 11 nulos en nuestros datos en la variable TotalCharges que posteriormente hay que corregir en nuestros datos.

2.3 Analisis de ceros

Examen preliminar de nuestros datos con la funcion considerada, para comprobar si existe un problema de calidad en nuestros datos respecto al numero de ceros considerados en las variables.

contar_ceros <- function(variable) {
  temp <- transmute(dt,if_else(variable==0,1,0))
  sum(temp)
}

num_ceros <- sapply(dt,contar_ceros)
num_ceros <- data.frame(VARIABLE=names(num_ceros),CEROS = as.numeric(num_ceros),stringsAsFactors = F) 
num_ceros <- num_ceros %>%
  arrange(desc(CEROS)) %>%
  mutate(PORCENTAJE = CEROS / nrow(dt) * 100)
num_ceros
##            VARIABLE CEROS PORCENTAJE
## 1     SeniorCitizen  5901 83.7853188
## 2            tenure    11  0.1561834
## 3        customerID     0  0.0000000
## 4            gender     0  0.0000000
## 5           Partner     0  0.0000000
## 6        Dependents     0  0.0000000
## 7      PhoneService     0  0.0000000
## 8     MultipleLines     0  0.0000000
## 9   InternetService     0  0.0000000
## 10   OnlineSecurity     0  0.0000000
## 11     OnlineBackup     0  0.0000000
## 12 DeviceProtection     0  0.0000000
## 13      TechSupport     0  0.0000000
## 14      StreamingTV     0  0.0000000
## 15  StreamingMovies     0  0.0000000
## 16         Contract     0  0.0000000
## 17 PaperlessBilling     0  0.0000000
## 18    PaymentMethod     0  0.0000000
## 19   MonthlyCharges     0  0.0000000
## 20            Churn     0  0.0000000
## 21     TotalCharges    NA         NA

CONCLUSION: No apreciamos problemas de calidad de datos salvando los NA’s citados de TotalCharges que corregiremos con posterioridad.

2.4 Analisis de Atipicos

2.4.1 Analisis Atipicos tipo numerico

out <- function(variable){
  t(t(head(sort(variable,decreasing = T),20))) 
}
lapply(dt,function(x){
  if(is.double(x)) out(x)
})
## $customerID
## NULL
## 
## $gender
## NULL
## 
## $SeniorCitizen
## NULL
## 
## $Partner
## NULL
## 
## $Dependents
## NULL
## 
## $tenure
## NULL
## 
## $PhoneService
## NULL
## 
## $MultipleLines
## NULL
## 
## $InternetService
## NULL
## 
## $OnlineSecurity
## NULL
## 
## $OnlineBackup
## NULL
## 
## $DeviceProtection
## NULL
## 
## $TechSupport
## NULL
## 
## $StreamingTV
## NULL
## 
## $StreamingMovies
## NULL
## 
## $Contract
## NULL
## 
## $PaperlessBilling
## NULL
## 
## $PaymentMethod
## NULL
## 
## $MonthlyCharges
##         [,1]
##  [1,] 118.75
##  [2,] 118.65
##  [3,] 118.60
##  [4,] 118.60
##  [5,] 118.35
##  [6,] 118.20
##  [7,] 117.80
##  [8,] 117.60
##  [9,] 117.50
## [10,] 117.45
## [11,] 117.35
## [12,] 117.20
## [13,] 117.15
## [14,] 116.95
## [15,] 116.85
## [16,] 116.80
## [17,] 116.75
## [18,] 116.60
## [19,] 116.60
## [20,] 116.55
## 
## $TotalCharges
##          [,1]
##  [1,] 8684.80
##  [2,] 8672.45
##  [3,] 8670.10
##  [4,] 8594.40
##  [5,] 8564.75
##  [6,] 8547.15
##  [7,] 8543.25
##  [8,] 8529.50
##  [9,] 8496.70
## [10,] 8477.70
## [11,] 8477.60
## [12,] 8476.50
## [13,] 8468.20
## [14,] 8456.75
## [15,] 8443.70
## [16,] 8436.25
## [17,] 8425.30
## [18,] 8425.15
## [19,] 8424.90
## [20,] 8405.00
## 
## $Churn
## NULL

CONCLUSION: No apreciamos nada relevante en lo referente a los atipicos numericos.

2.4.2 Analisis Atipicos tipo integer

out <- function(variable){
  t(t(table(variable))) 
}
lapply(dt,function(x){
  if(is.integer(x)) out(x)
})
## $customerID
## NULL
## 
## $gender
## NULL
## 
## $SeniorCitizen
##         
## variable [,1]
##        0 5901
##        1 1142
## 
## $Partner
## NULL
## 
## $Dependents
## NULL
## 
## $tenure
##         
## variable [,1]
##       0    11
##       1   613
##       2   238
##       3   200
##       4   176
##       5   133
##       6   110
##       7   131
##       8   123
##       9   119
##       10  116
##       11   99
##       12  117
##       13  109
##       14   76
##       15   99
##       16   80
##       17   87
##       18   97
##       19   73
##       20   71
##       21   63
##       22   90
##       23   85
##       24   94
##       25   79
##       26   79
##       27   72
##       28   57
##       29   72
##       30   72
##       31   65
##       32   69
##       33   64
##       34   65
##       35   88
##       36   50
##       37   65
##       38   59
##       39   56
##       40   64
##       41   70
##       42   65
##       43   65
##       44   51
##       45   61
##       46   74
##       47   68
##       48   64
##       49   66
##       50   68
##       51   68
##       52   80
##       53   70
##       54   68
##       55   64
##       56   80
##       57   65
##       58   67
##       59   60
##       60   76
##       61   76
##       62   70
##       63   72
##       64   80
##       65   76
##       66   89
##       67   98
##       68  100
##       69   95
##       70  119
##       71  170
##       72  362
## 
## $PhoneService
## NULL
## 
## $MultipleLines
## NULL
## 
## $InternetService
## NULL
## 
## $OnlineSecurity
## NULL
## 
## $OnlineBackup
## NULL
## 
## $DeviceProtection
## NULL
## 
## $TechSupport
## NULL
## 
## $StreamingTV
## NULL
## 
## $StreamingMovies
## NULL
## 
## $Contract
## NULL
## 
## $PaperlessBilling
## NULL
## 
## $PaymentMethod
## NULL
## 
## $MonthlyCharges
## NULL
## 
## $TotalCharges
## NULL
## 
## $Churn
## NULL

2.4.3 Analisis Boxplot.

a1 <- ggplot(dt, aes(y = tenure)) +
  geom_boxplot(fill="blue")+ggtitle("TENURE")

a2 <- ggplot(dt, aes(y = MonthlyCharges)) +
  geom_boxplot(fill="orange")+ggtitle("MONTTHLYCHARGES")

a3 <- ggplot(dt,aes(y = TotalCharges)) +
  geom_boxplot(fill="green")+ggtitle("TOTALCHARGES")


plot_grid(a1, a2, a3,labels = "AUTO")

No se muestran visualmente en los graficos datos atipicos.

CONCLUSION: Nada a destacar fuera de lo normal.

2.5 Analisis de Coherencia

2.5.1 Analisis Longitudinal

CONCLUSION: Nada a destacar puesto que no hay mismas variables con datos a lo largo del tiempo para analizarlas de forma longitudinal.

2.5.2 Analisis de Coherencia entre variables

Hemos detectado en el analisis exploratorio preliminar variables categoricas y niveles de las mismas poco coherentes con su contenido a pares o inclusive entre mas variables:

Por ejemplo se detecta en la variable ‘Multiplelines’ que existe un nivel que ya esta definido como variable categorica anteriormente como es PhoneService con el indicativo de ‘No Phone Service’ cuando no deberia esta asociado a la variable categorica ‘Multiplelines’ sino a ‘PhoneService’.

Por otro lado las variables categoricas’OnlineSecurity’,‘OnlineBackup’, ‘DeviceProtection’, ‘TechSupport, ’StreamingMovies’,StreamingTV, hacen referencia a productos que el cliente solo puede contratar, no tiene mucho sentido tener un nivel que diga ‘No internet Service’ en todas ellas.

CONCLUSION:Lo comentado produce multicolinealidad entre las variables categoricas respectivas, para evitarla habria que corregir el nivel de ‘Multiplelines’ y cambiar No Phone Service a un nivel ‘No’ de forma aditiva al nivel de ‘No’ ya existente,y exactamente lo mismo con ‘No internet Service’ para las variables categoricas ‘OnlineSecurity’, ‘OnlineBackup’,‘DeviceProtection’,‘TechSupport’,‘StreamingMovies’, ‘StreamingTV’.

Se considera pasar la variable factor SeniorCitizen de caracter(numerico) a string, para igualarla con las demas variables independientes.

2.5.3 Comprobacion de correlacion de variables continuas.

Comprobaremos con posterioridad la correlacion entre variables continuas para ver si nos sobra alguna variable que ya este bien informada y realizaremos en el caso de haberla las correcciones oportunas.

2.6 Acciones correctoras resultado del analisis de calidad de datos

2.6.1 Accion sobre NA’s

En primer lugar vamos a ocuparnos de los NA’s suscitados en la variable ‘TotalCharges’, concretamente 11 NA’s.

En este caso podemos hacer 2 cosas: eliminar las observaciones que son 11 respecto al total de los datos que son 7043 observaciones y seria bastante licito.

o bien imputarlos como voy a hacer yo sustituyendo estos valores perdidos NA’s en este caso con la mediana, de tal forma que:

dt$TotalCharges[is.na(dt$TotalCharges)] <- median(dt$TotalCharges,na.rm = TRUE)

Si comprobamos nuevamente, vemos que ya no tenemos datos NA´s

plot_missing(dt)

2.6.2 Conversion a factores variables susceptibles a ello

Vamos a convertir a factores las variables independientes y la variable objetivo que hemos # apartado previamente.

Ademas de cambiar la variable factor ‘SeniorCitizen’ de caracter(numerico) a string. De tal manera que:

dt <- dt %>%
  mutate_at(dt_factor,.funs = factor)

levels(dt$SeniorCitizen) <- c("No","Yes")
levels(dt$Churn) <- c("0","1")

# Comprobamos los cambios a factor:

glimpse(dt)
## Rows: 7,043
## Columns: 21
## $ customerID       <chr> "7590-VHVEG", "5575-GNVDE", "3668-QPYBK", "7795-CFOCW~
## $ gender           <fct> Female, Male, Male, Male, Female, Female, Male, Femal~
## $ SeniorCitizen    <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, N~
## $ Partner          <fct> Yes, No, No, No, No, No, No, No, Yes, No, Yes, No, Ye~
## $ Dependents       <fct> No, No, No, No, No, No, Yes, No, No, Yes, Yes, No, No~
## $ tenure           <int> 1, 34, 2, 45, 2, 8, 22, 10, 28, 62, 13, 16, 58, 49, 2~
## $ PhoneService     <fct> No, Yes, Yes, No, Yes, Yes, Yes, No, Yes, Yes, Yes, Y~
## $ MultipleLines    <fct> No phone service, No, No, No phone service, No, Yes, ~
## $ InternetService  <fct> DSL, DSL, DSL, DSL, Fiber optic, Fiber optic, Fiber o~
## $ OnlineSecurity   <fct> No, Yes, Yes, Yes, No, No, No, Yes, No, Yes, Yes, No ~
## $ OnlineBackup     <fct> Yes, No, Yes, No, No, No, Yes, No, No, Yes, No, No in~
## $ DeviceProtection <fct> No, Yes, No, Yes, No, Yes, No, No, Yes, No, No, No in~
## $ TechSupport      <fct> No, No, No, Yes, No, No, No, No, Yes, No, No, No inte~
## $ StreamingTV      <fct> No, No, No, No, No, Yes, Yes, No, Yes, No, No, No int~
## $ StreamingMovies  <fct> No, No, No, No, No, Yes, No, No, Yes, No, No, No inte~
## $ Contract         <fct> Month-to-month, One year, Month-to-month, One year, M~
## $ PaperlessBilling <fct> Yes, No, Yes, No, Yes, Yes, Yes, No, Yes, No, Yes, No~
## $ PaymentMethod    <fct> Electronic check, Mailed check, Mailed check, Bank tr~
## $ MonthlyCharges   <dbl> 29.85, 56.95, 53.85, 42.30, 70.70, 99.65, 89.10, 29.7~
## $ TotalCharges     <dbl> 29.85, 1889.50, 108.15, 1840.75, 151.65, 820.50, 1949~
## $ Churn            <fct> 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,~

2.6.3 Correccion de incoherencia entre variables categoricas

dt$MultipleLines <- recode(dt$MultipleLines, "No phone service"="No")

vect <- c('OnlineSecurity', 'OnlineBackup','DeviceProtection','TechSupport', 'StreamingMovies','StreamingTV')

for (i in vect){
  
dt[[i]] <- recode(dt[[i]], "No internet service"="No")}

# Comprobamos que ya esta corregidos los niveles en nuestros datos, y nos quedarian como aparecen:

dt %>%
  inspect_cat() %>% 
  show_plot()

glimpse(dt)
## Rows: 7,043
## Columns: 21
## $ customerID       <chr> "7590-VHVEG", "5575-GNVDE", "3668-QPYBK", "7795-CFOCW~
## $ gender           <fct> Female, Male, Male, Male, Female, Female, Male, Femal~
## $ SeniorCitizen    <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, N~
## $ Partner          <fct> Yes, No, No, No, No, No, No, No, Yes, No, Yes, No, Ye~
## $ Dependents       <fct> No, No, No, No, No, No, Yes, No, No, Yes, Yes, No, No~
## $ tenure           <int> 1, 34, 2, 45, 2, 8, 22, 10, 28, 62, 13, 16, 58, 49, 2~
## $ PhoneService     <fct> No, Yes, Yes, No, Yes, Yes, Yes, No, Yes, Yes, Yes, Y~
## $ MultipleLines    <fct> No, No, No, No, No, Yes, Yes, No, Yes, No, No, No, Ye~
## $ InternetService  <fct> DSL, DSL, DSL, DSL, Fiber optic, Fiber optic, Fiber o~
## $ OnlineSecurity   <fct> No, Yes, Yes, Yes, No, No, No, Yes, No, Yes, Yes, No,~
## $ OnlineBackup     <fct> Yes, No, Yes, No, No, No, Yes, No, No, Yes, No, No, N~
## $ DeviceProtection <fct> No, Yes, No, Yes, No, Yes, No, No, Yes, No, No, No, Y~
## $ TechSupport      <fct> No, No, No, Yes, No, No, No, No, Yes, No, No, No, No,~
## $ StreamingTV      <fct> No, No, No, No, No, Yes, Yes, No, Yes, No, No, No, Ye~
## $ StreamingMovies  <fct> No, No, No, No, No, Yes, No, No, Yes, No, No, No, Yes~
## $ Contract         <fct> Month-to-month, One year, Month-to-month, One year, M~
## $ PaperlessBilling <fct> Yes, No, Yes, No, Yes, Yes, Yes, No, Yes, No, Yes, No~
## $ PaymentMethod    <fct> Electronic check, Mailed check, Mailed check, Bank tr~
## $ MonthlyCharges   <dbl> 29.85, 56.95, 53.85, 42.30, 70.70, 99.65, 89.10, 29.7~
## $ TotalCharges     <dbl> 29.85, 1889.50, 108.15, 1840.75, 151.65, 820.50, 1949~
## $ Churn            <fct> 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,~

2.6.4 Comprobacion correlacion variables continuas

corr_cont <- dt%>%select(c("tenure", "MonthlyCharges", "TotalCharges"))
dat.cor <- cor(corr_cont, method = "pearson")

corrplot(dat.cor, method = "shade",
         shade.col = NA, tl.col = "black",
         tl.srt = 90,
         addCoef.col = "black", addcolorlabel = "no",
         order = "AOE")

CONCLUSION: Vemos sobre todo que hay una correlacion alta entre ’ ‘TotalCharges’ y las demas variables continuas ‘MonthlyCharges’ y ‘tenure’. El motivo es que la variable ‘TotalCharges’ parece ser combinacion lineal de las otras dos ‘Tenure’ y ‘MonthlyCharges’, ya que su producto da como resultado la variable ‘TotalCharges’. Vamos a tomar la accion de eliminar la variable ‘TotalCharges’ de nuestros datos porque consideramos que no la necesitamos por estar ya suficientemente informada.

dte <- dt # Copiamos el dataframe antes de hacer cambios, para posterior   
          # visualizacion grafica de variables al completo.

# Eliminamos variable

dt$TotalCharges=NULL

# Comprobamos los cambios:

glimpse(dt)
## Rows: 7,043
## Columns: 20
## $ customerID       <chr> "7590-VHVEG", "5575-GNVDE", "3668-QPYBK", "7795-CFOCW~
## $ gender           <fct> Female, Male, Male, Male, Female, Female, Male, Femal~
## $ SeniorCitizen    <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, N~
## $ Partner          <fct> Yes, No, No, No, No, No, No, No, Yes, No, Yes, No, Ye~
## $ Dependents       <fct> No, No, No, No, No, No, Yes, No, No, Yes, Yes, No, No~
## $ tenure           <int> 1, 34, 2, 45, 2, 8, 22, 10, 28, 62, 13, 16, 58, 49, 2~
## $ PhoneService     <fct> No, Yes, Yes, No, Yes, Yes, Yes, No, Yes, Yes, Yes, Y~
## $ MultipleLines    <fct> No, No, No, No, No, Yes, Yes, No, Yes, No, No, No, Ye~
## $ InternetService  <fct> DSL, DSL, DSL, DSL, Fiber optic, Fiber optic, Fiber o~
## $ OnlineSecurity   <fct> No, Yes, Yes, Yes, No, No, No, Yes, No, Yes, Yes, No,~
## $ OnlineBackup     <fct> Yes, No, Yes, No, No, No, Yes, No, No, Yes, No, No, N~
## $ DeviceProtection <fct> No, Yes, No, Yes, No, Yes, No, No, Yes, No, No, No, Y~
## $ TechSupport      <fct> No, No, No, Yes, No, No, No, No, Yes, No, No, No, No,~
## $ StreamingTV      <fct> No, No, No, No, No, Yes, Yes, No, Yes, No, No, No, Ye~
## $ StreamingMovies  <fct> No, No, No, No, No, Yes, No, No, Yes, No, No, No, Yes~
## $ Contract         <fct> Month-to-month, One year, Month-to-month, One year, M~
## $ PaperlessBilling <fct> Yes, No, Yes, No, Yes, Yes, Yes, No, Yes, No, Yes, No~
## $ PaymentMethod    <fct> Electronic check, Mailed check, Mailed check, Bank tr~
## $ MonthlyCharges   <dbl> 29.85, 56.95, 53.85, 42.30, 70.70, 99.65, 89.10, 29.7~
## $ Churn            <fct> 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,~

Guardamos sesion para poder recuperar el trabajo:

save.image(file='sesion1.RData')
#load(file='sesion1.RData')--para cargar

3. ANALISIS GRAFICO TIPO DE VARIABLES.

3.1 Variables Continuas

dtc <- dte%>%select(c(tenure,MonthlyCharges,TotalCharges))

s1 <-ggplot(data = dtc, aes(x=tenure))+
  geom_histogram(color="white", fill="blue")+
  theme_bw() +
  theme(legend.position = "none")

s2 <- ggplot(data = dtc, aes(x=MonthlyCharges))+
  geom_histogram(color="white", fill="orange") +
  theme_bw() +
  theme(legend.position = "none")

s3 <- ggplot(data = dtc, aes(x=TotalCharges))+
  geom_histogram(color="white", fill="green") +
  theme_bw() +
  theme(legend.position = "none")

plot_grid(s1, s2, s3,labels = "AUTO")

plot_density(
  dtc , 
  ncol    = 2,
  title   = "Distribucion variables continuas",
  ggtheme = theme_bw(),
  theme_config = list(
    plot.title = element_text(size = 16, face = "bold"),
    strip.text = element_text(colour = "black", size = 12, face = 2)
  )
)

Podemos destacar en nuestros datos que en la variable ‘tenure’ que representa la permanencia en meses del cliente, que la mayor frecuencia de clientes la encontramos en los primeros meses de permanencia y en lo ultimos meses que son los que llevan mas tiempo (apareciendo como picos maximos en las graficas).

Tambien podemos destacar la larga cola de asimetria a la derecha de la variable ‘TotalCharges’ que indica la cantidad total cargada al cliente y que nos indica que hay pocos clientes con cantidades cargadas elevadas o muy elevadas en su totalidad, manteniendose la mayoria en un primer nivel de rango de cantidad cargada.

3.2 Variables Categoricas

dts <- dt%>%select(-customerID)

plot_bar(
  dts,
  nrow  = 5, ncol = 3,
  title = "Numero de observaciones por grupo",
  ggtheme = theme_bw(),
  theme_config = list(plot.background = element_rect(fill = "yellow"),
    plot.title = element_text(size = 16, face = "bold"),
    strip.text = element_text(colour = "black", size = 10, face = 2),
    legend.position = "none"
  )
)

3.3 Penetracion de la variable objetivo en nuestras variables.

plot_bar(dts, by="Churn", nrow=5, ncol=3)

Un dato a destacar que podria ser interesante viendo las graficas de penetracion de la target, es que hay un abandono notable de los clientes que utilizan internet de fibra optica, con respecto a los que deciden no abandonar y respecto a otros servicios de internet como DSL. Seria un componente importante a retener tratando de fidelizar al cliente.

4. TRANSFORMACION DE DATOS

4.1 Preseleccion de variables con RandomForest

dt_ind <- dt %>%select(-Churn,-customerID)%>%names()

pre_rf <- randomForest(formula = reformulate(dt_ind,'Churn'), data= dt,mtry=2,ntree=50, importance = T)
imp_rf <- randomForest::importance(pre_rf)[,4] 
imp_rf <- data.frame(VARIABLE = names(imp_rf), IMP_RF = imp_rf) 
imp_rf <- imp_rf %>% arrange(desc(IMP_RF)) %>% mutate(RANKING_RF = 1:nrow(imp_rf)) 

4.2 Preseleccion con Information Value

temp <- mutate(dt,Churn = as.numeric(as.character(Churn))) %>% as.data.frame()  
imp_iv <- smbinning.sumiv(temp[c(dt_ind,'Churn')],y="Churn")
##  
## 
  |                                                        
  |                                                  |   0%
  |                                                        
  |---                                               |   5%
  |                                                        
  |-----                                             |  11%
  |                                                        
  |--------                                          |  16%
  |                                                        
  |-----------                                       |  21%
  |                                                        
  |-------------                                     |  26%
  |                                                        
  |----------------                                  |  32%
  |                                                        
  |------------------                                |  37%
  |                                                        
  |---------------------                             |  42%
  |                                                        
  |------------------------                          |  47%
  |                                                        
  |--------------------------                        |  53%
  |                                                        
  |-----------------------------                     |  58%
  |                                                        
  |--------------------------------                  |  63%
  |                                                        
  |----------------------------------                |  68%
  |                                                        
  |-------------------------------------             |  74%
  |                                                        
  |---------------------------------------           |  79%
  |                                                        
  |------------------------------------------        |  84%
  |                                                        
  |---------------------------------------------     |  89%
  |                                                        
  |-----------------------------------------------   |  95%
  |                                                        
  |--------------------------------------------------| 100%
## 
imp_iv <- imp_iv %>% mutate(Ranking = 1:nrow(imp_iv), IV = ifelse(is.na(.$IV),0,IV)) %>% select(-Process)
names(imp_iv) <- c('VARIABLE','IMP_IV','RANKING_IV')

4.3 Preseleccion Final

imp_final <- inner_join(imp_rf,imp_iv,by='VARIABLE') %>% 
  select(VARIABLE,IMP_RF,IMP_IV,RANKING_RF,RANKING_IV) %>% 
  mutate(RANKING_TOT = RANKING_RF + RANKING_IV) %>% 
  arrange(RANKING_TOT)
imp_final
##            VARIABLE    IMP_RF IMP_IV RANKING_RF RANKING_IV RANKING_TOT
## 1            tenure 264.66089 0.8659          1          2           3
## 2          Contract 221.47715 1.2386          2          1           3
## 3    MonthlyCharges 206.25964 0.4842          3          4           7
## 4   InternetService 116.67450 0.6179          5          3           8
## 5     PaymentMethod 119.62190 0.4571          4          5           9
## 6    OnlineSecurity  48.17474 0.1719          6          7          13
## 7  PaperlessBilling  40.64527 0.2030          7          6          13
## 8       TechSupport  40.47674 0.1575          8          8          16
## 9        Dependents  29.04432 0.1555         10          9          19
## 10    SeniorCitizen  30.67409 0.1057          9         11          20
## 11          Partner  28.48917 0.1188         11         10          21
## 12     OnlineBackup  27.87895 0.0359         12         12          24
## 13 DeviceProtection  26.54301 0.0230         13         13          26
## 14      StreamingTV  25.29591 0.0202         15         14          29
## 15    MultipleLines  25.64500 0.0082         14         16          30
## 16  StreamingMovies  23.53669 0.0191         17         15          32
## 17           gender  24.78697 0.0004         16         18          34
## 18     PhoneService  13.63273 0.0008         18         17          35

Vamos a ver si los metodos son fiables para ello hacemos una correlacion entre ellos para ver si dan cosas similares.

cor(imp_final$IMP_RF,imp_final$IMP_IV)
## [1] 0.9145929

Vemos que da una correlacion alta por tanto lo consideramos confiable.

Vamos ahora a descartar las variables que no hayan salido entre las 11 primeras y la extraeremos añadiendo a mayores la variable objetivo y customerID:

dt1 <- dt %>% select(c('tenure','Contract', 'MonthlyCharges', 'InternetService', 'PaymentMethod' , 'PaperlessBilling', 'TechSupport', 'OnlineSecurity', 'Dependents', 'Partner', 'OnlineBackup', 'Churn', 'customerID'))

y guardamos nuestro dataframe.

saveRDS(dt1,'cache2.rds')
# dt1 <- readRDS('cache2.rds')--carga

5. CREACION DE VARIABLES SINTETICAS

5.1 Discretizacion

discretizar <- function(vi,target){
  temp_df <- data.frame(vi = vi, target = target)
  temp_df$target <- as.numeric(as.character(temp_df$target))
  disc <- smbinning(temp_df, y = 'target', x = 'vi')
  return(disc)
} 

5.1.1 Discretizacion automatica

# MONTHLYCHARGES:

disc_temp_MonthlyCharges <- discretizar(dt1$MonthlyCharges,dt1$Churn)
dt1_temp <- select(dt1,MonthlyCharges,Churn) 
dt1_temp <- smbinning.gen(dt1_temp,disc_temp_MonthlyCharges,chrname = 'MONTHLYCHARGES_DISC')
dt1 <- cbind(dt1,dt1_temp[,3]) %>% select(-MonthlyCharges)

5.1.2 Discretizacion manual

Escogemos para discretizacion manual la variable ‘tenure’ por su alto nivel de importancia y significaco de negocio.

Vemos la distribucion que tiene la variable:

OBJETIVOS:

-Conseguir una discretizacion de la variable lo mas similar a la distribucion original.

-Que la penetracion de la target con respecto a la variable sea monotonica(crece o decrece en todo su dominio)

# Densidad:

ggplot(dt1,aes(tenure))+geom_density()+scale_x_continuous(limits = c(0, 72))

#Histograma:

ggplot(dt1,aes(tenure))+geom_histogram(bins=25)+scale_x_continuous(limits = c(0, 72))

La variable decrece hasta un cierto punto, donde se mantiene de forma continua y al final vuelve a crecer.

Comportamiento de la variable con respecto a la target:

ggplot(dt1,aes(tenure,fill=Churn))+geom_histogram(bins=50, position='fill')+scale_x_continuous(limits = c(0, 72))

Aqui vemos al analizar en el grafico el comportamiento de la variable con respecto a la target tiene un comportamiento monotonico decreciente,vemos que cuando mas aumenta la permanencia(tenure), mas disminuye el abandono(churn).

Una vez hemos verificado ambas cosas vamos a apoyarnos en los deciles para selecionar los mejores cortes.

as.data.frame(quantile(dt1$tenure,prob = seq(0, 1, length = 11)))
##      quantile(dt1$tenure, prob = seq(0, 1, length = 11))
## 0%                                                     0
## 10%                                                    2
## 20%                                                    6
## 30%                                                   12
## 40%                                                   20
## 50%                                                   29
## 60%                                                   40
## 70%                                                   50
## 80%                                                   60
## 90%                                                   69
## 100%                                                  72

Hacemos los cortes en funcion de lo anteriormente visto:

dt1 <- dt1 %>% mutate(TENURE_DISC = as.factor(case_when(
  between(tenure,-Inf,10) ~ '01_MENOR_10',
  between(tenure,10,30) ~ '02_DE_10_A_30',
  between(tenure,30,45) ~ '03_DE_30_A_45',
  between(tenure,45,60) ~ '04_DE_45_A_60',
  between(tenure,60,Inf) ~ '05_MAYOR_60',
  TRUE ~ '00_ERROR'))
)

Eliminamos la variable ‘tenure’ original sobrante.

dt1$tenure <- NULL

Despues de probar diferentes cortes halle distribucion semejante en los cortes a la original.

ggplot(dt1,aes(TENURE_DISC)) + geom_bar()

Comprobamos si la penetracion de la target es monotonica.

ggplot(dt1,aes(TENURE_DISC, fill=Churn))+geom_bar(position='fill')

y efectivamente es monotona decreciente.

Ahora inspeccionamos visualmente todas las variables que son factores:

dt1 %>% 
  select_if(is.factor) %>% 
  gather() %>% 
  ggplot(aes(value)) +
  geom_bar() +
  facet_wrap(~ key, scales = "free") +
  theme(axis.text=element_text(size=4))

Vemos ahora la penetracion de la target en cada categoria para ver las variables que han salido monotonicas.

a <- function(var1,var2) {
  dt1_temp <- data.frame(var1 = dt1[[var1]],var2 = dt1[[var2]])
  dt1_temp %>% 
    group_by(var1) %>% 
    summarise(Conteo = n(), Porc = mean(as.numeric(as.character(var2)))) %>% 
    ggplot(aes(var1,Porc)) + geom_bar(stat='identity') + xlab(var1)
}
dt1_nombres <- dt1 %>% select_if(is.factor) %>% names()
lapply(dt1_nombres,function(x){a(x,'Churn')})
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

Como podemos ver las unicas que no tienen una penetracion de la target monotonica son ‘PaymentMethod’, ‘MonthlyCharges’ y ‘InternetService’.

Lo dejamos aqui aunque se podria hacer el mismo proceso que hicimos con tenure para las otras variables en intentar ajustarlas para que sean monotonicas.

Ahora guardamos las discretizaciones en un objeto R,por que las necesitaremos despues si queremos poner el modelo en produccion.

discretizaciones <- list(
  disc_temp_MonthlyCharges = disc_temp_MonthlyCharges
)

saveRDS(discretizaciones,'02_CortesDiscretizaciones.rds')

Reordenamos variables por simplemente por estetica y guardamos el dataframe con el que vamos a trabajar en el modelado:

centrales <- setdiff(names(dt1),c('customerID', 'Churn'))
dt1 <- dt1%>%select(
  customerID,
  one_of(centrales),
  Churn)

Vemos como ha quedado el fichero antes de la modelizacion y guardamos.

glimpse(dt1)
## Rows: 7,043
## Columns: 13
## $ customerID          <chr> "7590-VHVEG", "5575-GNVDE", "3668-QPYBK", "7795-CF~
## $ Contract            <fct> Month-to-month, One year, Month-to-month, One year~
## $ InternetService     <fct> DSL, DSL, DSL, DSL, Fiber optic, Fiber optic, Fibe~
## $ PaymentMethod       <fct> Electronic check, Mailed check, Mailed check, Bank~
## $ PaperlessBilling    <fct> Yes, No, Yes, No, Yes, Yes, Yes, No, Yes, No, Yes,~
## $ TechSupport         <fct> No, No, No, Yes, No, No, No, No, Yes, No, No, No, ~
## $ OnlineSecurity      <fct> No, Yes, Yes, Yes, No, No, No, Yes, No, Yes, Yes, ~
## $ Dependents          <fct> No, No, No, No, No, No, Yes, No, No, Yes, Yes, No,~
## $ Partner             <fct> Yes, No, No, No, No, No, No, No, Yes, No, Yes, No,~
## $ OnlineBackup        <fct> Yes, No, Yes, No, No, No, Yes, No, No, Yes, No, No~
## $ MONTHLYCHARGES_DISC <fct> 02 <= 55.95, 03 <= 68.8, 02 <= 55.95, 02 <= 55.95,~
## $ TENURE_DISC         <fct> 01_MENOR_10, 03_DE_30_A_45, 01_MENOR_10, 03_DE_30_~
## $ Churn               <fct> 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,~
saveRDS(dt1,'cache3.rds')
#dt1 <- readRDS('cache3.rds')

6. MODELIZACION

6.1 Funciones para aplicar al modelo.

Funcion para crear una matriz de confusion

confusion<-function(real,scoring,umbral){ 
  conf<-table(real,scoring>=umbral)
  if(ncol(conf)==2) return(conf) else return(NULL)
}

Funcion para calcular las metricas de los modelos: acierto, precision, cobertura y F1

metricas<-function(matriz_conf){
  acierto <- (matriz_conf[1,1] + matriz_conf[2,2]) / sum(matriz_conf) *100
  precision <- matriz_conf[2,2] / (matriz_conf[2,2] + matriz_conf[1,2]) *100
  cobertura <- matriz_conf[2,2] / (matriz_conf[2,2] + matriz_conf[2,1]) *100
  F1 <- 2*precision*cobertura/(precision+cobertura)
  salida<-c(acierto,precision,cobertura,F1)
  return(salida)
}

Funcion para probar distintos umbrales y ver el efecto sobre precision y cobertura.

umbrales<-function(real,scoring){
  umbrales<-data.frame(umbral=rep(0,times=19),acierto=rep(0,times=19),precision=rep(0,times=19),cobertura=rep(0,times=19),F1=rep(0,times=19))
  cont <- 1
  for (cada in seq(0.05,0.95,by = 0.05)){
    datos<-metricas(confusion(real,scoring,cada))
    registro<-c(cada,datos)
    umbrales[cont,]<-registro
    cont <- cont + 1
  }
  return(umbrales)
}

Funciones que calculan la curva ROC y el AUC.

roc<-function(prediction){
  r<-performance(prediction,'tpr','fpr')
  plot(r)
}

auc<-function(prediction){
  a<-performance(prediction,'auc')
  return(a@y.values[[1]])
}

6.2 Desarrollo de la Modelizacion

En primer lugar vamos a ver los datos que nos arroja en proporcion nuestra variable objetivo, por si baja de un nivel minimo aceptable alguno de sus niveles o factores y hubiera que balancear los datos (entendiendo cuestionarselo por debajo del 10% de algun nivel).

barplot(prop.table(table(dt1$Churn)),col=c("orange","blue"))

y comprobamos que no se da el caso, lo consideramos suficiente y seguimos con nuestro proceso.

Vamos a utilizar para ello el paquete tidymodels para ayudarnos en esta tarea.

Vamos a eliminar de nuestros datos ‘CustomerID’, puesto que es un identificador que no tiene ningun poder predictivo ni descriptivo.

dt3 <- dt1%>%select(-customerID)

y vamos a dividir nuestros datos en un 70 % para entrenamiento y un 30% para test. ademas para todos los modelos consideraremos un cross-validation con vfolds=5.

6.2.1 Regresion Logistica

set.seed(1971)
split_inicial <- initial_split(
  data   = dt3,
  prop   = 0.7,
  strata = Churn
)
train <- training(split_inicial)
test  <- testing(split_inicial)


mod_rl <- logistic_reg(mode = 'classification') %>%
  set_engine(engine = "glm",family=binomial(link='logit'))


mod_rcp <- recipe(Churn~.,
                  data = train)


set.seed(1981)
cv_folds <- vfold_cv(
  data    = train,
  v       = 5,
  strata  = Churn
)

wf_m <- workflow()%>%
  add_recipe(mod_rcp)%>%
  add_model(mod_rl)


my_metrics <- metric_set(accuracy,roc_auc,kap)
set.seed(1981)
rlog <- wf_m%>%
  fit_resamples(
    resamples    = cv_folds,
    metrics      = my_metrics,
    control      = control_resamples(save_pred = TRUE)
  )
collect_metrics(rlog)
## # A tibble: 3 x 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy binary     0.799     5 0.00270 Preprocessor1_Model1
## 2 kap      binary     0.445     5 0.00972 Preprocessor1_Model1
## 3 roc_auc  binary     0.843     5 0.00417 Preprocessor1_Model1
show_best(rlog, metric = "roc_auc")
## # A tibble: 1 x 6
##   .metric .estimator  mean     n std_err .config             
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1 roc_auc binary     0.843     5 0.00417 Preprocessor1_Model1
best_rlog1 <- rlog%>%
  select_best("roc_auc")

wf_final <- 
  wf_m%>%
  finalize_workflow(best_rlog1)

set.seed(2100)
final_model_rl <- fit(object = wf_final, data = train)

model_extr<- final_model_rl%>%extract_fit_engine()
summary(model_extr)
## 
## Call:
## stats::glm(formula = ..y ~ ., family = ~binomial(link = "logit"), 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8189  -0.6882  -0.2902   0.6885   3.1644  
## 
## Coefficients:
##                                      Estimate Std. Error z value
## (Intercept)                          -0.29338    0.30450  -0.964
## ContractOne year                     -0.66256    0.12655  -5.236
## ContractTwo year                     -1.80472    0.22776  -7.924
## InternetServiceFiber optic            0.75268    0.20826   3.614
## InternetServiceNo                    -1.11629    0.31148  -3.584
## PaymentMethodCredit card (automatic) -0.26035    0.13568  -1.919
## PaymentMethodElectronic check         0.36789    0.11147   3.300
## PaymentMethodMailed check            -0.13329    0.13697  -0.973
## PaperlessBillingYes                   0.42361    0.08852   4.785
## TechSupportYes                       -0.34383    0.10520  -3.268
## OnlineSecurityYes                    -0.47564    0.10211  -4.658
## DependentsYes                        -0.16178    0.10590  -1.528
## PartnerYes                            0.06647    0.09184   0.724
## OnlineBackupYes                      -0.10219    0.09253  -1.104
## MONTHLYCHARGES_DISC02 <= 55.95       -0.02616    0.29964  -0.087
## MONTHLYCHARGES_DISC03 <= 68.8        -0.47830    0.34362  -1.392
## MONTHLYCHARGES_DISC04 <= 106.75       0.06797    0.35656   0.191
## MONTHLYCHARGES_DISC05 > 106.75        0.63719    0.43667   1.459
## TENURE_DISC02_DE_10_A_30             -1.02598    0.10322  -9.939
## TENURE_DISC03_DE_30_A_45             -1.05110    0.13295  -7.906
## TENURE_DISC04_DE_45_A_60             -1.27203    0.15178  -8.381
## TENURE_DISC05_MAYOR_60               -1.69245    0.19899  -8.505
##                                                  Pr(>|z|)    
## (Intercept)                                      0.335296    
## ContractOne year                      0.00000016453097794 ***
## ContractTwo year                      0.00000000000000231 ***
## InternetServiceFiber optic                       0.000301 ***
## InternetServiceNo                                0.000339 ***
## PaymentMethodCredit card (automatic)             0.055009 .  
## PaymentMethodElectronic check                    0.000965 ***
## PaymentMethodMailed check                        0.330488    
## PaperlessBillingYes                   0.00000170711003889 ***
## TechSupportYes                                   0.001081 ** 
## OnlineSecurityYes                     0.00000318916428306 ***
## DependentsYes                                    0.126595    
## PartnerYes                                       0.469234    
## OnlineBackupYes                                  0.269420    
## MONTHLYCHARGES_DISC02 <= 55.95                   0.930439    
## MONTHLYCHARGES_DISC03 <= 68.8                    0.163947    
## MONTHLYCHARGES_DISC04 <= 106.75                  0.848821    
## MONTHLYCHARGES_DISC05 > 106.75                   0.144508    
## TENURE_DISC02_DE_10_A_30             < 0.0000000000000002 ***
## TENURE_DISC03_DE_30_A_45              0.00000000000000266 ***
## TENURE_DISC04_DE_45_A_60             < 0.0000000000000002 ***
## TENURE_DISC05_MAYOR_60               < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5703.8  on 4928  degrees of freedom
## Residual deviance: 4101.7  on 4907  degrees of freedom
## AIC: 4145.7
## 
## Number of Fisher Scoring iterations: 6

Podemos observar a la hora de interpretar los coeficientes de estimacion y signos de los mismos que en general todos se encuentran en un mismo rango de valores sin coeficientes muy desproporcionados que se salgan de ese mismo rango.

La clase predictora son los que abandonan como clientes, por tanto podemos decir que segun los resultados los que presentan menor tasa de abandono son los clientes que tienen contrato de 2 años y los clientes con una permanencia mayor de 60 meses. Entendemos que esto tiene sentido debido a la fidelidad del cliente con mas tiempo de permanencia y mayor contrato en la empresa.

Asi mismo se interpreta que los que mayor tasa de abandono presentan son los que tienen mayores cargos mensuales, esto, tambien tiene sentido, que los que soportan mayores gastos consideren abandonar por considerar caro el servicio u considerar ofertas mejores de la competencia.

y tambien los clientes con servicio de internet de fibra optica, por algun motivo de negocio que desconocemos.

Despues de analizar los p-valores vamos a escoger las siguientes variables basandonos en el criterio de significacion del 99,9 % de intervalo de confianza(***) considerando las mas predictivas de todas y eliminar las restantes, de tal forma nos quedarian:

‘Contract’,‘InternetService’,‘PaymentMethod’,‘PaperlessBilling’,‘OnlineSecurity’,‘TENURE_DISC’

Ahora aplicamos las variables escogidas en un nuevo modelo de regresion logistica.

dt4 <- dt3%>%select(-TechSupport,-Dependents,-Partner,-OnlineBackup,-MONTHLYCHARGES_DISC)
set.seed(1981)
split_inicial <- initial_split(
  data   = dt4,
  prop   = 0.7,
  strata = Churn
)
train <- training(split_inicial)
test  <- testing(split_inicial)

mod_rl <- logistic_reg(mode = 'classification') %>%
  set_engine(engine = "glm",family=binomial(link='logit'))


mod_rcp <- recipe(Churn~.,
                  data = train)


set.seed(1981)
cv_folds <- vfold_cv(
  data    = train,
  v       = 5,
  strata  = Churn
)

wf_m <- workflow()%>%
        add_recipe(mod_rcp)%>%
        add_model(mod_rl)
        

my_metrics <- metric_set(accuracy,roc_auc,kap)
set.seed(1981)
rlog <- wf_m%>%
  fit_resamples(
  resamples    = cv_folds,
  metrics      = my_metrics,
  control      = control_resamples(save_pred = TRUE)
)
collect_metrics(rlog)
## # A tibble: 3 x 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy binary     0.799     5 0.00463 Preprocessor1_Model1
## 2 kap      binary     0.433     5 0.0158  Preprocessor1_Model1
## 3 roc_auc  binary     0.839     5 0.00532 Preprocessor1_Model1
show_best(rlog, metric = "roc_auc")
## # A tibble: 1 x 6
##   .metric .estimator  mean     n std_err .config             
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1 roc_auc binary     0.839     5 0.00532 Preprocessor1_Model1
best_rlog2 <- rlog%>%
  select_best("roc_auc")

wf_final <- 
  wf_m%>%
  finalize_workflow(best_rlog2)


set.seed(2100)
final_model_rl <- fit(object = wf_final, data = train)

model_extr1 <- final_model_rl%>%extract_fit_engine()
summary(model_extr1)
## 
## Call:
## stats::glm(formula = ..y ~ ., family = ~binomial(link = "logit"), 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7274  -0.6874  -0.2986   0.7139   3.0315  
## 
## Coefficients:
##                                      Estimate Std. Error z value
## (Intercept)                          -0.51965    0.14146  -3.673
## ContractOne year                     -0.80412    0.12543  -6.411
## ContractTwo year                     -1.67381    0.20787  -8.052
## InternetServiceFiber optic            0.94559    0.09213  10.263
## InternetServiceNo                    -0.98636    0.14874  -6.632
## PaymentMethodCredit card (automatic) -0.06093    0.13529  -0.450
## PaymentMethodElectronic check         0.43820    0.11364   3.856
## PaymentMethodMailed check            -0.08777    0.13830  -0.635
## PaperlessBillingYes                   0.37299    0.08783   4.247
## OnlineSecurityYes                    -0.40244    0.09967  -4.038
## TENURE_DISC02_DE_10_A_30             -0.98493    0.09959  -9.890
## TENURE_DISC03_DE_30_A_45             -1.19957    0.12972  -9.247
## TENURE_DISC04_DE_45_A_60             -1.31725    0.14502  -9.083
## TENURE_DISC05_MAYOR_60               -1.76452    0.18698  -9.437
##                                                  Pr(>|z|)    
## (Intercept)                                      0.000239 ***
## ContractOne year                     0.000000000144884445 ***
## ContractTwo year                     0.000000000000000813 ***
## InternetServiceFiber optic           < 0.0000000000000002 ***
## InternetServiceNo                    0.000000000033227057 ***
## PaymentMethodCredit card (automatic)             0.652437    
## PaymentMethodElectronic check                    0.000115 ***
## PaymentMethodMailed check                        0.525678    
## PaperlessBillingYes                  0.000021669630671992 ***
## OnlineSecurityYes                    0.000054004096979827 ***
## TENURE_DISC02_DE_10_A_30             < 0.0000000000000002 ***
## TENURE_DISC03_DE_30_A_45             < 0.0000000000000002 ***
## TENURE_DISC04_DE_45_A_60             < 0.0000000000000002 ***
## TENURE_DISC05_MAYOR_60               < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5703.8  on 4928  degrees of freedom
## Residual deviance: 4157.8  on 4915  degrees of freedom
## AIC: 4185.8
## 
## Number of Fisher Scoring iterations: 6

y vemos que ahora las variables de la relacion tienen al menos una categoria con 3 estrellas(***).

La Prediccion del scoring del modelo (conjunto test) sera:

rl_predict <- predict(model_extr1, test, type = "response")
head(rl_predict,20)
##          2          7         10         12         15         17         18 
## 0.04681705 0.43857893 0.02957676 0.01440497 0.46250901 0.02377296 0.03183484 
##         19         24         32         33         34         35         36 
## 0.35879783 0.03921674 0.67655702 0.05738301 0.18153105 0.37293345 0.03183484 
##         41         42         49         55         57         58 
## 0.19599097 0.01715155 0.02656879 0.06114540 0.09681750 0.20890767

El resultado es el scoring de la probabilidad de un cliente de que abandone la empresa.

Veamos graficamente la pinta que tiene:

plot(rl_predict~test$Churn)

Calculamos ahora el pseudo R cuadrado, visualizando con anterioridad los coeficientes y significancia de acuerdo con lo que hemos decidido escoger para nuestro modelo:

psdo_R <- 1-(model_extr1$deviance/model_extr1$null.deviance)

psdo_R
## [1] 0.2710531

El valor que nos da no es muy elevado, este valor nos indica que con las variables que tenemos somos capaces de explicar el 27% de los motivos por los que un cliente abandona, hay mas de un 70% que no somos capaces de explicar, pero con la información disponible es lo mejor hemos podido obtener.

Ahora tenemos que transformar la probabilidad en una decision de si el cliente va a abandonar o se mantiene.

Con la funcion umbrales probamos diferentes cortes y para este caso vamos a determinar el punto de corte que maximiza a la F1(Punto intermedio entre precision y cobertura)

umb_rl<-umbrales(test$Churn,rl_predict)
umb_rl
##    umbral  acierto precision cobertura       F1
## 1    0.05 52.69631  35.79288  98.57398 52.51662
## 2    0.10 62.25166  40.91954  95.18717 57.23473
## 3    0.15 66.55629  43.77133  91.44385 59.20369
## 4    0.20 71.14475  47.58621  86.09626 61.29442
## 5    0.25 73.93567  50.56306  80.03565 61.97378
## 6    0.30 76.96310  54.75578  75.93583 63.62957
## 7    0.35 77.38884  55.72414  72.01426 62.83048
## 8    0.40 78.85525  59.28339  64.88414 61.95745
## 9    0.45 78.76064  60.64639  56.86275 58.69365
## 10   0.50 79.09177  64.76427  46.52406 54.14938
## 11   0.55 79.23368  66.66667  43.49376 52.64293
## 12   0.60 77.90918  70.79646  28.52050 40.66074
## 13   0.65 77.57805  70.61611  26.55971 38.60104
## 14   0.70 77.01041  74.50980  20.32086 31.93277
## 15   0.75 76.44276  75.20000  16.75579 27.40525
## 16   0.80  0.80000   0.80000   0.80000  0.80000
## 17   0.85  0.85000   0.85000   0.85000  0.85000
## 18   0.90  0.90000   0.90000   0.90000  0.90000
## 19   0.95  0.95000   0.95000   0.95000  0.95000

Seleccionamos la formula para hallarlo, aunque tambien podemos deducirlo de la tabla que seria el numero mas alto de la tabla F1 del cual el siguiente baja, el que de el umbral corresponderia a esa fila.

umbral_final_rl<-umb_rl[which.max(umb_rl$F1),1]
umbral_final_rl
## [1] 0.3

Hallamos la matriz de confusion del umbral optimizado.

confusion(test$Churn,rl_predict,0.3)
##     
## real FALSE TRUE
##    0  1201  352
##    1   135  426
rl_metricas<-filter(umb_rl,umbral==umbral_final_rl)
rl_metricas
##   umbral acierto precision cobertura       F1
## 1    0.3 76.9631  54.75578  75.93583 63.62957

Evaluamos la ROC

#creamos el objeto prediction
rl_prediction<-prediction(rl_predict,test$Churn)
#visualizamos la ROC
roc(rl_prediction)

#Sacamos las metricas definitivas incluyendo el AUC

rl_metricas<-cbind(rl_metricas,AUC=round(auc(rl_prediction),2)*100)
print(t(rl_metricas))
##               [,1]
## umbral     0.30000
## acierto   76.96310
## precision 54.75578
## cobertura 75.93583
## F1        63.62957
## AUC       84.00000

6.2.2 Arbol de decision

set.seed(3001)
split_inicial <- initial_split(
  data   = dt3,
  prop   = 0.7,
  strata = Churn
)
train <- training(split_inicial)
test  <- testing(split_inicial)


mod_rp <- decision_tree(
  mode = 'classification',
  cost_complexity=0.00001,
  tree_depth = tune(),
  min_n = tune()
  )%>%
  set_engine(engine = "rpart",parms = list(split="information"))

mod_rcp <- recipe(Churn~.,
                  data = train)

set.seed(3001)
cv_folds <- vfold_cv(
               data = train,
               v = 5,
               strata = Churn
             )

wf_rp <- workflow()%>%
         add_recipe(mod_rcp)%>%
         add_model(mod_rp)
registerDoParallel(cores = parallel::detectCores() - 1)

set.seed(2341)
rp_fit <- wf_rp%>%tune_bayes(
                     resamples = cv_folds,
                     initial = 15,
                     iter = 20,
                     metrics = metric_set(roc_auc),
                     control = control_bayes(save_pred=TRUE,no_improve = 10, 
                     verbose = FALSE))

stopImplicitCluster()
show_best(rp_fit, metric = "roc_auc")
## # A tibble: 5 x 9
##   tree_depth min_n .metric .estimator  mean     n std_err .config          .iter
##        <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>            <int>
## 1          8    28 roc_auc binary     0.809     5 0.00708 Iter7                7
## 2          7    11 roc_auc binary     0.808     5 0.00448 Preprocessor1_M~     0
## 3          8    38 roc_auc binary     0.808     5 0.00732 Iter6                6
## 4         10    24 roc_auc binary     0.808     5 0.00403 Preprocessor1_M~     0
## 5         13    23 roc_auc binary     0.808     5 0.00362 Iter17              17
mejor_hiper1 <- select_best(rp_fit, metric="roc_auc")
mejor_hiper1
## # A tibble: 1 x 3
##   tree_depth min_n .config
##        <int> <int> <chr>  
## 1          8    28 Iter7
wf_finalrp <- 
  wf_rp%>%
  finalize_workflow(mejor_hiper1)


set.seed(3100)
final_model_rp <- fit(object = wf_finalrp, data = train)

model_ext <- final_model_rp%>%
              extract_fit_engine()

Visualizamos el arbol del modelo.

fancyRpartPlot(model_ext)

Ahora revisamos la tabla donde podemos ver en primer termino el parametro de complejidad, error relativo,error absoluto, y desviacion estandar.

printcp(model_ext)
## 
## Classification tree:
## rpart::rpart(formula = ..y ~ ., data = data, parms = ~list(split = "information"), 
##     cp = ~0.00001, maxdepth = ~8, minsplit = min_rows(28, data))
## 
## Variables actually used in tree construction:
##  [1] Contract            Dependents          InternetService    
##  [4] MONTHLYCHARGES_DISC OnlineBackup        OnlineSecurity     
##  [7] PaperlessBilling    Partner             PaymentMethod      
## [10] TechSupport         TENURE_DISC        
## 
## Root node error: 1308/4929 = 0.26537
## 
## n= 4929 
## 
##            CP nsplit rel error  xerror     xstd
## 1  0.06090724      0   1.00000 1.00000 0.023699
## 2  0.01490826      3   0.81728 0.81728 0.022121
## 3  0.00324924      5   0.78746 0.79281 0.021877
## 4  0.00267584     10   0.76988 0.82798 0.022224
## 5  0.00229358     12   0.76453 0.81498 0.022098
## 6  0.00152905     15   0.75765 0.81040 0.022053
## 7  0.00114679     17   0.75459 0.81728 0.022121
## 8  0.00076453     19   0.75229 0.81957 0.022143
## 9  0.00050968     22   0.75000 0.82187 0.022165
## 10 0.00045872     25   0.74847 0.82569 0.022202
## 11 0.00025484     36   0.74312 0.83028 0.022247
## 12 0.00001000     39   0.74235 0.83410 0.022283

En ella nos fijamos en el parametro de complejidad desde el que tomamos de partida 0.00001 abajo del todo hasta arriba de esa relacion, tenemos que mirar en la columna xerror(error absoluto), que empieza a bajar desde arriba, En el momento que encontramos un numero donde empieza a cambiar la tendencia (a aumentar) ese seria, el parametro de complejidad que tenemos que tomar correspondiente a ese valor, y por ahi vamos a podar el arbol y ejecutar otra vez el modelo.

En este caso el parametro de complejidad que corresponderia a lo dicho anteriormente parece que minimiza aproximadamente en 0.001.

Tambien lo podemos interpretar ayudandonos viendo la representacion grafica como vemos a continuacion.

plotcp(model_ext)

Cambiamos el parametro de complejidad de partida a 0.001 e incluimos en el arbol que no tenga mas de 5 niveles.

set.seed(3001)
split_inicial <- initial_split(
  data   = dt3,
  prop   = 0.7,
  strata = Churn)

mod_rp <- decision_tree(
  mode = 'classification',
  cost_complexity=0.001,
  tree_depth =5,
  min_n = tune()
)%>%
  set_engine(engine = "rpart",parms = list(split="information"))

mod_rcp <- recipe(Churn~.,
                  data = train)

set.seed(3001)
cv_folds <- vfold_cv(
  data = train,
  v = 5,
  strata = Churn)

wf_rp <- workflow()%>%
  add_recipe(mod_rcp)%>%
  add_model(mod_rp)
registerDoParallel(cores = parallel::detectCores() - 1)

set.seed(2201)
rp_fit <- wf_rp%>%tune_bayes(
  resamples = cv_folds,
  initial = 15,
  iter = 30,
  metrics = metric_set(roc_auc),
  control = control_bayes(save_pred=TRUE,no_improve = 10, verbose = FALSE))

stopImplicitCluster()
show_best(rp_fit, metric = "roc_auc")
## # A tibble: 5 x 8
##   min_n .metric .estimator  mean     n std_err .config               .iter
##   <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                 <int>
## 1    31 roc_auc binary     0.793     5 0.00922 Preprocessor1_Model10     0
## 2    29 roc_auc binary     0.793     5 0.00922 Preprocessor1_Model13     0
## 3    30 roc_auc binary     0.793     5 0.00922 Iter1                     1
## 4    27 roc_auc binary     0.793     5 0.00922 Iter2                     2
## 5    28 roc_auc binary     0.793     5 0.00922 Iter3                     3
mejor_hiper2 <- select_best(rp_fit, metric="roc_auc")
mejor_hiper2
## # A tibble: 1 x 2
##   min_n .config              
##   <int> <chr>                
## 1    31 Preprocessor1_Model10
wf_finalrp <- 
  wf_rp%>%
  finalize_workflow(mejor_hiper2)


set.seed(3100)
final_model_rp <- fit(object = wf_finalrp, data = train)

model_extr2 <- final_model_rp%>%
  extract_fit_engine()

Visualizamos tabla complejidad-errores.

printcp(model_extr2)
## 
## Classification tree:
## rpart::rpart(formula = ..y ~ ., data = data, parms = ~list(split = "information"), 
##     cp = ~0.001, maxdepth = ~5, minsplit = min_rows(31, data))
## 
## Variables actually used in tree construction:
## [1] Contract        InternetService OnlineSecurity  PaymentMethod  
## [5] TENURE_DISC    
## 
## Root node error: 1308/4929 = 0.26537
## 
## n= 4929 
## 
##          CP nsplit rel error  xerror     xstd
## 1 0.0609072      0   1.00000 1.00000 0.023699
## 2 0.0149083      3   0.81728 0.81728 0.022121
## 3 0.0022936      5   0.78746 0.80352 0.021985
## 4 0.0010000      7   0.78287 0.80581 0.022007

Podemos ver una clara tendencia descendente de xerror(error absoluto) hasta llegar al parametro de complejidad 0.0007 de partida y tanto su error relativo como el absoluto son bastane proximos en dicho parametro y en los demas.

Por tanto consideramos como buena dicha aproximacion.

Visualizamos grafica error validacion cruzada y complejidad.

plotcp(model_extr2)

Visualizamos el arbol del modelo.

fancyRpartPlot(model_extr2)

Reglas que podrian ser utiizadas para hacer implantacion del arbol.(y es Churn–la variable objetivo(dependiente) en las reglas del arbol)

rpart.rules(model_extr2, style = 'tall', cover = T,roundint=FALSE)
## ..y is 0.07 with cover 45% when
##     Contract is One year or Two year
## 
## ..y is 0.25 with cover 0% when
##     Contract is Month-to-month
##     InternetService is Fiber optic
##     TENURE_DISC is 01_MENOR_10
##     PaymentMethod is Mailed check
##     OnlineSecurity is Yes
## 
## ..y is 0.29 with cover 25% when
##     Contract is Month-to-month
##     InternetService is DSL or No
## 
## ..y is 0.34 with cover 8% when
##     Contract is Month-to-month
##     InternetService is Fiber optic
##     TENURE_DISC is 02_DE_10_A_30 or 03_DE_30_A_45 or 04_DE_45_A_60 or 05_MAYOR_60
##     PaymentMethod is Bank transfer (automatic) or Credit card (automatic) or Mailed check
## 
## ..y is 0.42 with cover 5% when
##     Contract is Month-to-month
##     InternetService is Fiber optic
##     TENURE_DISC is 03_DE_30_A_45 or 04_DE_45_A_60 or 05_MAYOR_60
##     PaymentMethod is Electronic check
## 
## ..y is 0.54 with cover 1% when
##     Contract is Month-to-month
##     InternetService is Fiber optic
##     TENURE_DISC is 01_MENOR_10
##     PaymentMethod is Bank transfer (automatic) or Credit card (automatic) or Electronic check
##     OnlineSecurity is Yes
## 
## ..y is 0.58 with cover 5% when
##     Contract is Month-to-month
##     InternetService is Fiber optic
##     TENURE_DISC is 02_DE_10_A_30
##     PaymentMethod is Electronic check
## 
## ..y is 0.73 with cover 11% when
##     Contract is Month-to-month
##     InternetService is Fiber optic
##     TENURE_DISC is 01_MENOR_10
##     OnlineSecurity is No

Vemos el numero de nodos final por cliente.

rp_numnod <- rpart.predict(model_extr2,test,nn=T)

head(rp_numnod)
##            0          1 nn
## 4  0.9279238 0.07207616  2
## 9  0.4235294 0.57647059 59
## 12 0.9279238 0.07207616  2
## 13 0.9279238 0.07207616  2
## 15 0.4235294 0.57647059 59
## 18 0.9279238 0.07207616  2

Calculo de prediccion del modelo (scorings).

rp_predict <- predict(model_extr2, test, type = 'prob')[,2]
head(rp_predict,20)
##          4          9         12         13         15         18         21 
## 0.07207616 0.57647059 0.07207616 0.07207616 0.57647059 0.07207616 0.28662930 
##         22         23         24         28         33         35         38 
## 0.07207616 0.28662930 0.07207616 0.28662930 0.07207616 0.28662930 0.33587786 
##         41         42         53         54         57         60 
## 0.07207616 0.07207616 0.28662930 0.72916667 0.07207616 0.07207616

Veamos graficamente como se ve:

plot(rp_predict~test$Churn)

Probamos cortes con la funcion umbrales, visualizamos la tabla respectiva.

umb_rp<-umbrales(test$Churn,rp_predict)
umb_rp
##    umbral  acierto precision cobertura       F1
## 1    0.05  0.05000   0.05000   0.05000  0.05000
## 2    0.10 66.84011  43.92361  90.19608 59.07764
## 3    0.15 66.84011  43.92361  90.19608 59.07764
## 4    0.20 66.84011  43.92361  90.19608 59.07764
## 5    0.25 66.84011  43.92361  90.19608 59.07764
## 6    0.30 77.43614  56.46154  65.41889 60.61107
## 7    0.35 79.37559  63.10273  53.65419 57.99615
## 8    0.40 79.37559  63.10273  53.65419 57.99615
## 9    0.45 79.89593  67.43590  46.88057 55.31020
## 10   0.50 79.89593  67.43590  46.88057 55.31020
## 11   0.55 79.80132  67.91444  45.27629 54.33155
## 12   0.60 78.80795  74.04255  31.01604 43.71859
## 13   0.65 78.80795  74.04255  31.01604 43.71859
## 14   0.70 78.80795  74.04255  31.01604 43.71859
## 15   0.75  0.75000   0.75000   0.75000  0.75000
## 16   0.80  0.80000   0.80000   0.80000  0.80000
## 17   0.85  0.85000   0.85000   0.85000  0.85000
## 18   0.90  0.90000   0.90000   0.90000  0.90000
## 19   0.95  0.95000   0.95000   0.95000  0.95000

Lo hallamos directamente con la siguiente sintaxis.

umbral_final_rp<-umb_rp[which.max(umb_rp$F1),1]
umbral_final_rp
## [1] 0.3

Hallamos la matriz de confusion del umbral optimizado.

confusion(test$Churn,rp_predict,0.3)
##     
## real FALSE TRUE
##    0  1270  283
##    1   194  367
rp_metricas<-filter(umb_rp,umbral==umbral_final_rp)
rp_metricas
##   umbral  acierto precision cobertura       F1
## 1    0.3 77.43614  56.46154  65.41889 60.61107

Evaluamos la ROC

#creamos el objeto prediction
rp_prediction<-prediction(rp_predict,test$Churn)
#visualizamos la ROC
roc(rp_prediction)

Sacamos las metricas definitivas incluyendo el AUC

rp_metricas<-cbind(rp_metricas,AUC=round(auc(rp_prediction),2)*100)
print(t(rp_metricas))
##               [,1]
## umbral     0.30000
## acierto   77.43614
## precision 56.46154
## cobertura 65.41889
## F1        60.61107
## AUC       82.00000

6.2.3 Random Forest

mod_rf <- rand_forest(
  mode  = "classification",
  mtry  = tune(),
  trees = 500,
  min_n = tune()
) %>%
  set_engine(engine = "randomForest",importance=T)

set.seed(3001)
cv_folds <- vfold_cv(
  data = train,
  v = 5,
  strata = Churn)

mod_rcp <- recipe(Churn~.,
                  data = train)

wf_rf <- workflow()%>%
  add_recipe(mod_rcp)%>%
  add_model(mod_rf)

hiperpar_grid=grid_max_entropy(
  # Rango de busqueda para cada hiperparametro
  mtry(range = c(1L, 10L), trans = NULL),
  min_n(range = c(2L, 100L), trans = NULL),
  # Numero de combinaciones totales
  size = 100
)
registerDoParallel(cores = parallel::detectCores() - 1)

set.seed(2051)
rf_fit <- wf_rf%>%tune_grid(
  resamples = cv_folds,
  metrics = metric_set(roc_auc),
  control = control_resamples(save_pred = TRUE),
  grid = hiperpar_grid)

stopImplicitCluster()
show_best(rf_fit, metric = "roc_auc")
## # A tibble: 5 x 8
##    mtry min_n .metric .estimator  mean     n std_err .config              
##   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
## 1     1    42 roc_auc binary     0.827     5 0.00325 Preprocessor1_Model70
## 2     1    19 roc_auc binary     0.826     5 0.00346 Preprocessor1_Model80
## 3     1    92 roc_auc binary     0.826     5 0.00347 Preprocessor1_Model89
## 4     1    81 roc_auc binary     0.825     5 0.00318 Preprocessor1_Model77
## 5     2    46 roc_auc binary     0.824     5 0.00354 Preprocessor1_Model96
mejor_hiper3 <- select_best(rf_fit, metric="roc_auc")
mejor_hiper3
## # A tibble: 1 x 3
##    mtry min_n .config              
##   <int> <int> <chr>                
## 1     1    42 Preprocessor1_Model70
wf_finalrf <- 
  wf_rf%>%
  finalize_workflow(mejor_hiper3)


set.seed(3100)
final_model_rf <- fit(object = wf_finalrf, data = train)

model_xtr <- final_model_rf%>%
  extract_fit_engine()

Visualizamos las variables de nuestros datos mas importantes en un grafico.

varImpPlot(model_xtr)

Como hay 2 criterios, creamos una unica variable agregada para tener una mejor idea de la importancia de cada variable.

importancia <- randomForest::importance(model_xtr)[,3:4]
#normalizamos para poner las dos variables en la misma escala. los valores negativos son las que menos predicen y los positivos las que mas
importancia_norm <- as.data.frame(scale(importancia))
#creamos una única variable como suma de las otras
importancia_norm <- importancia_norm %>% mutate(
  Variable = rownames(importancia_norm),
  Imp_tot = MeanDecreaseAccuracy + MeanDecreaseGini) %>%
  mutate(Imp_tot = Imp_tot + abs(min(Imp_tot))) %>% 
  arrange(desc(Imp_tot)) %>% 
  select(Variable,Imp_tot,MeanDecreaseAccuracy,MeanDecreaseGini)
#hacemos un gráfico para ver la curva de caída de importancia
ggplot(importancia_norm, aes(reorder(Variable,-Imp_tot),Imp_tot)) + geom_bar(stat = "identity") + theme(axis.text.x = element_text(angle = 90,size = 7))

importancia_norm
##                                Variable   Imp_tot MeanDecreaseAccuracy
## TENURE_DISC                 TENURE_DISC 5.4236211           1.49575077
## Contract                       Contract 5.2713893           0.67705055
## PaymentMethod             PaymentMethod 4.5010238           0.79539834
## InternetService         InternetService 4.2146883           0.71651601
## MONTHLYCHARGES_DISC MONTHLYCHARGES_DISC 3.5369080           0.26410281
## OnlineSecurity           OnlineSecurity 2.3577771           0.28371930
## PaperlessBilling       PaperlessBilling 2.0165402          -0.19271414
## TechSupport                 TechSupport 2.0141544          -0.09635295
## OnlineBackup               OnlineBackup 1.0460069          -0.66157242
## Dependents                   Dependents 0.6261723          -1.37779260
## Partner                         Partner 0.0000000          -1.90410567
##                     MeanDecreaseGini
## TENURE_DISC                1.1089356
## Contract                   1.7754041
## PaymentMethod              0.8866908
## InternetService            0.6792376
## MONTHLYCHARGES_DISC        0.4538705
## OnlineSecurity            -0.7448769
## PaperlessBilling          -0.6096803
## TechSupport               -0.7084273
## OnlineBackup              -1.1113553
## Dependents                -0.8149698
## Partner                   -0.9148290

Vemos el grafico y nuestros datos en la tabla de forma decreciente segun nivel de importancia, y vemos donde puede estar un corte claro diferencia de una con respecto a la anterior, y vemos de la variable ‘TechSupport’ a ‘OnlineBackup’ un corte interesante.

Vamos a deseleccionar de nuestros datos esas 3 variables sobrantes por debajo del corte de ‘TechSupport’ y nos quedaremos con esta y las de la parte superior.

dt5 <- dt3%>%select(-OnlineBackup,-Dependents,-Partner)

y ahora vamos a realizar de nuevo el modelo con las nuevas variables.

set.seed(5001)
split_inicial <- initial_split(
  data   = dt5,
  prop   = 0.7,
  strata = Churn
)
train <- training(split_inicial)
test  <- testing(split_inicial)

mod_rf <- rand_forest(
  mode  = "classification",
  mtry  = tune(),
  trees = 500,
  min_n = tune()
) %>%
  set_engine(engine = "randomForest",importance=T)

mod_rcp <- recipe(Churn~.,
                  data = train)

set.seed(3001)
cv_folds <- vfold_cv(
  data = train,
  v = 5,
  strata = Churn)

wf_rf <- workflow()%>%
  add_recipe(mod_rcp)%>%
  add_model(mod_rf)

hiperpar_grid=grid_max_entropy(
  # Rango de busqueda para cada hiperparametro
  mtry(range = c(1L, 10L), trans = NULL),
  min_n(range = c(2L, 100L), trans = NULL),
  # Numero de combinaciones totales
  size = 100
)
registerDoParallel(cores = parallel::detectCores() - 1)


set.seed(2051)
rf_fit <- wf_rf%>%tune_grid(
  resamples = cv_folds,
  metrics = metric_set(roc_auc),
  control = control_resamples(save_pred = TRUE),
  grid = hiperpar_grid)

stopImplicitCluster()
show_best(rf_fit, metric = "roc_auc")
## # A tibble: 5 x 8
##    mtry min_n .metric .estimator  mean     n std_err .config              
##   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
## 1     2    21 roc_auc binary     0.833     5 0.00664 Preprocessor1_Model60
## 2     2    37 roc_auc binary     0.831     5 0.00631 Preprocessor1_Model66
## 3     1    81 roc_auc binary     0.831     5 0.00611 Preprocessor1_Model77
## 4     2    29 roc_auc binary     0.831     5 0.00591 Preprocessor1_Model50
## 5     2    46 roc_auc binary     0.831     5 0.00552 Preprocessor1_Model96
mejor_hiper4 <- select_best(rf_fit, metric="roc_auc")
mejor_hiper4
## # A tibble: 1 x 3
##    mtry min_n .config              
##   <int> <int> <chr>                
## 1     2    21 Preprocessor1_Model60
wf_finalrf <- 
  wf_rf%>%
  finalize_workflow(mejor_hiper4)


set.seed(3100)
final_model_rf <- fit(object = wf_finalrf, data = train)

model_extr3 <- final_model_rf%>%
  extract_fit_engine()

Calculo de prediccion del modelo (scorings).

rf_predict <- predict(model_extr3, test, type = 'prob')[,2]
head(rf_predict,20)
##     1     3     5     7     8     9    13    15    16    17    25    26    27 
## 0.798 0.150 0.988 0.382 0.054 0.514 0.012 0.286 0.000 0.000 0.012 0.054 0.508 
##    28    31    34    37    40    42    45 
## 0.584 0.000 0.024 0.988 0.428 0.000 0.068

Veamos graficamente como se ve:

plot(rf_predict~test$Churn)

Probamos cortes con la funcion umbrales, visualizamos la tabla respectiva.

umb_rf<-umbrales(test$Churn,rf_predict)
umb_rf
##    umbral  acierto precision cobertura       F1
## 1    0.05 72.18543  48.53420  79.67914 60.32389
## 2    0.10 75.82781  53.22997  73.44029 61.72285
## 3    0.15 77.05771  55.41311  69.34046 61.59937
## 4    0.20 78.28761  57.79817  67.37968 62.22222
## 5    0.25 78.99716  59.73378  63.99287 61.79002
## 6    0.30 79.37559  60.98418  61.85383 61.41593
## 7    0.35 79.75402  62.29205  60.07130 61.16152
## 8    0.40 79.65941  63.50515  54.90196 58.89101
## 9    0.45 79.94324  64.73118  53.65419 58.67446
## 10   0.50 79.99054  65.06550  53.11943 58.48871
## 11   0.55 79.37559  66.06684  45.81105 54.10526
## 12   0.60 79.56481  68.69565  42.24599 52.31788
## 13   0.65 79.28098  68.35821  40.81996 51.11607
## 14   0.70 79.32829  68.67470  40.64171 51.06383
## 15   0.75 79.13907  69.10828  38.68093 49.60000
## 16   0.80 78.85525  69.93007  35.65062 47.22550
## 17   0.85 78.85525  70.21277  35.29412 46.97509
## 18   0.90 77.38884  75.77640  21.74688 33.79501
## 19   0.95 76.77389  75.73529  18.36007 29.55524

Lo hallamos directamente con la siguiente sintaxis.

umbral_final_rf<-umb_rf[which.max(umb_rf$F1),1]
umbral_final_rf
## [1] 0.2

Hallamos la matriz de confusion del umbral optimizado.

confusion(test$Churn,rf_predict,0.2)
##     
## real FALSE TRUE
##    0  1277  276
##    1   183  378
rf_metricas<-filter(umb_rf,umbral==umbral_final_rf)
rf_metricas
##   umbral  acierto precision cobertura       F1
## 1    0.2 78.28761  57.79817  67.37968 62.22222

Evaluamos la ROC

#creamos el objeto prediction
rf_prediction<-prediction(rf_predict,test$Churn)
#visualizamos la ROC
roc(rf_prediction)

Sacamos las metricas definitivas incluyendo el AUC

rf_metricas<-cbind(rf_metricas,AUC=round(auc(rf_prediction),2)*100)
print(t(rf_metricas))
##               [,1]
## umbral     0.20000
## acierto   78.28761
## precision 57.79817
## cobertura 67.37968
## F1        62.22222
## AUC       83.00000

6.2.4 XGBOOST (extreme gradient boosting)

set.seed(6000)
split_inicial <- initial_split(
  data   = dt3,
  prop   = 0.7,
  strata = Churn
)
train <- training(split_inicial)
test  <- testing(split_inicial)



mod_xg <- boost_tree(mtry = tune(),
                     min_n = tune(),
                     tree_depth = tune(),
                     trees = 500,
                     learn_rate = 0.1,
                     sample_size = 0.8) %>% 
  set_engine("xgboost") %>% 
  set_mode("classification")


set.seed(3001)
cv_folds <- vfold_cv(
  data = train,
  v = 5,
  strata = Churn)


wf_xg<- workflow()%>%
  add_formula(Churn~.)%>%
  add_model(mod_xg)


# Parametros para tuning
params_xgb <- parameters(
  finalize(mtry(), x = train[, -1]),
  min_n(range = c(2L, 50L)),
  tree_depth(range = c(3L, 8L))
)

# Grid
set.seed(2000)
grid_xgb <- params_xgb %>% 
  grid_max_entropy(size = 10)
registerDoParallel(cores = parallel::detectCores() - 1)

set.seed(2631)
xg_fit <- wf_xg%>%tune_grid(
  resamples = cv_folds,
  metrics = metric_set(roc_auc),
  grid = grid_xgb,
  control = control_grid(save_pred = TRUE))


stopImplicitCluster()
show_best(xg_fit, metric = "roc_auc")
## # A tibble: 5 x 9
##    mtry min_n tree_depth .metric .estimator  mean     n std_err .config         
##   <int> <int>      <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>           
## 1     1    37          7 roc_auc binary     0.840     5 0.00210 Preprocessor1_M~
## 2     2    36          4 roc_auc binary     0.838     5 0.00241 Preprocessor1_M~
## 3     1     8          3 roc_auc binary     0.838     5 0.00240 Preprocessor1_M~
## 4     7    48          8 roc_auc binary     0.837     5 0.00186 Preprocessor1_M~
## 5     9    29          3 roc_auc binary     0.835     5 0.00188 Preprocessor1_M~
mejor_hiper5 <- select_best(xg_fit, metric="roc_auc")
mejor_hiper5
## # A tibble: 1 x 4
##    mtry min_n tree_depth .config              
##   <int> <int>      <int> <chr>                
## 1     1    37          7 Preprocessor1_Model10
wf_finalxg <- 
  wf_rf%>%
  finalize_workflow(mejor_hiper5)


set.seed(3300)
final_model_xg <- fit(object = wf_finalxg, data = train)

model_extr4 <- final_model_xg%>%
  extract_fit_engine()

Calculo de prediccion del modelo (scorings).

xg_predict <- predict(model_extr4, test, type = 'prob')[,2]
head(xg_predict,20)
##     1     2     3     4     6     8     9    17    19    24    25    28    31 
## 0.486 0.000 0.118 0.000 0.766 0.058 0.316 0.000 0.128 0.002 0.002 0.344 0.000 
##    32    37    39    47    50    51    58 
## 0.534 0.766 0.472 0.174 0.000 0.472 0.036

Veamos graficamente como se ve:

plot(xg_predict~test$Churn)

Probamos cortes con la funcion umbrales, visualizamos la tabla respectiva.

umb_xg<-umbrales(test$Churn,xg_predict)
umb_xg
##    umbral  acierto precision cobertura       F1
## 1    0.05 69.77294  46.14625  83.24421 59.37699
## 2    0.10 75.07096  52.16837  72.90553 60.81784
## 3    0.15 77.57805  56.50224  67.37968 61.46341
## 4    0.20 78.94986  60.35714  60.24955 60.30330
## 5    0.25 78.61873  60.96579  54.01070 57.27788
## 6    0.30 79.13907  63.27434  50.98039 56.46594
## 7    0.35 79.47020  65.91479  46.88057 54.79167
## 8    0.40 79.51750  66.84211  45.27629 53.98512
## 9    0.45 79.61211  68.36158  43.13725 52.89617
## 10   0.50 79.32829  72.79412  35.29412 47.53902
## 11   0.55 78.66604  72.17742  31.90731 44.25216
## 12   0.60 77.15232  75.65789  20.49911 32.25806
## 13   0.65 76.72658  76.74419  17.64706 28.69565
## 14   0.70 76.72658  76.74419  17.64706 28.69565
## 15   0.75 76.72658  76.74419  17.64706 28.69565
## 16   0.80  0.80000   0.80000   0.80000  0.80000
## 17   0.85  0.85000   0.85000   0.85000  0.85000
## 18   0.90  0.90000   0.90000   0.90000  0.90000
## 19   0.95  0.95000   0.95000   0.95000  0.95000

Lo hallamos directamente con la siguiente sintaxis.

umbral_final_xg<-umb_xg[which.max(umb_xg$F1),1]
umbral_final_xg
## [1] 0.15

Hallamos la matriz de confusion del umbral optimizado.

confusion(test$Churn,xg_predict,0.15)
##     
## real FALSE TRUE
##    0  1261  292
##    1   183  378
xg_metricas<-filter(umb_xg,umbral==umbral_final_xg)
xg_metricas
##   umbral  acierto precision cobertura       F1
## 1   0.15 77.57805  56.50224  67.37968 61.46341

Evaluamos la ROC

#creamos el objeto prediction
xg_prediction<-prediction(xg_predict,test$Churn)
#visualizamos la ROC
roc(xg_prediction)

Sacamos las metricas definitivas incluyendo el AUC

xg_metricas<-cbind(xg_metricas,AUC=round(auc(xg_prediction),2)*100)
print(t(xg_metricas))
##               [,1]
## umbral     0.15000
## acierto   77.57805
## precision 56.50224
## cobertura 67.37968
## F1        61.46341
## AUC       82.00000

6.3 Comparativa de metricas y curvas ROC.

6.3.1 Comparativa de metricas

comparativa <- rbind(rl_metricas,rp_metricas,rf_metricas,xg_metricas)
rownames(comparativa) <- c('Regresion Logistica', 'Arbol Decision', 'Random Forest', 'Extreme Gradient Boosting')
t(comparativa)
##           Regresion Logistica Arbol Decision Random Forest
## umbral                0.30000        0.30000       0.20000
## acierto              76.96310       77.43614      78.28761
## precision            54.75578       56.46154      57.79817
## cobertura            75.93583       65.41889      67.37968
## F1                   63.62957       60.61107      62.22222
## AUC                  84.00000       82.00000      83.00000
##           Extreme Gradient Boosting
## umbral                      0.15000
## acierto                    77.57805
## precision                  56.50224
## cobertura                  67.37968
## F1                         61.46341
## AUC                        82.00000

6.3.2 Comparativas de curvas ROC de los modelos.

rl_auc <- rlog%>%collect_predictions(summarize=TRUE,parameters=best_rlog2)%>%
  roc_curve(Churn,.pred_0)%>%
  mutate(model = "Regresion Logistica")

rp_auc <- rp_fit%>%collect_predictions(summarize=TRUE,parameters=mejor_hiper2)%>%
  roc_curve(Churn,.pred_0)%>%
  mutate(model = "Arbol de Decision") 

rf_auc <- rf_fit%>%collect_predictions(summarize=TRUE,parameters=mejor_hiper4)%>%
  roc_curve(Churn,.pred_0)%>%
  mutate(model = "Random Forest") 

xg_auc <- xg_fit%>%collect_predictions(summarize=TRUE,parameters=mejor_hiper5)%>%
  roc_curve(Churn,.pred_0)%>%
  mutate(model = "Xgboost")

bind_rows(rl_auc,rp_auc,rf_auc,xg_auc)%>% # Dibuja las 4 curvas AUC juntas
  ggplot(aes(x=1-specificity, y = sensitivity, col = model))+ # Especifica el eje x e y dibuja la columna para usar el nombre de la metrica
  geom_path(lwd = 1.5, alpha = 0.8) + # Conecta a las 4 AUC, lwd = anchura de linea, alpha = Color transparencia del valor.
  geom_abline(lty = 3) + # abline del dibujo, lty= linea tipo
  coord_equal() + # Asegura que los rangos de los ejes sean iguales
  scale_color_viridis_d(option = "plasma", end = .6)

En base a la comparativa el modelo de Regresion Logistica ha obtenido levemente mejores resultados, asi que seleccionamos este modelo ganador en la comparativa.

6.4 Guardado del modelo seleccionado.

Añadimos la variable scoring a nuestro dataframe en base a la prediccion del modelo seleccionado y guardamos el modelo final.

dt1$SCORING_CHURN <- predict(model_extr1,dt1,type='response')

# Visualizamos tabla con el scoring reflejado:

kable(head(dt1,6), booktabs = T) %>%
  kable_styling(font_size=7)
customerID Contract InternetService PaymentMethod PaperlessBilling TechSupport OnlineSecurity Dependents Partner OnlineBackup MONTHLYCHARGES_DISC TENURE_DISC Churn SCORING_CHURN
7590-VHVEG Month-to-month DSL Electronic check Yes No No No Yes Yes 02 <= 55.95 01_MENOR_10 0 0.5723726
5575-GNVDE One year DSL Mailed check No No Yes No No No 03 <= 68.8 03_DE_30_A_45 0 0.0468170
3668-QPYBK Month-to-month DSL Mailed check Yes No Yes No No Yes 02 <= 55.95 01_MENOR_10 1 0.3459541
7795-CFOCW One year DSL Bank transfer (automatic) No Yes Yes No No No 02 <= 55.95 03_DE_30_A_45 0 0.0508932
9237-HQITU Month-to-month Fiber optic Electronic check Yes No No No No No 04 <= 106.75 01_MENOR_10 1 0.7750630
9305-CDSKC Month-to-month Fiber optic Electronic check Yes No No No No No 04 <= 106.75 01_MENOR_10 1 0.7750630
#Hacemos los guardados que incluyen al modelo final.

saveRDS(model_extr1,'modelo_final.rds')

saveRDS(dt1,'cache4.rds')

7.EVALUACION Y ANALISIS DE NEGOCIO.

Vamos a visualizar el abandono real por tramos de scoring. Este grafico es muy potente para ver que el modelo es consistente, ya que debe presentar una linea descendente en la tasa de abandono conforme se desciende en el scoring.

Creamos una funcion para visualizar la contratacion real por percentiles de scoring.

vis <- function(scoring,real) {
  #Preparar el dataframe de visualizacion
  vis_df <- data.frame(Scoring = scoring, Perc_Scoring = cut_number(scoring, 10), Real = real)
  levels(vis_df$Perc_Scoring) <- seq(from = 100,to = 10,by = -10)
  vis_gr <- vis_df %>% group_by(Perc_Scoring) %>% summarise(Tasa_Contr = mean(as.numeric(as.character(Real)))) %>% arrange(Perc_Scoring)
  #Ordenar el factor para el grafico
  vis_gr$Perc_Scoring <- factor(vis_gr$Perc_Scoring, levels = vis_gr$Perc_Scoring[order(vis_gr$Perc_Scoring, decreasing = T)])
  #hacemos el grafico
  ggplot(vis_gr,aes(Perc_Scoring, Tasa_Contr)) + 
    geom_col(fill='grey') + 
    geom_hline(aes(yintercept = mean(as.numeric(as.character(vis_df$Real)))), col = 'black') +
    labs(title = 'Abandono real por tramo de scoring', x = 'Tramo de Scoring', y = 'Abandono Real')
}
vis(dt1$SCORING_CHURN,dt1$Churn)

Como vemos el modelo esta bien calibrado pues se aprecia esa tendencia descendente de la tasa de abandono conforme disminuye el scoring.

Vemos tambien en la grafica que a partir del 40% el abandono es inferior a la media, a partir de este punto no nos interesaria lanzar las campañas, aunque para asegurarnos un mejor ROI, se tratara de coger tramos anteriores a ese 40%, para mayor confiabilidad de los resultados.

7.1 CASO: LANZAMIENTO DE CAMPAÑA CLIENTES TELCO.

Nuestro caso sera lanzar una campaña de fidelizacion para los clientes haciendole diversas promociones y descuentos incentivando el engagement con el cliente.

La campaña tiene un presupuesto de 23000 euros, un coste unitario de 25 euros por cliente y sera realizada por call center.

7.1.1 OPCION 1: Criterio de Negocio-presupuesto total asignado a campaña.

Haciendo las cuentas el numero de clientes a contactar sera 23000/25=920

Filtramos pues los clientes a los que vamos a lanzar la campaña para intentar fidelizarlos.

y escojemos estos 920 clientes primeros ordenados por scoring que nos han salido al ajustar el coste del cliente en funcion del tamaño de la campaña.

y visualizamos los 50 primeros.

tamaño_campaña <- 920
bote_campaña <- dt1%>%
    filter(Churn==0)%>%
    arrange(desc(SCORING_CHURN))%>%
    slice(1:tamaño_campaña)%>%
    select(customerID,SCORING_CHURN)
#Previsualizamos la salida
head(bote_campaña,50)
##     customerID SCORING_CHURN
##  1: 4482-EWFMI      0.775063
##  2: 4847-TAJYI      0.775063
##  3: 8266-VBFQL      0.775063
##  4: 2799-ARNLO      0.775063
##  5: 0021-IKXGC      0.775063
##  6: 4115-NZRKS      0.775063
##  7: 5168-MSWXT      0.775063
##  8: 1452-VOQCH      0.775063
##  9: 4234-XTNEA      0.775063
## 10: 8270-RKSAP      0.775063
## 11: 6786-OBWQR      0.775063
## 12: 2123-AGEEN      0.775063
## 13: 6630-UJZMY      0.775063
## 14: 8242-SOQUO      0.775063
## 15: 0455-XFASS      0.775063
## 16: 7439-DKZTW      0.775063
## 17: 7363-QTBIW      0.775063
## 18: 6077-BDPXA      0.775063
## 19: 4291-SHSBH      0.775063
## 20: 9957-YODKZ      0.775063
## 21: 7321-VGNKU      0.775063
## 22: 1219-NNDDO      0.775063
## 23: 6435-SRWBJ      0.775063
## 24: 2081-VEYEH      0.775063
## 25: 3030-YDNRM      0.775063
## 26: 3318-NMQXL      0.775063
## 27: 2018-QKYGT      0.775063
## 28: 1960-UYCNN      0.775063
## 29: 8087-LGYHQ      0.775063
## 30: 1393-IMKZG      0.775063
## 31: 4813-HQMGZ      0.775063
## 32: 7245-NIIWQ      0.775063
## 33: 8739-XNIKG      0.775063
## 34: 2789-HQBOU      0.775063
## 35: 9603-OAIHC      0.775063
## 36: 1475-VWVDO      0.775063
## 37: 1197-BVMVG      0.775063
## 38: 8622-ZLFKO      0.775063
## 39: 4927-WWOOZ      0.775063
## 40: 0722-SVSFK      0.775063
## 41: 8325-QRPZR      0.775063
## 42: 5542-TBBWB      0.775063
## 43: 4770-UEZOX      0.775063
## 44: 6860-YRJZP      0.775063
## 45: 5150-ITWWB      0.775063
## 46: 6818-DJXAA      0.775063
## 47: 9643-AVVWI      0.775063
## 48: 1640-PLFMP      0.775063
## 49: 0840-DFEZH      0.775063
## 50: 2545-EBUPK      0.775063
##     customerID SCORING_CHURN

Vamos a ver graficamente si de esta forma estamos aprovechando el potencial de nuestro modelo

penetracion_target <- mean(as.numeric(as.character(dt1$Churn)))
dt1 %>% 
  arrange(desc(SCORING_CHURN)) %>% 
  ggplot(aes(y = SCORING_CHURN, x = seq_along(SCORING_CHURN))) +
  geom_line() + 
  geom_vline(xintercept = tamaño_campaña, col = 'orange') +
  geom_hline(yintercept = penetracion_target,col='blue') +
  labs(x = 'CLIENTES ORDENADOS POR SCORING', y = 'SCORING')

Vemos que tenemos clientes cercanos a 3000 con un scorring mayor de la media(punto de corte linea azul hacia la izquierda), a partir de ahi hacia la derecha tendrian un scoring inferior a la media, si tomamos el triangulo que dibuja donde intersectan las rectas naranja del tamaño de la campaña(en nuestro caso 920 clientes) y azul de penetracion de la media de la target con nuestra linea del grafico de descenso del scoring tenemos el recorrido de clientes posibles a los cuales lanzar la campaña y tambien vemos que podriamos tratar de ampliar a mayor tamaño de campaña porque tendriamos margen para ello.

7.1.2 OPCION 2: Establecemos un criterio de negocio (Scoring x2 o x3)

Suele ser normal coger aquellos clientes cuyo scoring este x2 o x3 por encima de la penetracion de la target

scoring_x2 <- penetracion_target * 2
scoring_x3 <- penetracion_target * 3
bote_campaña_x2 <- dt1 %>% 
  filter(Churn==0 & SCORING_CHURN > scoring_x2) %>% 
  select(customerID,SCORING_CHURN)
#Tamaño del bote_campaña_x2
nrow(bote_campaña_x2)
## [1] 401
bote_campaña_x3 <- dt1 %>% 
  filter(Churn==0 & SCORING_CHURN > scoring_x3) %>% 
  select(customerID,SCORING_CHURN)
#Tamaño del bote_campaña_x3
nrow(bote_campaña_x3)
## [1] 0

Vamos a ver graficamente si de esta forma estamos aprovechando el potencial de nuestro modelo

dt1 %>% 
  arrange(desc(SCORING_CHURN)) %>% 
  ggplot(aes(y = SCORING_CHURN, x = seq_along(SCORING_CHURN))) +
  geom_line() + 
  geom_hline(yintercept = scoring_x2,col='blue') +
  geom_hline(yintercept = scoring_x3,col='orange') +
  geom_hline(yintercept = penetracion_target,col='black') +
  labs(x = 'CLIENTES ORDENADOS POR SCORING', y = 'SCORING')

Al hacer el estudio suponiendo que no tenemos restricciones presupuestarias para esta opcion concreta, vemos la recomendacion del numero de clientes(tamaño de campaña) que dariamos para lanzar la campaña analizando el scoring 2 veces por encima de la media de penetracion de la target (scoringx2) y el scoring 3 veces por encima de esa misma media(scoringx3), y al hacerlo vemos que tenemos 401 clientes para lanzar la campaña con resultados optimos con un nivel de scoringx2, en cambio no somos capaces de seleccionar clientes por encima de la media 3 veces(scoringx3), nos da cero(numero de clientes) el resultado, y se puede observar esto visualmente tambien en el grafico en el que por encima de la linea naranja (ni tan siquiera llega) no supera la linea de descenso del scoring esa franja.

Por lo tanto los maximos resultados optimos que obtenemos son los que ofrece el nivel de scoringx2.

7.1.3 OPCION 3: Calculamos el tamaño maximo de campaña que resultara rentable

Para ello tenemos en cuenta el ingreso medio previsto y el coste medio por accion comercial.

Asi pues calculamos el ingreso medio del cliente retenido.

clientes_retenidos <- dt%>%filter(Churn==0)

ingreso_medio_cliente <- mean(clientes_retenidos$MonthlyCharges)

ingreso_medio_cliente
## [1] 61.26512

suponiendo un margen neto de un 15%

margen_medio_anual=61.26*0.15*12

margen_medio_anual
## [1] 110.268
mar_medio <- margen_medio_anual

Supuesto de margen medio anual por cliente = 110.26€

Supuesto de coste medio por accion comercial (call center cliente contactado) = 25€

Definicion de margen esperado = probabilidad de evento * margen evento

Definicion de margen neto = margen esperado - coste medio

mar_medio <- 110.26
coste_medio <- 25
#Calculamos el margen esperado de cada cliente
df_campaña <- dt1 %>% 
  filter(Churn==0) %>% 
  mutate(
    ME = SCORING_CHURN * mar_medio,
    MN = ME - coste_medio) %>% 
  arrange(desc(MN)) %>% 
  mutate(INDICE = 1:nrow(.)) %>% 
  select(customerID,INDICE,ME,MN)
head(df_campaña,50)
##     customerID INDICE       ME       MN
##  1: 4482-EWFMI      1 85.45845 60.45845
##  2: 4847-TAJYI      2 85.45845 60.45845
##  3: 8266-VBFQL      3 85.45845 60.45845
##  4: 2799-ARNLO      4 85.45845 60.45845
##  5: 0021-IKXGC      5 85.45845 60.45845
##  6: 4115-NZRKS      6 85.45845 60.45845
##  7: 5168-MSWXT      7 85.45845 60.45845
##  8: 1452-VOQCH      8 85.45845 60.45845
##  9: 4234-XTNEA      9 85.45845 60.45845
## 10: 8270-RKSAP     10 85.45845 60.45845
## 11: 6786-OBWQR     11 85.45845 60.45845
## 12: 2123-AGEEN     12 85.45845 60.45845
## 13: 6630-UJZMY     13 85.45845 60.45845
## 14: 8242-SOQUO     14 85.45845 60.45845
## 15: 0455-XFASS     15 85.45845 60.45845
## 16: 7439-DKZTW     16 85.45845 60.45845
## 17: 7363-QTBIW     17 85.45845 60.45845
## 18: 6077-BDPXA     18 85.45845 60.45845
## 19: 4291-SHSBH     19 85.45845 60.45845
## 20: 9957-YODKZ     20 85.45845 60.45845
## 21: 7321-VGNKU     21 85.45845 60.45845
## 22: 1219-NNDDO     22 85.45845 60.45845
## 23: 6435-SRWBJ     23 85.45845 60.45845
## 24: 2081-VEYEH     24 85.45845 60.45845
## 25: 3030-YDNRM     25 85.45845 60.45845
## 26: 3318-NMQXL     26 85.45845 60.45845
## 27: 2018-QKYGT     27 85.45845 60.45845
## 28: 1960-UYCNN     28 85.45845 60.45845
## 29: 8087-LGYHQ     29 85.45845 60.45845
## 30: 1393-IMKZG     30 85.45845 60.45845
## 31: 4813-HQMGZ     31 85.45845 60.45845
## 32: 7245-NIIWQ     32 85.45845 60.45845
## 33: 8739-XNIKG     33 85.45845 60.45845
## 34: 2789-HQBOU     34 85.45845 60.45845
## 35: 9603-OAIHC     35 85.45845 60.45845
## 36: 1475-VWVDO     36 85.45845 60.45845
## 37: 1197-BVMVG     37 85.45845 60.45845
## 38: 8622-ZLFKO     38 85.45845 60.45845
## 39: 4927-WWOOZ     39 85.45845 60.45845
## 40: 0722-SVSFK     40 85.45845 60.45845
## 41: 8325-QRPZR     41 85.45845 60.45845
## 42: 5542-TBBWB     42 85.45845 60.45845
## 43: 4770-UEZOX     43 85.45845 60.45845
## 44: 6860-YRJZP     44 85.45845 60.45845
## 45: 5150-ITWWB     45 85.45845 60.45845
## 46: 6818-DJXAA     46 85.45845 60.45845
## 47: 9643-AVVWI     47 85.45845 60.45845
## 48: 1640-PLFMP     48 85.45845 60.45845
## 49: 0840-DFEZH     49 85.45845 60.45845
## 50: 2545-EBUPK     50 85.45845 60.45845
##     customerID INDICE       ME       MN

Visualizamos las curvas

Localizamos el punto en el que el margen neto pasa a ser cero o menos

MN_cero <- df_campaña %>% filter(MN <= 0 ) %>% slice(1) %>% select(INDICE)
MN_cero <- MN_cero$INDICE

Hacemos el grafico

ggplot(df_campaña,aes(x = INDICE)) + 
  geom_line(aes(y = ME, col = "ME")) + 
  geom_line(aes(y = MN, col = "MN")) + 
  geom_hline(aes(yintercept = coste_medio, col = 'COSTE MEDIO')) + 
  geom_vline(aes(xintercept = MN_cero, col = 'MARGEN NETO CERO')) +
  labs(x = 'TAMAÑO DE CAMPAÑA', y = 'MARGEN', colour = 'KPI')

print(paste('Tamaño maximo de campaña rentable: ',MN_cero))
## [1] "Tamaño maximo de campaña rentable:  1688"

Vemos que segun los parametros que nosotros le hemos definido donde la linea vertical verde se cruza hacia abajo con la azul(ME) y la rosa(MN) que coincide con el valor ‘0’ de la linea del MN(rosa) al cruzarse, a partir de ese punto (toda la linea verde que marca el equilibrio) hacia la derecha perdemos dinero, hacia la izquierda de la linea verde ganamos dinero, asi pues determinamos ese punto donde intersecta la linea verde como tamaño maximo de campaña rentable, que luego imprimimos con la sentencia print para verlo, ya que no se muestra en numero en la grafica y seria ese el tamaño maximo de la campaña que resultara rentable para nuestro caso que serian 1688 clientes. Lanzar la campaña a mas de 1688 clientes seria empezar a perder dinero.Mientras que cortar en ese punto (1688 clientes), seria lo que nos optimizaria el tamaño de campaña(punto de equilibrio).

7.1.4 OPCION 4: Calculamos el punto optimo de retorno de la inversion

Vamos a calcular 2 nuevas variables que sean un agregado de los ingresos agregados y de los gastos agregados en cada potencial tamaño de campaña, y el ROI como diferencia de las anteriores, y vamos a localizar el tamaño de la campaña que va a maximizar ese ROI y tambien cuanto vamos a ganar previsiblemente

Vamos a usar la funcion cumsum(), que hace una suma acumulada secuencial de la variable que le pases como parametro

df_campaña <- df_campaña %>% 
  mutate(
    INGRESOS_AGRE = cumsum(ME),
    COSTES_AGRE = INDICE * coste_medio,
    ROI = INGRESOS_AGRE - COSTES_AGRE)
head(df_campaña,50)
##     customerID INDICE       ME       MN INGRESOS_AGRE COSTES_AGRE        ROI
##  1: 4482-EWFMI      1 85.45845 60.45845      85.45845          25   60.45845
##  2: 4847-TAJYI      2 85.45845 60.45845     170.91689          50  120.91689
##  3: 8266-VBFQL      3 85.45845 60.45845     256.37534          75  181.37534
##  4: 2799-ARNLO      4 85.45845 60.45845     341.83378         100  241.83378
##  5: 0021-IKXGC      5 85.45845 60.45845     427.29223         125  302.29223
##  6: 4115-NZRKS      6 85.45845 60.45845     512.75067         150  362.75067
##  7: 5168-MSWXT      7 85.45845 60.45845     598.20912         175  423.20912
##  8: 1452-VOQCH      8 85.45845 60.45845     683.66756         200  483.66756
##  9: 4234-XTNEA      9 85.45845 60.45845     769.12601         225  544.12601
## 10: 8270-RKSAP     10 85.45845 60.45845     854.58445         250  604.58445
## 11: 6786-OBWQR     11 85.45845 60.45845     940.04290         275  665.04290
## 12: 2123-AGEEN     12 85.45845 60.45845    1025.50135         300  725.50135
## 13: 6630-UJZMY     13 85.45845 60.45845    1110.95979         325  785.95979
## 14: 8242-SOQUO     14 85.45845 60.45845    1196.41824         350  846.41824
## 15: 0455-XFASS     15 85.45845 60.45845    1281.87668         375  906.87668
## 16: 7439-DKZTW     16 85.45845 60.45845    1367.33513         400  967.33513
## 17: 7363-QTBIW     17 85.45845 60.45845    1452.79357         425 1027.79357
## 18: 6077-BDPXA     18 85.45845 60.45845    1538.25202         450 1088.25202
## 19: 4291-SHSBH     19 85.45845 60.45845    1623.71046         475 1148.71046
## 20: 9957-YODKZ     20 85.45845 60.45845    1709.16891         500 1209.16891
## 21: 7321-VGNKU     21 85.45845 60.45845    1794.62736         525 1269.62736
## 22: 1219-NNDDO     22 85.45845 60.45845    1880.08580         550 1330.08580
## 23: 6435-SRWBJ     23 85.45845 60.45845    1965.54425         575 1390.54425
## 24: 2081-VEYEH     24 85.45845 60.45845    2051.00269         600 1451.00269
## 25: 3030-YDNRM     25 85.45845 60.45845    2136.46114         625 1511.46114
## 26: 3318-NMQXL     26 85.45845 60.45845    2221.91958         650 1571.91958
## 27: 2018-QKYGT     27 85.45845 60.45845    2307.37803         675 1632.37803
## 28: 1960-UYCNN     28 85.45845 60.45845    2392.83647         700 1692.83647
## 29: 8087-LGYHQ     29 85.45845 60.45845    2478.29492         725 1753.29492
## 30: 1393-IMKZG     30 85.45845 60.45845    2563.75336         750 1813.75336
## 31: 4813-HQMGZ     31 85.45845 60.45845    2649.21181         775 1874.21181
## 32: 7245-NIIWQ     32 85.45845 60.45845    2734.67026         800 1934.67026
## 33: 8739-XNIKG     33 85.45845 60.45845    2820.12870         825 1995.12870
## 34: 2789-HQBOU     34 85.45845 60.45845    2905.58715         850 2055.58715
## 35: 9603-OAIHC     35 85.45845 60.45845    2991.04559         875 2116.04559
## 36: 1475-VWVDO     36 85.45845 60.45845    3076.50404         900 2176.50404
## 37: 1197-BVMVG     37 85.45845 60.45845    3161.96248         925 2236.96248
## 38: 8622-ZLFKO     38 85.45845 60.45845    3247.42093         950 2297.42093
## 39: 4927-WWOOZ     39 85.45845 60.45845    3332.87937         975 2357.87937
## 40: 0722-SVSFK     40 85.45845 60.45845    3418.33782        1000 2418.33782
## 41: 8325-QRPZR     41 85.45845 60.45845    3503.79626        1025 2478.79626
## 42: 5542-TBBWB     42 85.45845 60.45845    3589.25471        1050 2539.25471
## 43: 4770-UEZOX     43 85.45845 60.45845    3674.71316        1075 2599.71316
## 44: 6860-YRJZP     44 85.45845 60.45845    3760.17160        1100 2660.17160
## 45: 5150-ITWWB     45 85.45845 60.45845    3845.63005        1125 2720.63005
## 46: 6818-DJXAA     46 85.45845 60.45845    3931.08849        1150 2781.08849
## 47: 9643-AVVWI     47 85.45845 60.45845    4016.54694        1175 2841.54694
## 48: 1640-PLFMP     48 85.45845 60.45845    4102.00538        1200 2902.00538
## 49: 0840-DFEZH     49 85.45845 60.45845    4187.46383        1225 2962.46383
## 50: 2545-EBUPK     50 85.45845 60.45845    4272.92227        1250 3022.92227
##     customerID INDICE       ME       MN INGRESOS_AGRE COSTES_AGRE        ROI

Visualizamos las curvas

ggplot(df_campaña,aes(x = INDICE)) +
  geom_line(aes(y = INGRESOS_AGRE, col='INGRESOS_AGRE')) + 
  geom_line(aes(y = COSTES_AGRE, col='COSTES_AGRE')) +
  geom_line(aes(y = ROI, col='ROI')) + 
  labs(y='EUROS', x = 'TAMAÑO DE CAMPAÑA')

Lo que hago a posteriori es amplificar la visualizacion de la linea del ROI azul hasta el punto de la curva en el que el ROI=0 a partir de ahi hacia la derecha cada accion comercial que realicemos estaremos perdiendo dinero. Con ello pretendemos dilucidar el pto optimo de la curva (que va a coincidir con el hallado en la opcion 3), la diferencia es que de forma añadida podremos entre otros datos sacar el ROI de la campaña en ese punto optimo que es lo que estamos buscando y no tenemos.

Vamos por tanto a visualizar un zoom sobre el ROI solo en los tamaños de campaña que son positivos para localizar el punto optimo

df_campaña %>% 
  filter(ROI > 0) %>% 
  ggplot(aes(x = INDICE)) +
  geom_line(aes(y = ROI, col='ROI')) +
  geom_vline(aes(xintercept = MN_cero, col = 'PUNTO OPTIMO')) +
  labs(x = 'TAMAÑO DE CAMPAÑA', y = 'ROI', colour = 'KPI')

Esta curva(linea azul del ROI)interceptara con la linea roja vertical que corresponde al MN_cero para darnos el punto optimo descrito.

Imprimimos un pequeño informe donde nos dara el resultado final con los datos totales requeridos de forma resumida.

cat(
  paste(
    'El tamaño optimo de campaña para el ROI es de:', MN_cero, 'clientes',
    '\nCon unos ingresos esperados de margen neto acumulado de:',   round(df_campaña[which(df_campaña$INDICE == MN_cero),'INGRESOS_AGRE']), '€',
    '\nY unos costes agregados de:',
    df_campaña[which(df_campaña$INDICE == MN_cero),'COSTES_AGRE'], '€',
    '\nQue van a generar un Retorno Neto de la Inversion de:',
    round(df_campaña[which(df_campaña$INDICE == MN_cero),'ROI']),'€'
  )
)
## El tamaño optimo de campaña para el ROI es de: 1688 clientes 
## Con unos ingresos esperados de margen neto acumulado de: 78985 € 
## Y unos costes agregados de: 42200 € 
## Que van a generar un Retorno Neto de la Inversion de: 36785 €