Synopsis

In order to helo Regork with their telecommunication efforts, we have prepared a report and a detailed analysis. We explored patterns within customer status to determine what could cause customers to leave and what could cause them to stay. We also asses three different algorithms to compare model performance. We came up with some logical reasons why we feel some variables are predictors.

Library

These following packages will help reproduce the results contained within this report

library(tidymodels)
## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
## ✔ broom        1.0.6     ✔ recipes      1.1.0
## ✔ dials        1.3.0     ✔ rsample      1.2.1
## ✔ dplyr        1.1.4     ✔ tibble       3.2.1
## ✔ ggplot2      3.5.1     ✔ tidyr        1.3.1
## ✔ infer        1.0.7     ✔ tune         1.2.1
## ✔ modeldata    1.4.0     ✔ workflows    1.1.4
## ✔ parsnip      1.2.1     ✔ workflowsets 1.1.0
## ✔ purrr        1.0.2     ✔ yardstick    1.3.1
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ recipes::step()  masks stats::step()
## • Search for functions across packages at https://www.tidymodels.org/find/
library(tidyverse) 
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ lubridate 1.9.3     ✔ stringr   1.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ purrr::discard()    masks scales::discard()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ stringr::fixed()    masks recipes::fixed()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ readr::spec()       masks yardstick::spec()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(baguette)
library(vip) 
## 
## Attaching package: 'vip'
## 
## The following object is masked from 'package:utils':
## 
##     vi
library(pdp) 
## 
## Attaching package: 'pdp'
## 
## The following object is masked from 'package:purrr':
## 
##     partial
library(here)
## here() starts at /Users/trungphan
library(kernlab)
## 
## Attaching package: 'kernlab'
## 
## The following object is masked from 'package:purrr':
## 
##     cross
## 
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## 
## The following object is masked from 'package:scales':
## 
##     alpha
library(rpart.plot)
## Loading required package: rpart
## 
## Attaching package: 'rpart'
## 
## The following object is masked from 'package:dials':
## 
##     prune
library(ranger)
library(ggplot2)
library(earth)
## Loading required package: Formula
## Loading required package: plotmo
## Loading required package: plotrix
## 
## Attaching package: 'plotrix'
## 
## The following object is masked from 'package:scales':
## 
##     rescale
library(readr)

Data Preparation

The following data will be stored within the environment for code throughout.

retention <- read_csv("Downloads/customer_retention.csv")
## Rows: 6999 Columns: 20
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (16): Gender, Partner, Dependents, PhoneService, MultipleLines, Internet...
## dbl  (4): SeniorCitizen, Tenure, MonthlyCharges, TotalCharges
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
retention <- retention %>%
  dplyr::mutate(Status = as.factor(Status))

retention <- drop_na(retention)
retention %>% 
  is.na() %>%
  sum()
## [1] 0

Exploratory Analysis

Percentage of churn.

left <- sum(retention$Status == 'Left')

churn <- left/6988
# Contract Type and Customer Status
contract_type <- retention %>%
  group_by(Contract, Status) %>%
  summarise(count = n()) %>%
  ungroup()
## `summarise()` has grouped output by 'Contract'. You can override using the
## `.groups` argument.
ggplot(contract_type, aes(x = Contract, y = count, fill = Status)) +
  geom_bar(position = "dodge", stat = "identity") +
  scale_y_continuous(name = "Number of Customers", labels = comma_format()) +
  xlab("Contract Type") +
  ggtitle("Contract Type vs Customer Status") +
  scale_fill_manual(values = c("darkblue", "lightblue")) +
  theme(legend.position = "bottom")

It appears that month-to-month is the most common contract type so therefore, it contains the most amount of customers that are current and have left. One thing to note, is that the one year contract has the lowest number of current customers and a higher number of customers that have left than the two year contract.

library(ggplot2)
library(dplyr)
library(scales)

