Problema de análisis

Este análisis se centra en el comportamiento de los clientes de telecomunicaciones que tienen más probabilidades de abandonar la plataforma. El objetivo radica en descubrir el comportamiento más sorprendente de los clientes a través de EDA y luego usar algunas de las técnicas de análisis predictivo para determinar los clientes que tienen más probabilidades de abandonar.

Sobre los datos

El dataset comprende información sobre:

  • Clientes que se fueron en el último mes: la columna se llama Renuncia
  • Servicios a los que se ha suscrito cada cliente: teléfono, varias líneas, Internet, seguridad en línea, respaldo en línea, protección de dispositivos, soporte técnico y transmisión de TV y películas
  • Información de la cuenta del cliente: cuánto tiempo ha sido cliente, contrato, método de pago, facturación electrónica, cargos mensuales y cargos totales.
  • Información demográfica sobre los clientes: sexo, rango de edad y si tienen socios y dependientes

Análisis descriptivo

Se muestran a continuación las primeras seis filas del conjunto de datos:

customerID gender SeniorCitizen Partner Dependents tenure PhoneService MultipleLines InternetService OnlineSecurity OnlineBackup DeviceProtection TechSupport StreamingTV StreamingMovies Contract PaperlessBilling PaymentMethod MonthlyCharges TotalCharges Churn
7590-VHVEG Female 0 Yes No 1 No No phone service DSL No Yes No No No No Month-to-month Yes Electronic check 29.85 29.85 No
5575-GNVDE Male 0 No No 34 Yes No DSL Yes No Yes No No No One year No Mailed check 56.95 1889.50 No
3668-QPYBK Male 0 No No 2 Yes No DSL Yes Yes No No No No Month-to-month Yes Mailed check 53.85 108.15 Yes
7795-CFOCW Male 0 No No 45 No No phone service DSL Yes No Yes Yes No No One year No Bank transfer (automatic) 42.30 1840.75 No
9237-HQITU Female 0 No No 2 Yes No Fiber optic No No No No No No Month-to-month Yes Electronic check 70.70 151.65 Yes
9305-CDSKC Female 0 No No 8 Yes Yes Fiber optic No No Yes No Yes Yes Month-to-month Yes Electronic check 99.65 820.50 Yes

Hay 7043 filas de datos para 21 variables de clientes de telecomunicaciones.

Limpieza del dataset

Se verifica la existencia de missing cases en las distintas variables:

options(repr.plot.width = 6, repr.plot.height = 4)
missing_data <- telco %>% summarise_all(funs(sum(is.na(.))/n()))
missing_data <- gather(missing_data, key = "variables", value = "percent_missing")

#Gráfico
ggplot(missing_data, aes(x = reorder(variables, percent_missing), y = percent_missing)) +
  geom_bar(stat = "identity", fill = "#5817F5", aes(color = I('white')), size = 0.3)+
  xlab('Variables')+
  ylab('Valores perdidos') + 
  coord_flip()+ 
  theme_bw() +
  theme(plot.title = element_text(size=12, face = "bold.italic", color = "#2e2959"),
                          axis.text.y = element_text(size = 8, color = "#2e2959"),
                          axis.text.x = element_text(size = 8, color = "#2e2959"),
                          axis.title.x = element_text(color = "#2e2959"),
                          axis.title.y = element_text(color = "#2e2959"),
                          plot.caption = element_text(color = "#2e2959"),
                          legend.text = element_text(color = "#2e2959"),
                          legend.title = element_text(color = "#2e2959"),
                          plot.subtitle = element_text(color = "#2e2959"))

Se encuentra que en la variable TotalCharges existen filas sin casos, por lo que se suprimen dichas observaciones.

telco <- telco[complete.cases(telco),]

También se modifica la escala de medición de la variable SeniorCitizen (de factor, a interger):

telco$SeniorCitizen <- as.factor(ifelse(telco$SeniorCitizen==1, 'YES', 'NO'))

Para hacer el modelo dicriminante, se usarán solo variables cuantitativas (tenure, MonthlyCharges, TotalCharges) como variables independientes y variable “churn (renuncia)” como variable dependiente.

telco_reducido <- subset(telco, select = c(tenure, MonthlyCharges, TotalCharges, Churn))

Para evaluar los resultados, se divide el conjunto de datos en conjunto de entrenamiento y conjunto de prueba:

