O que é um churn?
Podemos definir brevemente o churn do cliente (mais comumente chamado de “churn”) como clientes que param de fazer negócios com uma empresa ou um serviço. Existem problemas de clientes em diferentes áreas de negócios e neste estudo vamos nos concentrar na área de telecomunicações. Aqui, queremos prever quais clientes deixarão seu provedor de telecomunicações atual. Esse tipo de estudo é chamado de Análise de Churn (Churn Analysis). É um tópico moderno em departamentos de gerenciamento de relacionamento com os clientes (CRM - Customer Relationship Management) porque é mais custoso encontrar novos clientes do que manter os existentes. Logo, as empresas querem melhorar a retenção de seus clientes.
Como identificar clientes mais propensos a sair?
A identificação é realizada por meio de um banco de dados sobre os clientes que já saíram do negócio. Com esses dados, é possível desenvolver um modelo que identifique os clientes que possuem um perfil semelhante dos que já deixaram de usar o serviço.
Para simular um experimento em que queremos prever se nossos clientes vão dar “churn”, precisamos trabalhar com um banco de dados particionado, seguindo a metodologia de validação cruzada (cross validation). O banco de dados tem duas partes, uma será o conjunto de treinamento que será usado para criar o modelo, e a segunda parte será o conjunto de testes que será usado para avaliar o nosso modelo.
Nesse caso, conhecidas as respostas dos clientes do conjunto de teste, poderemos comparar a previsão do modelo com as respostas verdadeiras. No entanto, na realidade, não sabemos quais serão as respostas verdadeiras. Portanto, devemos segmentar principalmente clientes com alta probabilidade de churn. Essa probabilidade será dada pelo nosso modelo.
A avaliação do modelo, com base nas previsões no conjunto de testes, é calculada pela:
Taxa de acerto (%): número de previsões classificadas corretamente/total da amostra. Em valores brutos conhecidos como verdadeiros positivos;
Taxa de erro (%): número de previsões classificadas incorretamente/total da amostra. Em valores brutos conhecidos como falsos positivos.
Obs: As taxas são referentes às previsões das observações de interesse, em nosso caso, quando o cliente dá “churn” no serviço.
Parte 1 - Seleção do modelo
Primeiros passos
Importamos e visualizamos rapidamente os dados.
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse")
# Importe a base da IBM através do site mencionado no começo do material
path = "WA_Fn-UseC_-Telco-Customer-Churn.csv"
db_churn = read_csv(path)Ou você pode ler direto do Github
path = "https://raw.githubusercontent.com/treselle-systems/customer_churn_analysis/master/WA_Fn-UseC_-Telco-Customer-Churn.csv"
db_churn = read_csv(path)
O conjunto de dados inclui informações sobre:
Clientes (Customers) que saíram no último mês (“Churn” )
Serviços contratados pelos clientes - telefone (“PhoneService” ), multiplas linhas(“MultipleLines”), internet (“InternetService”), online security (“OnlineSecurity”), online backup (“OnlineBackup”), proteção do aparelho (“DeviceProtection”), tech support (“TechSupport” ), streaming TV (“StreamingTV”) e streaming movies (“StreamingMovies” ).
Informações sobre as contas dos clientes - Há quanto tempo são clientes (“Contract” ), método de pagamento (“PaymentMethod”), faturamento sem papel (“PaperlessBilling”), Cobranças mensais (“MonthlyCharges”), Custos totais (“TotalCharges”), tempo que a pessoa é cliente da operadora, em meses - customer life (“tenure”)2.
Informações demográficas: gênero (“gender” ), se tem parceiros (“Partner”), se tem dependentes ( “Dependents” ), idade (“SeniorCitizen”).
A variável SeniorCitizen está classificada como inteiro, assumindo valores 1 para quem é idoso e 0 caso contrário. Assim, codificamos-a como uma variável classificada como fator.
# Mudando a classe da variável para fator
db_churn$SeniorCitizen = as.factor(db_churn$SeniorCitizen)
summary(db_churn)## customerID gender SeniorCitizen Partner
## Length:7043 Length:7043 0:5901 Length:7043
## Class :character Class :character 1:1142 Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## Dependents tenure PhoneService MultipleLines
## Length:7043 Min. : 0.00 Length:7043 Length:7043
## Class :character 1st Qu.: 9.00 Class :character Class :character
## Mode :character Median :29.00 Mode :character Mode :character
## Mean :32.37
## 3rd Qu.:55.00
## Max. :72.00
##
## InternetService OnlineSecurity OnlineBackup
## Length:7043 Length:7043 Length:7043
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## DeviceProtection TechSupport StreamingTV
## Length:7043 Length:7043 Length:7043
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## StreamingMovies Contract PaperlessBilling
## Length:7043 Length:7043 Length:7043
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## PaymentMethod MonthlyCharges TotalCharges Churn
## Length:7043 Min. : 18.25 Min. : 18.8 Length:7043
## Class :character 1st Qu.: 35.50 1st Qu.: 401.4 Class :character
## Mode :character Median : 70.35 Median :1397.5 Mode :character
## Mean : 64.76 Mean :2283.3
## 3rd Qu.: 89.85 3rd Qu.:3794.7
## Max. :118.75 Max. :8684.8
## NA's :11
Note que após a codificação da variável SeniorCitizen e do uso da função summary(), os dados apresentaram 3 variáveis númericas e 18 variáveis categóricas.
Agora, observemos a taxa de churn na base de dados.
if (!require("janitor")) install.packages("janitor")
library("janitor")
prop = tabyl(db_churn$Churn)
prop
Percebe-se que temos apenas 11 valores faltantes (missing) em nosso conjunto de dados e todos estão vinculados à variável TotalCharges (total de taxas) e nenhum deles está com a variável Churn igual a não, ou seja, são clientes que permanceram no negócio. Eles representam apenas 0,16% de nossas observações totais. Então, decidimos removê-los.
library(dplyr)
# Filtrando os missings
n_NA = db_churn %>%
filter(is.na(TotalCharges)) %>%
select(Churn)
# Porcentagem de valores faltantes (NA) na base de dados
100*nrow(n_NA)/nrow(db_churn)## [1] 0.1561834# 11 NA, 0,16% do nosso banco de dados e nennhum deles representa churn positivo
# Removendo as observações faltantes
db_churn = db_churn %>%
filter(!is.na(TotalCharges))
Agora podemos dividir o banco de dados em duas partes, uma para treinamento e outra para teste.
if (!require("caret")) install.packages("caret")
library("caret")
set.seed(7)
trainId = createDataPartition(db_churn$Churn,
p = 0.7, list = FALSE, times = 1)
# Dados de treinamento e teste respectivamente
db_train = db_churn[trainId,]
db_test = db_churn[-trainId,]
Como observamos anteriormente (summary()), existem 3 variáveis numéricas e que não possuem a mesma escala nas observações. Portanto, precisamos usar um método para redimensioná-las. Primeiro verificamos a distribuição dessas variáveis.
gather_train = gather(db_train %>%
select(customerID, MonthlyCharges,TotalCharges, tenure),
variable, value, -customerID)
ggplot(gather_train , aes(value)) + facet_wrap(~variable, scales = 'free_x') +
geom_histogram()
Destaca-se que nenhuma das variáveis no gráfico anterior se assemelha a distribuição normal, portanto, redimensionamos-as sem normalização. Para isso, foi criada uma função de padronização.
normalize = function(x) {
result = (x - min(x, na.rm = TRUE)
) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
return(result)
}
#norm.train = lapply(db_train %>%
# select(MonthlyCharges, TotalCharges, tenure), normalize)
#norm.train = do.call(cbind, norm.train) %>%
# as.data.frame()
norm.train = sapply(db_train %>%
select(MonthlyCharges, TotalCharges, tenure), normalize)
Ao visualizar os dados é possível identificar que algumas variáveis têm os seguintes fatores:
“Sim”;
“Não”;
“Nenhum serviço de internet” (ou “nenhum serviço de telefone”).
O terceiro fator não nos dá nenhuma informação, então recodificamos-o para “Não”.
factor.train = lapply(db_train %>%
select(-customerID,-MonthlyCharges,
-TotalCharges, -tenure),
function(x){
x = gsub("No internet service", "No", x)
x = gsub("No phone service", "No", x)
return(x)
})
factor.train = do.call(cbind, factor.train) %>%
as.data.frame()
db_train = cbind(customerID = db_train[,1],
factor.train, norm.train)
Agora que o conjunto de treinamento foi organizado e limpo, podemos estimar nosso train com 3 modelos diferentes:
1. Decision Tree
2. Random Forest
3. Logistic RegressionModelagem
Muitas são as opções de modelos para este tipo de problema e para o nosso caso foram escolhidos três dos modelos mais usados, são eles: Decision Tree, Random Florest e Logistic Regression. Ao final os modelos serão comparados entre si e avaliado qual apresentou o melhor poder preditivo.
Decision Tree
As árvores de decisão são um dos modelos mais utilizados em inferência indutiva. É um instrumento de apoio à tomada de decisão que consiste em uma representação gráfica, a partir de uma decisão inicial, de todas alternativas possíveis.
Para criar e visualizar uma árvore de decisão utilizamos o pacote rpart. Veja um exemplo a seguir. No gráfico observamos a árvore de decisão que agrupa as observações de acordo com suas variáveis. Ela possui 2 componentes principais, folhas e nós.
As folhas representam um grupo de observações. Para cada folha, uma resposta é dada, “Sim” ou “Não”. Abaixo dessas respostas, os números representam a porcentagem de churn em uma folha e, finalmente, vemos a porcentagem de observações totais na folha;
Os nós mostram qual variável foi usada para separar uma folha em duas sub-folhas.
O contrato é a principal variável na decisão de churn. Faz sentido, porque é mais difícil mudar os provedores de telecomunicações se os clientes tiverem um contrato de longo prazo (1 ou 2 anos, neste caso) do que um contrato de mês a mês (deve haver cláusula contratual que não permita que o cliente saia). Assim, clientes com contratos mais curtos estão mais propensos a sair.
if (!require("rpart")) install.packages("rpart")
library("rpart")
if (!require("rpart.plot")) install.packages("rpart.plot")
library("rpart.plot")
tree = rpart(Churn ~., data = db_train %>%
select(-customerID), method = "class")
rpart.plot(tree)
Random Florest
Imagine que você quer comprar o novo álbum do Zeca Pagodinho. Você pergunta a um amigo se você deve comprá-lo. Seu amigo lhe diz: “Compre”. Mas você não tem certeza de que seu amigo e você compartilhem o mesmo gosto musical, então você decide fazer a mesma pergunta para outros 29 amigos. No final, você tem 20 profissionais a favor e 10 contras. Na maioria, você decide comprá-lo. Floresta aleatória funciona da mesma maneira!
As florestas aleatórias ou as florestas de decisão aleatórias são um método em modelagem estatística tanto para classificação quanto para regressão. Como o próprio nome já diz, essa ferramenta utiliza várias árvores de decisão com diferentes opções de váriaveis de entrada e cria um ranking para as categorias “mais votadas”. Assim, o modelo final é determinado a partir da categoria mais significante do ranking criado.
Uma grande vantagem desse método é que ele corrige o overfitting 3 no conjunto de treinamento, algo comum nas árvores de decisão e muito ruim no aspecto de previsão dos dados. Além disso, o método se mostra bastante eficiente para grandes bases de dados.
Aqui, usamos um método de validação cruzada de 5 vezes. Não temos tantas variáveis, então escolhemos 5 valores diferentes para otimizar o número de variáveis via o argumento da função tuneLenght.
if (!require("randomForest")) install.packages("randomForest")
library("randomForest")
ctrl = trainControl(method = "cv", number=5,
classProbs = TRUE, summaryFunction = twoClassSummary)
model.rf = train(Churn ~., data = db_train %>% select(-customerID),
method = "rf",
ntree = 75,
tuneLength = 5,
metric = "ROC",
trControl = ctrl)
model.rf## Random Forest
##
## 4924 samples
## 19 predictor
## 2 classes: 'No', 'Yes'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 3939, 3939, 3939, 3940, 3939
## Resampling results across tuning parameters:
##
## mtry ROC Sens Spec
## 2 0.8282532 0.9219917 0.4682899
## 7 0.8186202 0.9004149 0.4965576
## 12 0.8158146 0.8899032 0.4919803
## 17 0.8123449 0.8946058 0.4889298
## 23 0.8100009 0.8882434 0.5003831
##
## ROC was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.
Logistic Regression
O modelo de regressão logística é semelhante ao modelo de regressão linear. No entanto, no modelo logístico a variável resposta \(Y_i\) é binária. Uma variável binária assume dois valores, como por exemplo, \(Y_i=0\) e \(Y_i=1\), denominados “fracasso” e “sucesso”, respectivamente.
Assim, precisamos recodificar a variável “churn” onde 0 seja “Não” e 1 seja “Sim”. Além disso, temos que escolher as combinações de variáveis que melhor explicarão a decisão de churn.
db_train$ChurnNum = ifelse(db_train$Churn == "Yes",1,0)
good_model = step(glm(ChurnNum ~.,data = db_train %>%
select(-customerID, -Churn ),
family=binomial(link='logit')),
direction="both")## Start: AIC=4115.6
## ChurnNum ~ gender + SeniorCitizen + Partner + Dependents + PhoneService +
## MultipleLines + InternetService + OnlineSecurity + OnlineBackup +
## DeviceProtection + TechSupport + StreamingTV + StreamingMovies +
## Contract + PaperlessBilling + PaymentMethod + MonthlyCharges +
## TotalCharges + tenure
##
## Df Deviance AIC
## - Partner 1 4067.6 4113.6
## - TechSupport 1 4067.7 4113.7
## - Dependents 1 4067.9 4113.9
## - gender 1 4068.0 4114.0
## - OnlineSecurity 1 4068.1 4114.1
## - OnlineBackup 1 4068.2 4114.2
## - PhoneService 1 4068.3 4114.3
## <none> 4067.6 4115.6
## - DeviceProtection 1 4069.8 4115.8
## - MonthlyCharges 1 4070.1 4116.1
## - SeniorCitizen 1 4070.5 4116.5
## - StreamingMovies 1 4071.0 4117.0
## - StreamingTV 1 4071.1 4117.1
## - InternetService 2 4073.7 4117.7
## - MultipleLines 1 4075.5 4121.5
## - PaperlessBilling 1 4076.3 4122.3
## - TotalCharges 1 4081.4 4127.4
## - PaymentMethod 3 4088.7 4130.7
## - Contract 2 4120.3 4164.3
## - tenure 1 4144.6 4190.6
##
## Step: AIC=4113.62
## ChurnNum ~ gender + SeniorCitizen + Dependents + PhoneService +
## MultipleLines + InternetService + OnlineSecurity + OnlineBackup +
## DeviceProtection + TechSupport + StreamingTV + StreamingMovies +
## Contract + PaperlessBilling + PaymentMethod + MonthlyCharges +
## TotalCharges + tenure
##
## Df Deviance AIC
## - TechSupport 1 4067.8 4111.8
## - gender 1 4068.0 4112.0
## - Dependents 1 4068.1 4112.1
## - OnlineSecurity 1 4068.2 4112.2
## - OnlineBackup 1 4068.2 4112.2
## - PhoneService 1 4068.3 4112.3
## <none> 4067.6 4113.6
## - DeviceProtection 1 4069.8 4113.8
## - MonthlyCharges 1 4070.1 4114.1
## - SeniorCitizen 1 4070.5 4114.5
## - StreamingMovies 1 4071.1 4115.1
## - StreamingTV 1 4071.1 4115.1
## + Partner 1 4067.6 4115.6
## - InternetService 2 4073.8 4115.8
## - MultipleLines 1 4075.5 4119.5
## - PaperlessBilling 1 4076.4 4120.4
## - TotalCharges 1 4081.4 4125.4
## - PaymentMethod 3 4088.7 4128.7
## - Contract 2 4120.3 4162.3
## - tenure 1 4145.4 4189.4
##
## Step: AIC=4111.75
## ChurnNum ~ gender + SeniorCitizen + Dependents + PhoneService +
## MultipleLines + InternetService + OnlineSecurity + OnlineBackup +
## DeviceProtection + StreamingTV + StreamingMovies + Contract +
## PaperlessBilling + PaymentMethod + MonthlyCharges + TotalCharges +
## tenure
##
## Df Deviance AIC
## - gender 1 4068.2 4110.2
## - Dependents 1 4068.2 4110.2
## - OnlineSecurity 1 4068.3 4110.3
## <none> 4067.8 4111.8
## - OnlineBackup 1 4070.6 4112.6
## - SeniorCitizen 1 4070.7 4112.7
## + TechSupport 1 4067.6 4113.6
## + Partner 1 4067.7 4113.7
## - PhoneService 1 4072.8 4114.8
## - DeviceProtection 1 4075.6 4117.6
## - PaperlessBilling 1 4076.5 4118.5
## - TotalCharges 1 4081.4 4123.4
## - MonthlyCharges 1 4083.9 4125.9
## - StreamingTV 1 4084.4 4126.4
## - StreamingMovies 1 4084.8 4126.8
## - PaymentMethod 3 4088.9 4126.9
## - MultipleLines 1 4093.8 4135.8
## - InternetService 2 4102.4 4142.4
## - Contract 2 4121.1 4161.1
## - tenure 1 4145.5 4187.5
##
## Step: AIC=4110.16
## ChurnNum ~ SeniorCitizen + Dependents + PhoneService + MultipleLines +
## InternetService + OnlineSecurity + OnlineBackup + DeviceProtection +
## StreamingTV + StreamingMovies + Contract + PaperlessBilling +
## PaymentMethod + MonthlyCharges + TotalCharges + tenure
##
## Df Deviance AIC
## - Dependents 1 4068.6 4108.6
## - OnlineSecurity 1 4068.7 4108.7
## <none> 4068.2 4110.2
## - OnlineBackup 1 4071.0 4111.0
## - SeniorCitizen 1 4071.1 4111.1
## + gender 1 4067.8 4111.8
## + TechSupport 1 4068.0 4112.0
## + Partner 1 4068.1 4112.1
## - PhoneService 1 4073.2 4113.2
## - DeviceProtection 1 4075.9 4115.9
## - PaperlessBilling 1 4077.0 4117.0
## - TotalCharges 1 4081.8 4121.8
## - MonthlyCharges 1 4084.3 4124.3
## - StreamingTV 1 4084.8 4124.8
## - StreamingMovies 1 4085.2 4125.2
## - PaymentMethod 3 4089.3 4125.3
## - MultipleLines 1 4094.2 4134.2
## - InternetService 2 4102.8 4140.8
## - Contract 2 4121.4 4159.4
## - tenure 1 4145.9 4185.9
##
## Step: AIC=4108.63
## ChurnNum ~ SeniorCitizen + PhoneService + MultipleLines + InternetService +
## OnlineSecurity + OnlineBackup + DeviceProtection + StreamingTV +
## StreamingMovies + Contract + PaperlessBilling + PaymentMethod +
## MonthlyCharges + TotalCharges + tenure
##
## Df Deviance AIC
## - OnlineSecurity 1 4069.1 4107.1
## <none> 4068.6 4108.6
## - OnlineBackup 1 4071.5 4109.5
## - SeniorCitizen 1 4072.1 4110.1
## + Dependents 1 4068.2 4110.2
## + gender 1 4068.2 4110.2
## + Partner 1 4068.4 4110.4
## + TechSupport 1 4068.5 4110.5
## - PhoneService 1 4073.7 4111.7
## - DeviceProtection 1 4076.4 4114.4
## - PaperlessBilling 1 4077.4 4115.4
## - TotalCharges 1 4082.5 4120.5
## - MonthlyCharges 1 4084.8 4122.8
## - StreamingTV 1 4085.3 4123.3
## - StreamingMovies 1 4085.8 4123.8
## - PaymentMethod 3 4090.0 4124.0
## - MultipleLines 1 4094.8 4132.8
## - InternetService 2 4103.5 4139.5
## - Contract 2 4123.1 4159.1
## - tenure 1 4147.3 4185.3
##
## Step: AIC=4107.13
## ChurnNum ~ SeniorCitizen + PhoneService + MultipleLines + InternetService +
## OnlineBackup + DeviceProtection + StreamingTV + StreamingMovies +
## Contract + PaperlessBilling + PaymentMethod + MonthlyCharges +
## TotalCharges + tenure
##
## Df Deviance AIC
## <none> 4069.1 4107.1
## - SeniorCitizen 1 4072.5 4108.5
## + OnlineSecurity 1 4068.6 4108.6
## + Dependents 1 4068.7 4108.7
## + gender 1 4068.7 4108.7
## + Partner 1 4068.9 4108.9
## + TechSupport 1 4069.0 4109.0
## - OnlineBackup 1 4074.4 4110.4
## - PaperlessBilling 1 4078.1 4114.1
## - PhoneService 1 4080.8 4116.8
## - DeviceProtection 1 4081.9 4117.9
## - TotalCharges 1 4083.0 4119.0
## - PaymentMethod 3 4090.4 4122.4
## - StreamingTV 1 4101.9 4137.9
## - StreamingMovies 1 4102.5 4138.5
## - MultipleLines 1 4106.8 4142.8
## - MonthlyCharges 1 4106.9 4142.9
## - Contract 2 4123.4 4157.4
## - InternetService 2 4146.5 4180.5
## - tenure 1 4148.1 4184.1
Parece que 13 variáveis são suficientes para o nosso modelo. Agora vejamos as estimativas dos coeficientes.
summary(good_model)##
## Call:
## glm(formula = ChurnNum ~ SeniorCitizen + PhoneService + MultipleLines +
## InternetService + OnlineBackup + DeviceProtection + StreamingTV +
## StreamingMovies + Contract + PaperlessBilling + PaymentMethod +
## MonthlyCharges + TotalCharges + tenure, family = binomial(link = "logit"),
## data = db_train %>% select(-customerID, -Churn))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9012 -0.6731 -0.2923 0.7290 3.3888
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.60150 0.22743 2.645 0.008174
## SeniorCitizen1 0.18366 0.09965 1.843 0.065337
## PhoneServiceYes 1.04493 0.30762 3.397 0.000682
## MultipleLinesYes 0.69549 0.11417 6.092 1.12e-09
## InternetServiceFiber optic 2.82341 0.32796 8.609 < 2e-16
## InternetServiceNo -2.87861 0.40080 -7.182 6.86e-13
## OnlineBackupYes 0.25915 0.11367 2.280 0.022619
## DeviceProtectionYes 0.41179 0.11563 3.561 0.000369
## StreamingTVYes 0.92831 0.16380 5.667 1.45e-08
## StreamingMoviesYes 0.92595 0.16194 5.718 1.08e-08
## ContractOne year -0.66307 0.12800 -5.180 2.22e-07
## ContractTwo year -1.24824 0.20212 -6.176 6.59e-10
## PaperlessBillingYes 0.26399 0.08847 2.984 0.002847
## PaymentMethodCredit card (automatic) -0.18326 0.13747 -1.333 0.182523
## PaymentMethodElectronic check 0.30024 0.11163 2.690 0.007152
## PaymentMethodMailed check -0.02901 0.13531 -0.214 0.830262
## MonthlyCharges -8.02745 1.32045 -6.079 1.21e-09
## TotalCharges 2.63684 0.72491 3.637 0.000275
## tenure -4.24828 0.52098 -8.154 3.51e-16
##
## (Intercept) **
## SeniorCitizen1 .
## PhoneServiceYes ***
## MultipleLinesYes ***
## InternetServiceFiber optic ***
## InternetServiceNo ***
## OnlineBackupYes *
## DeviceProtectionYes ***
## StreamingTVYes ***
## StreamingMoviesYes ***
## ContractOne year ***
## ContractTwo year ***
## PaperlessBillingYes **
## PaymentMethodCredit card (automatic)
## PaymentMethodElectronic check **
## PaymentMethodMailed check
## MonthlyCharges ***
## TotalCharges ***
## tenure ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5702.8 on 4923 degrees of freedom
## Residual deviance: 4069.1 on 4905 degrees of freedom
## AIC: 4107.1
##
## Number of Fisher Scoring iterations: 6
Considerando um nível de significância de 5% e a partir do valor p, verificou-se que a maioira das variáveis é significante. Apenas as variáveis SeniorCitizen = 1, PaymentMethod = 'Credit card (automatic)' e PaymentMethod = 'Mailed check' não são significantes.
Como já era esperado, as estimativas negativas das variáveis de contrato (contract) de um ou dois anos, confirmou que contratos mais longos diminuem as chances de churn do cliente. Foi verificado também, que quanto maior for as taxas mensais, maior é a probabilidade do cliente sair, uma vez que a estimativa do coeficiente da variável MonthlyCharges foi de -8.02745.
Seleção do modelo
Agora que temos nossos 3 modelos, decidimos testá-los em nosso conjunto de dados de teste. Antes disso, precisamos aplicar as mesmas transformações que aplicamos anteriormente ao conjunto de treinamento.
norm.test = lapply(db_test %>% select(MonthlyCharges,TotalCharges, tenure), normalize)
norm.test = do.call(cbind, norm.test) %>%
as.data.frame()
factor.test = lapply(db_test %>%
select(-customerID,-MonthlyCharges,
-TotalCharges, -tenure),
function(x){
x = gsub("No internet service", "No", x)
x = gsub("No phone service", "No", x)
return(x)
})
factor.test = do.call(cbind, factor.test) %>% as.data.frame()
db_test = cbind( customerID = db_test[,1], factor.test , norm.test)
db_test$ChurnNum = ifelse(db_test$Churn == "Yes", 1, 0)
Agora, realizamos a previsão nos dados de teste utilizando a função predict() para cada tipo de modelo e obtemos suas estimativas.
# Previsão
# Árvore de decição
pred_tree = predict(tree, db_test %>%
select(-customerID), type = c("prob"))[,2]
# Floresta aleatória
pred_rf = predict(object=model.rf, db_test %>%
select(-customerID), type='prob')[,2]
# Regressão logística
pred_logistic = predict(good_model, db_test, type="response")
Com as previsões calculadas, utilizaremos a seguir a curva de ROC para verificar qual modelo teve o melhor ajuste. Sabe-se que quão maior a área sob a curva, melhor foi a precisão do modelo.
if (!require("ROSE")) install.packages("ROSE")
library("ROSE")
roc_tree = roc.curve(response = db_test$ChurnNum, pred_tree,
col = "#0d84da", main = "Curva ROC",
ylab = "Taxa de verdadeiro positivo", xlab = "Taxa de falso positivo")
roc_rf = roc.curve(response = db_test$ChurnNum, pred_rf,
col = "#ef0a30", add.roc=TRUE)
roc_logistic = roc.curve(response = db_test$ChurnNum, pred_logistic,
col = "#45a163", add.roc=TRUE)
legend("bottomright", legend=c("Decision tree", "Randon Forest",
"Logistic regression"),
lwd = 2, col = c("#0d84da", "#ef0a30", "#45a163"))
Note que as curvas dos modelos de floresta aleatória e regressão logística estão bem próximas, porém, o logístico apresenta uma leve vantagem ao decorrer da curva. Dessa forma, e por se tratar de um modelo mais simples, escolhe-se o modelo de regressão logística.
Parte 2 - Indicadores
Na seção anterior foi escolhido o modelo de regressão logística como o mais apropriado para os nossos dados. Obteve-se a seguinte curva ROC com um área sob a curva (AUC) de 0,843, classificada na literatura como boa.
if (!require("ROSE")) install.packages("ROSE")
library("ROSE")
roc_logistic = roc.curve(response = db_test$ChurnNum, pred_logistic,
col = "#45a163")
Agora imaginamos dois cenários:
1. Supomos que temos um orçamento grande e queremos segmentar muitos clientes. Vamos nos
concentrar no valor limite. Queremos otimizar o número de boas classificações.
2. Temos um orçamento restrito. Por isso, queremos entrar em contato com apenas 25% dos
clientes em nosso banco de dados. Concentraremos-nos na probabilidade de churn dos clientes,
também conhecida como pontuação. Em seguida, selecionaremos os 25% melhores clientes com
base na sua probabilidade de churn usando uma curva de elevação.Cenário 1: Optimização de limites (Threshold Optimization)
Cada cliente tem uma pontuação que corresponde à sua probabilidade de churn. Vamos ver a pontuação dos cinco primeiros clientes do nosso banco de dados.
head(pred_logistic,5)## 1 2 3 4 5
## 0.62513394 0.32512336 0.52812027 0.31273422 0.01264661Originalmente, o limite é 0.5. Predizemos uma resposta positiva se a pontuação for maior do que 0.5 e uma resposta negativa se a pontuação for menor que 0.5. Significa que para os 5 clientes acima que o modelo prevê uma resposta positiva (churn) para o cliente 1 e para o cliente 3 e uma resposta negativa para os outros 3 clientes (não churn).
Para otimizar o limite, queremos calcular diferentes indicadores estatísticos (acurácia, precisão, sensibilidade, f1 e kappa4) para diferentes valores de limites de 0.05 para 0.847 com um valor de passo de 0.001. Não usaremos valores acima 0.847 pois não existem scores maiores que esse valor na base de teste (veja abaixo).
quantile(pred_logistic)## 0% 25% 50% 75% 100%
## 0.001383889 0.048880103 0.200493493 0.499503150 0.847223943
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse")
if (!require("caret")) install.packages("caret")
library("caret")
if (!require("e1071")) install.packages("e1071")
library("e1071")
comp = cbind.data.frame(answer = db_test$ChurnNum,
pred=pred_logistic) %>%
arrange(desc(pred))
indic_perf = function(x){
compare = comp %>%
mutate(pred = ifelse(pred>x,1,0))
if(ncol(table(compare))>1){
mat = confusionMatrix(table(compare), positive = "1")
#acuracy
acc = mat$overall["Accuracy"]
#Kappa
kap = mat$overall["Kappa"]
#sensitivity
sen = mat$byClass["Sensitivity"]
#F1
f1_stat = mat$byClass["F1"]
#Precision
prec = mat$byClass["Precision"]
}else{
acc = NA
prec = NA
sen = NA
kap = NA
f1_stat = NA
}
return(data.frame(threshold = x, accuracy = acc,
precision = prec, sensitivity = sen,
kappa = kap, f1= f1_stat))
}
indics = do.call(rbind, lapply(seq(0.05,0.95, by=0.001),
indic_perf)) %>%
filter(!is.na(accuracy))
if (!require("plotly")) install.packages("plotly")
library("plotly")
if (!require("IRdisplay")) install.packages("IRdisplay")
library("IRdisplay")
gather_indics = tidyr::gather(indics, variable,
value, -threshold) %>%
group_by(variable) %>%
mutate(color = (max(value) == value),
threshold = as.numeric(threshold) )
ggplot(gather_indics , aes(x= threshold, y=value)) +
ggtitle("Indicator values by thresholds")+
geom_point(aes(color = color), size=0.5) +
facet_wrap(~variable, scales = 'free_x') +
scale_color_manual(values = c(NA, "tomato")) +
labs(x="thresholds", y=" ") +
geom_line(color="navy") + theme_bw()+
theme( legend.position="none")
Nós desenhamos uma tabela onde mostramos o valor máximo para cada indicador [threshold entre 0.050 e 0.057 == max precisão; threshold igual a 0.339 == max f1; threshold entre 0.584 e 0.585 == max acurácia e assim sucessivamente).
max_indics = indics %>%
filter(accuracy == max(accuracy, na.rm=TRUE) |
precision == max(precision, na.rm = TRUE) | sensitivity == max(sensitivity, na.rm = TRUE) |
kappa == max(kappa, na.rm = TRUE) | f1 == max(f1, na.rm = TRUE) )
max_indics
Ainda temos muitos valores, então decidimos remover algumas linhas. Mantivemos os cinco limites mais relevantes: 0.057, 0.339, 0.584, 0.585 e 0.799.
max_indics %>%
filter( threshold %in% c("0.057", "0.339", "0.584", "0.585", "0.799"))
Nesta fase, podemos escolher um limite com base no que consideramos ser o mais importante. Por questões de orçamento vamos continuar com 0.339 . Por isso, temos a seguinte matriz de confusão.
compare = comp %>%
mutate(pred = ifelse(pred>0.339,1,0))
confusionMatrix(table(compare), positive = "1")## Confusion Matrix and Statistics
##
## pred
## answer 0 1
## 0 1206 342
## 1 139 421
##
## Accuracy : 0.7718
## 95% CI : (0.7533, 0.7896)
## No Information Rate : 0.638
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4758
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.5518
## Specificity : 0.8967
## Pos Pred Value : 0.7518
## Neg Pred Value : 0.7791
## Prevalence : 0.3620
## Detection Rate : 0.1997
## Detection Prevalence : 0.2657
## Balanced Accuracy : 0.7242
##
## 'Positive' Class : 1
##
Com um limite de 0.339, temos o melhor valor f1 e temos um nível aceitável de acurácia, precisão e sensibilidade. Agora vamos nos concentrar no segundo cenário.
Cenário 2: curva de elevação (Lift curve)
Agora queremos atingir um tipo de cliente (como maior probabilidade de sair) em nosso banco de dados. Vamos classificar os clientes por seus scores e então, observar a taxa de churn em diferentes percentuais. A melhor maneira de visualizar a taxa de churn por percentil é fazendo uma curva de elevação, isto é, classificar os clientes por sua probabilidade de sair e depois separá-los em percentis.
Quanto maior a probabilidade de sair, maior será a classificação. Isso significa que o 5º percentil representa os 5% de clientes que têm a maior probabilidade de sair de acordo com nosso modelo.
A propósito, sabemos que há 26,54% de churn em nosso banco de dados. Portanto, se tivermos um modelo aleatório, espera-se que haja 26,54% de churn em cada percentil. Assim, podemos comparar a taxa de churn no nosso modelo e no aleatório.
fivtile = nrow(comp)/20 # total de clientes por percentil 5% = 105
step = floor(fivtile * 1:20) # definiu os percentis na base de dados
pct = sapply(step, function(x){
return(mean(comp[1:x,1]))}) # calcula a média de churn por percentil.
#paste(seq(from = 5, to = 100, by=5), "%",sep=" ")
lift = data.frame(label= seq(from = 5, to = 100, by=5), score = pct*100)
q = ggplot(lift, aes(x=label, y=score))+
geom_bar(stat="identity",position="stack",color="navy", fill="navy")+
ggtitle("Taxa de churn pelo percentual acumulado de clientes
\n com alta probabilidade de sair")+
coord_cartesian(xlim = c(5,100), ylim = c(0,100))+
scale_y_continuous(breaks = seq(from = 0, to = 100, by=25),
labels = function(x) paste0(x,"%", sep = ""))+
scale_x_continuous(breaks = c(5, 25, 50, 75, 100),
labels = function(x) paste0(x,"%", sep = ""))+
labs(x="percentil acumulado", y="taxa de churn") +
geom_line(aes(y=score[20]),linetype=2, size=1, color="tomato")+
theme_minimal()+
theme(
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
plot.title = element_text(size=17,face="bold", hjust=0.5),
axis.text.x = element_text(vjust=1),
axis.ticks.x = element_blank(),
axis.title.x = element_text(size=13, face="bold", vjust=-0.4),
axis.title.y = element_text(size=13,face="bold")#,
#strip.text.x = element_text(face="italic", size=11)
)
print(q)
Nosso modelo é realmente eficiente comparado ao modelo aleatório. Detectamos pelo menos 2 vezes mais churn usando nosso modelo até o 40º percentil. Se segmentarmos 25% de nosso banco de dados, observamos que a taxa de cancelamentos é de 61,67%. Tal iniciativa permite que uma empresa com um orçamento baixo para combater a taxa de churn possa ter maior foco em um determinado tipo de cliente, mais que isso, com a curva de elevação os C-levels podem decidir o que desejam atacar.
Conclusão
Na primeira parte limpamos o banco de dados, fizemos uma análise exploratória de dados e testamos diferentes modelos para finalmente escolher aquele com a maior área sob a curva (AUC), a regressão logística.
Na segunda parte, imaginamos dois cenários diferentes. No primeiro, nos concentramos no valor do limite (threshold value) e, no segundo, queremos segmentar apenas 25% do banco de dados. Assim, na primeira parte, encontramos o valor limite otimizado: 0.339. Na segunda parte, vimos que nosso modelo é eficiente para 25% dos clientes da nossa base de dados. A taxa esperada de churn sem o modelo é de 26,54%, enquanto a taxa de churn é de 61,67% para o nosso modelo.
Neste estudo será feita uma análise de churn em uma base de dados fornecida pela IBM. Para obtenção da base e mais informações clique aqui. Para mais informações sobre o estudo original clique aqui e aqui↩
tenure = 1/churn. Para maiores detalhes, acesse: https://rowansimpson.com/2014/11/26/estimating-tenure/↩
overfitting é um termo usado quando um modelo se ajusta muito bem a um conjunto de dados previamente observado, mas se mostra ruim e ineficaz para previsão de novos dados.↩
Esse artigo mostra como usar o pacote Caret e explica a estatística Kappa - http://www.jstatsoft.org/article/view/v028i05/v28i05.pdf↩