retention_summary <- retention %>%
  filter(InternetService != "No") %>%
  group_by(OnlineSecurity, Status) %>%
  summarise(count = n(), .groups = 'drop')

ggplot(retention_summary, aes(x = OnlineSecurity, y = count, fill = Status)) +
  geom_bar(stat = "identity", position = "dodge") +  
  scale_y_continuous(name = "Number of Customers", labels = comma_format()) +
  xlab("Online Security") +
  ggtitle("Online Security vs Customer Status") +
  scale_fill_manual(values = c("darkblue", "lightblue")) +  
  theme(legend.position = "bottom")

More customers do not have online security. However, it does appear that the number of customers that have left drops significantly when they do have online security. This may be correlated because the number of current customers did not drop that much when they you go from no online security to having online security.

graph3 <- retention %>%
  filter(PaymentMethod != "No") %>%
  group_by(PaymentMethod, Status) %>%
  summarise(count = n(), .groups = 'drop')

ggplot(graph3, aes(x = PaymentMethod, y = count, fill = Status)) +
  geom_bar(stat = "identity", position = "dodge") +  
  ggtitle("Payment Method vs Customer Status") +
  theme(axis.text.x = element_text(angle = 50, size = 7, vjust = 0.5)) +
  labs(y = "Count of Customers", x = "Payment Method") +
  scale_fill_manual(values = c("darkblue", "lightblue")) +  
  theme(legend.position = "bottom")

From just looking at this graph it seems like electronic check customers leave the most. The other payment methods are pretty even in terms of number of current customers and those that have left.

graph4 <- retention %>%
  filter(InternetService != "No") %>%
  group_by(InternetService, Status) %>%
  summarise(count = n(), .groups = 'drop')

ggplot(graph4, aes(x = InternetService, y = count, fill = Status)) +
  geom_bar(stat = "identity", position = "dodge") +  
  ggtitle("Internet Service vs Customer Status") +
  labs(y = "Count of Customers", x = "Internet Service") +
  scale_fill_manual(values = c("darkblue", "lightblue")) +  
  theme(legend.position = "bottom")

From this graph we can see that significantly more customers left when they had Fiber obtic internet service.

Machine Learning

Logistic Regression

For the first machine learning model, we are doing logistic regression model

set.seed(123)
split <- initial_split(retention, prop = 0.7, strata = "Status")
retention_train <- training(split)
retention_test <- testing(split)
log_recipe <- recipe(Status ~ ., data = retention_train) %>%
  step_normalize(all_numeric_predictors()) %>%
  step_dummy(all_nominal_predictors())
logistic_spec <- logistic_reg() %>% 
  set_engine("glm")

workflow <- workflow() %>%
  add_model(logistic_spec) %>%
  add_recipe(log_recipe)
set.seed(123)
cv_folds <- vfold_cv(retention_train, v = 5, strata = "Status")
set.seed(123)
logistic_results <- workflow %>%
  fit_resamples(
    resamples = cv_folds,
    metrics = metric_set(roc_auc, accuracy),
    control = control_resamples(save_pred = TRUE)
  )
collect_metrics(logistic_results)
## # A tibble: 2 × 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy binary     0.799     5 0.00401 Preprocessor1_Model1
## 2 roc_auc  binary     0.845     5 0.00521 Preprocessor1_Model1

We can see that the Logistic Regression model have good results for accuracy and mean AUC. However, we are doing 3 more to have a full comparison.

Multivariate Adaptive Regression Splines (MARS)

For the second model, we ran a MARS (Multivariate Adaptive Regression Splines) algorithm to train a classification model on the customer retention data set.

mars_recipe <- recipe(Status ~ ., data = retention_train)
set.seed(123)
mars_kfolds <- vfold_cv(retention_train, v = 5, strata = "Status")
mars_mod <- mars(num_terms = tune(), prod_degree = tune()) %>%
  set_mode("classification")