set.seed(100)
test.index <- sample(1:nrow(telco_reducido), 1000)

#Conjunto de entrenamiento
training_set <- telco_reducido[-test.index,]

#conjunto de prueba
test_set <- telco_reducido[test.index,]

Se realizan modificaciones sobre el conjunto de datos de entrenamiento:

training_set_shapiro <- training_set

#Se modifica la variable dependiente (VD)
training_set_shapiro$Renuncia <- ifelse(training_set_shapiro$Churn=="Yes", "Renuncia", "No Renuncia")

#Se suprime la VD en su formato anterior
training_set_shapiro <- subset( training_set_shapiro, select = -Churn )

#Se toma una muestra 
training_set_shapiro <- sample_n(training_set_shapiro, 1500)

kable_styling(kable(head(training_set_shapiro)))
tenure MonthlyCharges TotalCharges Renuncia
52 89.45 4577.75 No Renuncia
15 19.75 311.60 No Renuncia
11 82.90 880.05 No Renuncia
15 80.20 1217.25 Renuncia
7 69.55 521.35 No Renuncia
59 96.65 5580.80 No Renuncia

También se modifica la variable dependiente para volverla factorial:

training_set_shapiro$Renuncia <- as.factor(training_set_shapiro$Renuncia)
class(training_set_shapiro$Renuncia)
## [1] "factor"

Distribución de variables

Habiendo preparado los datasets, se procede a indagar la distribución en cada variable:

options(repr.plot.width =6, repr.plot.height = 2)
#Gráfico 
ggplot(training_set_shapiro, aes(y= tenure, x = "", fill = Renuncia)) + 
geom_boxplot()+ 
theme_bw()+
xlab(" ")+
  theme(plot.title = element_text(size=12, face = "bold.italic", color = "#2e2959"),
                          axis.text.y = element_text(size = 8, color = "#2e2959"),
                          axis.text.x = element_text(size = 8, color = "#2e2959"),
                          axis.title.x = element_text(color = "#2e2959"),
                          axis.title.y = element_text(color = "#2e2959"),
                          plot.caption = element_text(color = "#2e2959"),
                          legend.text = element_text(color = "#2e2959"),
                          legend.title = element_text(color = "#2e2959"),
                          plot.subtitle = element_text(color = "#2e2959")) +
  scale_fill_brewer("Renuncia", palette = "Paired")

Se observa que la mediana de la antiguedad es mayor en el caso de los clientes que no han renunciado al servicio.

Cargas mensuales

ggplot(training_set_shapiro, aes(y= MonthlyCharges, x = "", fill = Renuncia)) + 
geom_boxplot()+ 
theme_bw()+
xlab(" ")+
  theme(plot.title = element_text(size=12, face = "bold.italic", color = "#2e2959"),
                          axis.text.y = element_text(size = 8, color = "#2e2959"),
                          axis.text.x = element_text(size = 8, color = "#2e2959"),
                          axis.title.x = element_text(color = "#2e2959"),
                          axis.title.y = element_text(color = "#2e2959"),
                          plot.caption = element_text(color = "#2e2959"),
                          legend.text = element_text(color = "#2e2959"),
                          legend.title = element_text(color = "#2e2959"),
                          plot.subtitle = element_text(color = "#2e2959")) +
  scale_fill_brewer("Renuncia", palette = "Paired")

Los clientes que abandonaron el servicio tienen una mediana más alta en la carga de servicios mensuales.

Cargas totales

ggplot(training_set_shapiro, aes(y= TotalCharges, x = "", fill = Renuncia)) + 
geom_boxplot()+ 
theme_bw()+
xlab(" ")+
  theme(plot.title = element_text(size=12, face = "bold.italic", color = "#2e2959"),
                          axis.text.y = element_text(size = 8, color = "#2e2959"),
                          axis.text.x = element_text(size = 8, color = "#2e2959"),
                          axis.title.x = element_text(color = "#2e2959"),
                          axis.title.y = element_text(color = "#2e2959"),
                          plot.caption = element_text(color = "#2e2959"),
                          legend.text = element_text(color = "#2e2959"),
                          legend.title = element_text(color = "#2e2959"),
                          plot.subtitle = element_text(color = "#2e2959")) +
  scale_fill_brewer("Renuncia", palette = "Paired")

Los clientes que renunciaron al servicio posee una mediana menor en sus cargas totales, respecto a los que no renunciaron.

Luego, se realiza una representación mediante histogramas de cada variable para cada categoría de la VD:

par(mfcol = c(2, 3))
for (k in 1:3) {
  j0 <- names(training_set_shapiro)[k]
  #br0 <- seq(min(datos[, k]), max(datos[, k]), le = 11)
  x0 <- seq(min(training_set_shapiro[, k]), max(training_set_shapiro[, k]), le = 50)
  for (i in 1:2) {
    i0 <- levels(training_set_shapiro$Renuncia)[i]
    x <- training_set_shapiro[training_set_shapiro$Renuncia == i0, j0]
    hist(x, proba = T, col = grey(0.8), main = paste("Caso ", i0),
    xlab = j0)
    lines(x0, dnorm(x0, mean(x), sd(x)), col = "#5817F5", lwd = 2)
  }
}

Se muestra lo equivalente mediante otra visualización:

plot1 <- ggplot(data = training_set_shapiro, aes(x = tenure)) +
         geom_density(aes(colour = Renuncia)) + theme_bw() +
  theme(plot.title = element_text(size=12, face = "bold.italic", color = "#2e2959"),
                          axis.text.y = element_text(size = 8, color = "#2e2959"),
                          axis.text.x = element_text(size = 8, color = "#2e2959"),
                          axis.title.x = element_text(color = "#2e2959"),
                          axis.title.y = element_text(color = "#2e2959"),
                          plot.caption = element_text(color = "#2e2959"),
                          legend.text = element_text(color = "#2e2959"),
                          legend.title = element_text(color = "#2e2959"),
                          plot.subtitle = element_text(color = "#2e2959")) +
  scale_color_brewer("Renuncia", palette = "Paired")
plot2 <- ggplot(data = training_set_shapiro, aes(x = MonthlyCharges)) +
         geom_density(aes(colour = Renuncia)) + theme_bw() +
  theme(plot.title = element_text(size=12, face = "bold.italic", color = "#2e2959"),
                          axis.text.y = element_text(size = 8, color = "#2e2959"),
                          axis.text.x = element_text(size = 8, color = "#2e2959"),
                          axis.title.x = element_text(color = "#2e2959"),
                          axis.title.y = element_text(color = "#2e2959"),
                          plot.caption = element_text(color = "#2e2959"),
                          legend.text = element_text(color = "#2e2959"),
                          legend.title = element_text(color = "#2e2959"),
                          plot.subtitle = element_text(color = "#2e2959")) +
  scale_color_brewer("Renuncia", palette = "Paired")
plot3 <- ggplot(data = training_set_shapiro, aes(x = TotalCharges)) +
         geom_density(aes(colour = Renuncia)) + theme_bw() +
  theme(plot.title = element_text(size=12, face = "bold.italic", color = "#2e2959"),
                          axis.text.y = element_text(size = 8, color = "#2e2959"),
                          axis.text.x = element_text(size = 8, color = "#2e2959"),
                          axis.title.x = element_text(color = "#2e2959"),
                          axis.title.y = element_text(color = "#2e2959"),
                          plot.caption = element_text(color = "#2e2959"),
                          legend.text = element_text(color = "#2e2959"),
                          legend.title = element_text(color = "#2e2959"),
                          plot.subtitle = element_text(color = "#2e2959")) +
  scale_color_brewer("Renuncia", palette = "Paired")

ggarrange(plot1, plot2, plot3, common.legend = TRUE, legend = "bottom")

Supuestos

Normalidad

Se realiza el contraste de normalidad Shapiro-Wilk para cada variable, según la categoría de abandono

library(reshape2)
library(knitr)
library(dplyr)
datos_tidy <- melt(training_set_shapiro, value.name = "valor")
kable(datos_tidy %>% group_by(Renuncia, variable) %>% summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,10)))
Renuncia variable p_value_Shapiro.test
No Renuncia tenure 0
No Renuncia MonthlyCharges 0
No Renuncia TotalCharges 0
Renuncia tenure 0
Renuncia MonthlyCharges 0
Renuncia TotalCharges 0

Se verifica en todos los casos un p-value < 0.05.

Homogeneidad de varianzas

Se procede a verificar la homogeneidad de varianzas mediante el test de Bartlett.

options(scipen = 999)