mars_grid <- grid_regular(num_terms(range = c(1,30)), prod_degree(), levels = 50)
mars_wf <- workflow() %>% add_recipe(mars_recipe) %>% add_model(mars_mod)
mars_results <- mars_wf %>% tune_grid(resamples = mars_kfolds, grid = mars_grid)
mars_results %>% collect_metrics() %>% filter(.metric == "roc_auc") %>%
arrange(desc(mean))
## # A tibble: 60 × 8
##    num_terms prod_degree .metric .estimator  mean     n std_err .config         
##        <int>       <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>           
##  1        20           1 roc_auc binary     0.850     5 0.00486 Preprocessor1_M…
##  2        19           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
##  3        21           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
##  4        22           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
##  5        23           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
##  6        24           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
##  7        25           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
##  8        26           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
##  9        27           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
## 10        28           1 roc_auc binary     0.849     5 0.00486 Preprocessor1_M…
## # ℹ 50 more rows
autoplot(mars_results)

This MARS model provides great accuracy and AUC mean values. We can find the best hyperparameter values, and the most important features in this model

mars_best <- select_best(mars_results, metric = "roc_auc")

mars_final_wf <- workflow() %>% 
  add_model(mars_mod) %>% add_formula(Status ~ .) %>% 
  finalize_workflow(mars_best)

mars_final_wf %>% 
  fit(data = retention_train) %>%
  extract_fit_parsnip() %>%
  vip(10, type = "rss")

This plot shows the most important feature in this model is Tenure. This means that the longer a customer stays as current customer, they are less likely to leave.

Decision Tree

For the next model, we have Decision Tree

dt_mod <- decision_tree(mode = 'classification') %>%
  set_engine("rpart")
model_recipe <- recipe(Status ~ ., data = retention_train) %>%
  step_normalize(all_numeric_predictors()) %>%
  step_dummy(all_nominal_predictors())
dt_fit <- workflow() %>%
  add_recipe(model_recipe) %>%
  add_model(dt_mod) %>%
  fit(data = retention_train)
rpart.plot::rpart.plot(dt_fit$fit$fit$fit)
## Warning: Cannot retrieve the data used to build the model (so cannot determine roundint and is.binary for the variables).
## To silence this warning:
##     Call rpart.plot with roundint=FALSE,
##     or rebuild the rpart model with model=TRUE.

set.seed(123)
kfold <- vfold_cv(retention_train, v = 5)
dt_results <- fit_resamples(dt_mod, model_recipe, kfold)
collect_metrics(dt_results)
## # A tibble: 3 × 6
##   .metric     .estimator  mean     n std_err .config             
##   <chr>       <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy    binary     0.787     5 0.00560 Preprocessor1_Model1
## 2 brier_class binary     0.163     5 0.00317 Preprocessor1_Model1
## 3 roc_auc     binary     0.703     5 0.00571 Preprocessor1_Model1

The result of Decision Tree shows that this is least accurate model for making our decisions.

TUNING

dt_mod <- decision_tree(
  mode = "classification",
  cost_complexity = tune(),
  tree_depth = tune(),
  min_n = tune()
 ) %>%
  set_engine("rpart")
dt_hyper_grid <- grid_regular(
  cost_complexity(),
  tree_depth(),
  min_n(),
  levels = 5
 )
set.seed(123)
dt_results <- tune_grid(dt_mod, model_recipe, resamples = kfold, grid = dt_hyper_grid)
show_best(dt_results, metric = "roc_auc", n = 5)
## # A tibble: 5 × 9
##   cost_complexity tree_depth min_n .metric .estimator  mean     n std_err
##             <dbl>      <int> <int> <chr>   <chr>      <dbl> <int>   <dbl>
## 1    0.0000000001          8    30 roc_auc binary     0.821     5 0.00894
## 2    0.0000000178          8    30 roc_auc binary     0.821     5 0.00894
## 3    0.00000316            8    30 roc_auc binary     0.821     5 0.00894
## 4    0.0000000001          8    40 roc_auc binary     0.819     5 0.00929
## 5    0.0000000178          8    40 roc_auc binary     0.819     5 0.00929
## # ℹ 1 more variable: .config <chr>
dt_best_model <- select_best(dt_results, metric = 'roc_auc')