bartlett.test(training_set_shapiro$tenure~training_set_shapiro$Renuncia)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  training_set_shapiro$tenure by training_set_shapiro$Renuncia
## Bartlett's K-squared = 19.146, df = 1, p-value = 0.00001211
bartlett.test(training_set_shapiro$MonthlyCharges~training_set_shapiro$Renuncia)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  training_set_shapiro$MonthlyCharges by training_set_shapiro$Renuncia
## Bartlett's K-squared = 25.619, df = 1, p-value = 0.000000416
bartlett.test(training_set_shapiro$TotalCharges~training_set_shapiro$Renuncia)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  training_set_shapiro$TotalCharges by training_set_shapiro$Renuncia
## Bartlett's K-squared = 25.418, df = 1, p-value = 0.0000004616

Se observa que en todos los casos el p-value < 0,05, por lo que se rechaza la H0 que afirma igualdad de varianzas.

Test M de Box

La hipótesis nula para esta prueba es que las matrices de covarianza observadas para las variables dependientes son iguales en todos los grupos. En otras palabras, un resultado de prueba no significativo (es decir, uno con un valor p grande) indica que las matrices de covarianza son iguales. La estadística de prueba generada se llama estadística M de Box.

Esta prueba informará un resultado estadísticamente significativo cuando en realidad no exista uno.

library(biotools)
boxM(data=training_set_shapiro[, 1:3], group=training_set_shapiro$Renuncia)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  training_set_shapiro[, 1:3]
## Chi-Sq (approx.) = 374.73, df = 6, p-value < 0.00000000000000022

Se rechaza la H0.

Modelo

El primer resultado es un modelo discriminante lineal que puede predecir qué personas se irían usando las variables independientes. Después de hacer este modelo, se calcula la relación de clasificación correcta.

#linear discrimanant model
lda1 <- lda(formula = Churn ~ tenure + MonthlyCharges + TotalCharges, data= training_set)

#training (0,78)
predicted_lda1_train <- predict(lda1, newdata = training_set)
table(training_set$Churn, predicted_lda1_train$class)
##      
##         No  Yes
##   No  3963  473
##   Yes  858  738
cat("Accuracy training: ",sum(diag(as.matrix(table(training_set$Churn, predicted_lda1_train$class)))) / sum(sum(as.matrix(table(training_set$Churn, predicted_lda1_train$class)))))
## Accuracy training:  0.7793435
#test (0,78)
predicted_lda1 <- predict(lda1, newdata = test_set)
table(test_set$Churn, predicted_lda1$class)
##      
##        No Yes
##   No  636  91
##   Yes 133 140
cat("Accuracy testing: ",sum(diag(as.matrix(table(test_set$Churn, predicted_lda1$class)))) / sum(sum(as.matrix(table(test_set$Churn, predicted_lda1$class)))))
## Accuracy testing:  0.776

Salida del modelo de discriminante lineal:

lda1
## Call:
## lda(Churn ~ tenure + MonthlyCharges + TotalCharges, data = training_set)
## 
## Prior probabilities of groups:
##        No       Yes 
## 0.7354111 0.2645889 
## 
## Group means:
##       tenure MonthlyCharges TotalCharges
## No  37.98873       61.36540     2577.293
## Yes 18.04950       74.24721     1538.819
## 
## Coefficients of linear discriminants:
##                          LD1
## tenure         -0.0161921776
## MonthlyCharges  0.0379090498
## TotalCharges   -0.0003956721

Se dibujan las áreas que se delimitan para cada una de las variables:

library(klaR)
partimat(Renuncia ~ tenure + MonthlyCharges + TotalCharges,
         data= training_set_shapiro, method = "lda", prec = 10,
         image.colors = c("darkgoldenrod1", "skyblue2"),
         col.mean = "firebrick")

Comparación de modelo lineal y modelo no lineal

El modelo discriminante lineal muestra aproximadamente un 78% de precisión. El siguiente modelo es el modelo discriminatorio no lineal.

#nonliner discrimanant model
qda1 <-  qda(formula = Churn ~ tenure + MonthlyCharges + TotalCharges, data= training_set)

#trainign (74%)
predicted_qda1_train <- predict(qda1, newdata = training_set)
table(training_set$Churn, predicted_qda1_train$class)
##      
##         No  Yes
##   No  3455  981
##   Yes  563 1033
cat("Accuracy training: ",sum(diag(as.matrix(table(training_set$Churn, predicted_qda1_train$class)))) / sum(sum(as.matrix(table(training_set$Churn, predicted_qda1_train$class)))))
## Accuracy training:  0.7440318
#test (72%)
predicted_qda1 <-  predict(qda1, newdata = test_set)
table(test_set$Churn, predicted_qda1$class)
##      
##        No Yes
##   No  541 186
##   Yes  93 180
#the ratio of correct ratio of correct classification.
cat("Accuracy training: ",sum(diag(as.matrix(table(test_set$Churn, predicted_qda1$class)))) / sum(sum(as.matrix(table(test_set$Churn, predicted_qda1$class)))))
## Accuracy training:  0.721

Salida del modelo de discriminante no lineal:

qda1
## Call:
## qda(Churn ~ tenure + MonthlyCharges + TotalCharges, data = training_set)
## 
## Prior probabilities of groups:
##        No       Yes 
## 0.7354111 0.2645889 
## 
## Group means:
##       tenure MonthlyCharges TotalCharges
## No  37.98873       61.36540     2577.293
## Yes 18.04950       74.24721     1538.819

Se dibujan las áreas que se delimitan para cada una de las variables:

partimat(formula = Renuncia ~ tenure + MonthlyCharges + TotalCharges,
         data= training_set_shapiro,
         method = "qda", prec = 5,
         image.colors = c("darkgoldenrod1", "skyblue2"),
         col.mean = "firebrick")

El modelo discriminante no lineal muestra aproximadamente el 72% de precisión. Es mejor la predicción del modelo discriminante lineal.

CURVA ROC

Por último, se analiza la comparación de modelos mediante una curva ROC para cada modelo.

#ROC curve
#linear discrimininant
par(mfrow = c(2,1))
lda1_p <- predict(lda1, newdata = test_set, type = "response")
lda1_pr <- prediction(as.numeric(lda1_p$class), test_set$Churn)
lda1_prf <- performance(lda1_pr, measure = "tpr", x.measure = "fpr")
plot(lda1_prf, main = 'ROC curve of linear discriminant model')

lda1_auc <- performance(lda1_pr, measure = "auc")
lda1_auc <- lda1_auc@y.values
cat("AUC Discriminante Lineal: ",round(as.numeric(lda1_auc), digits = 4))
## AUC Discriminante Lineal:  0.6938
#nonlinear discriminant model
qda1_p <- predict(qda1, newdata = test_set, type = "response")
qda1_pr <- prediction(as.numeric(qda1_p$class), test_set$Churn)
qda1_prf <- performance(qda1_pr, measure = "tpr", x.measure = "fpr")
plot(qda1_prf, main = 'ROC curve of nonlinear discriminant model')

qda1_auc <- performance(qda1_pr, measure = "auc")
qda1_auc <- qda1_auc@y.values
cat("Discriminante no lineal: ",round(as.numeric(qda1_auc), digits = 4))
## Discriminante no lineal:  0.7017

Resultado de comparación

Desde la curva ROC y AUC (Área bajo la curva), el AUC de los modelos es similar (Lineal = 0.69; No lineal = 0.70) lo que indica que el modelo es aceptable.

Predicción de nuevos clientes

Finalmente, se prueba realizar una estimación para un nuevo cliente para saber si va a abandonar los servicios o no.

Prueba 1

nuevas_observaciones <- data.frame(tenure = 10, MonthlyCharges = 57, TotalCharges = 1800)

predict(object = lda1, newdata = nuevas_observaciones)
## $class
## [1] No
## Levels: No Yes
## 
## $posterior
##         No      Yes
## 1 0.739147 0.260853
## 
## $x
##         LD1
## 1 0.2719118

Prueba 2

nuevas_observaciones <- data.frame(tenure = 50, MonthlyCharges = 127, TotalCharges = 2800)

predict(object = lda1, newdata = nuevas_observaciones)
## $class
## [1] Yes
## Levels: No Yes
## 
## $posterior
##          No       Yes
## 1 0.2836631 0.7163369
## 
## $x
##        LD1
## 1 1.882186

Conclusión

Se verifica que ante más antiguedad y más cargas mensuales, más probable es la predicción de abandono del servicio. Si bien las cargas totales siguen una tendencia alejada de la mediana observada en cada categoría, cabe considerar que este es el coeficiente con menor peso en el modelo. Estos resultados son similares a los obtenidos mediante la aplicación del método del regresión logística.