dt_final_wf <- workflow() %>%
  add_recipe(model_recipe) %>%
  add_model(dt_mod) %>%
  finalize_workflow(dt_best_model)

dt_final_fit <- dt_final_wf %>%
  fit(data = retention_train)

dt_final_fit %>%
  extract_fit_parsnip() %>%
  vip(20)

Even though this model shows that Total Charges is the most important feature, Tenure is still the second highest on this.

Bagging

Finally, we do Bagging for the last model

bag_mod <- bag_tree() %>%
  set_engine("rpart", times = 5) %>%
  set_mode("classification")
bag_results <- fit_resamples(bag_mod, model_recipe, kfold)
collect_metrics(bag_results)
## # A tibble: 3 × 6
##   .metric     .estimator  mean     n std_err .config             
##   <chr>       <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy    binary     0.769     5 0.00687 Preprocessor1_Model1
## 2 brier_class binary     0.167     5 0.00525 Preprocessor1_Model1
## 3 roc_auc     binary     0.778     5 0.0102  Preprocessor1_Model1

The result for this model are close to Decision Tree, meaning that these two are not the most accurate for this.

TUNING

bag_mod <- bag_tree() %>%
  set_engine("rpart", times = tune()) %>%
  set_mode("classification")
bag_hyper_grid <- expand.grid(times = c(5, 25, 50, 100, 200, 300))
set.seed(123)
bag_results <- tune_grid(bag_mod, model_recipe, resamples = kfold, grid = bag_hyper_grid)
show_best(bag_results, metric = "roc_auc")
## # A tibble: 5 × 7
##   times .metric .estimator  mean     n std_err .config             
##   <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1   200 roc_auc binary     0.819     5 0.0106  Preprocessor1_Model5
## 2   300 roc_auc binary     0.819     5 0.0109  Preprocessor1_Model6
## 3   100 roc_auc binary     0.815     5 0.0105  Preprocessor1_Model4
## 4    50 roc_auc binary     0.811     5 0.0110  Preprocessor1_Model3
## 5    25 roc_auc binary     0.807     5 0.00893 Preprocessor1_Model2

Confusion Matrix

Now we construct a Confusion Matrix to find the validity of the model

confus_matrix <- logistic_reg() %>%
fit(Status ~ ., data = retention_train)

confus_matrix %>% predict(retention_test) %>% 
  bind_cols(retention_test %>% select(Status)) %>%
  conf_mat(truth = Status, estimate = .pred_class)
##           Truth
## Prediction Current Left
##    Current    1362  225
##    Left        178  332

Conclusion

  • Our goal of this project is to identify customer patterns to improve retention of Regork.

  • Numbers of customers that are predicted to leave

# Calculate customers that are predicted to leave
current <- sum(retention$Status == 'Current') 

current*churn
## [1] 1363.05
  • Predicted loss in revenue if no action is taken
# Calculate loss in revenue
revenue_current <- sum(retention$TotalCharges[which(retention$Status == 'Current')])
revenue_current
## [1] 13106622
loss <- revenue_current*churn
loss
## [1] 3481095

  • We found that Tenure is the highest variable that influence the customer status with the company. Alongside, Total Charges and Monthly Charges also play an important role for the customer to whether stay or not.

  • Hence, here are some reccomendations we have after this research process:

    • Offer discounts and loyalty points for specific payment methods or recurring subscriptions.
    • Incentivize one- or two-year contracts with discounts and rewards to make them more attractive than month-to-month options.
    • Promote leveled incentives for different levels of tenure