Introduction: At Regork Telecom, we are currently tackling the crucial business challenge of predicting why customers leave and implementing strategies for retention. This imperative endeavor not only aims to bolster revenue streams but also underscores our commitment to prioritizing customer satisfaction and loyalty. By discerning which customers are likely to leave in the future, we can proactively engage with them, tailor solutions to address their needs, and ultimately fortify our relationship with them. To do this, we will do exploratory data analysis to find trends in the data, perform machine learning through logistic regression, MARS, decision trees, and random forest, in order to find the best model at predicting if someone will leave. Through this, we will be better able to make adjustments to the business in order to retain customers. This initiative exemplifies Regork Telecom’s dedication to delivering exceptional service and ensuring that our valued customers remain at the forefront of our business priorities.

Libraries Required

library(tidymodels)
library(tidyverse)
library(baguette)
library(vip)
library(pdp)
library(here)
library(kernlab)
library(ggplot2)
library(earth)

Importing the data

Import and clean the data:

retention <- read.csv('Labs/data/customer_retention.csv')
sum(is.na(retention))
## [1] 11

Cleaning up dataset:

retention <- retention %>%
  mutate(Status = factor(Status)) %>%
  na.omit 
sum(is.na(retention))
## [1] 0

Exploratory Analysis:

os_counts <- retention %>%
  count(OnlineSecurity) %>%
  arrange(desc(n))

ggplot(os_counts, aes(x = fct_reorder(OnlineSecurity, n), y = n)) +
  geom_col(fill = "skyblue", width = 0.7) +
  theme_minimal() +
  labs(title = "Presence of Online Security Among Users Who Left",
       x = "Online Security",
       y = "Number of Users") +
  coord_flip() +  
  theme(axis.text.y = element_text(size = 10))

contract_counts <- retention %>%
  count(Contract) %>%
  arrange(desc(n))

# Plotting
ggplot(contract_counts, aes(x = fct_reorder(Contract, n), y = n)) +
  geom_col(fill = "skyblue", width = 0.7) +
  theme_minimal() +
  labs(title = "Contract Types of Users Who Left",
       x = "Contract Type",
       y = "Number of Users") +
  coord_flip() +  # Flipping axes for better readability
  theme(axis.text.y = element_text(size = 10))

Here we see there is a trend in customers leaving when they have month to month contracts.

tenure_counts <- retention %>%
  count(Tenure) %>%
  arrange(desc(n))


ggplot(tenure_counts, aes(x = Tenure, y = n)) +
  geom_line(color = "skyblue", linewidth = 1) +
  geom_point(color = "skyblue", linewidth = 3) +  
  theme_minimal() +
  labs(title = "Distribution of Tenure Among Departed Users",
       x = "Months",
       y = "Number of Users")

As tenure increases, the customer leaving decreases until 60 months, then it increases again.

ggplot(retention, aes(Tenure)) + 
  geom_bar(fill = "pink") +
  ggtitle("Partner and Tenure") +
  facet_wrap(~Partner)

This graph shows that people who are married are less likely to leave.

Next, lets move onto machine learning.

Machine Learning:

For our first model, we will look at logistic regression.

set.seed(123)
log_split <- initial_split(retention, prop = .7, strata = Status)
log_train <- training(log_split)
log_test <- testing(log_split)
set.seed(123)
log_kfolds <- vfold_cv(log_train, v = 5, strata = Status)
results <- logistic_reg() %>%
  fit_resamples(Status ~ ., log_kfolds)
collect_metrics(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 found an roc_auc value of .8448, which is good, but should be able to be improved.

final_fit <- logistic_reg() %>%
  fit(Status ~ ., data = log_train)
tidy(final_fit)
## # A tibble: 31 × 5
##    term                          estimate std.error statistic   p.value
##    <chr>                            <dbl>     <dbl>     <dbl>     <dbl>
##  1 (Intercept)                     1.02     0.969      1.05    2.95e- 1
##  2 GenderMale                     -0.0461   0.0779    -0.592   5.54e- 1
##  3 SeniorCitizen                   0.258    0.101      2.55    1.08e- 2
##  4 PartnerYes                     -0.141    0.0925    -1.52    1.29e- 1
##  5 DependentsYes                  -0.0475   0.108     -0.441   6.59e- 1
##  6 Tenure                         -0.0646   0.00762   -8.47    2.35e-17
##  7 PhoneServiceYes                -0.0417   0.774     -0.0539  9.57e- 1
##  8 MultipleLinesNo phone service  NA       NA         NA      NA       
##  9 MultipleLinesYes                0.442    0.211      2.10    3.56e- 2
## 10 InternetServiceFiber optic      1.48     0.950      1.56    1.18e- 1
## # ℹ 21 more rows
final_fit %>%
  predict(log_test) %>%
  bind_cols(log_test %>% select(Status)) %>%
  conf_mat(truth = Status, estimate = .pred_class)
##           Truth
## Prediction Current Left
##    Current    1362  225
##    Left        178  332

We find that our model is usually predicting correctly, but we are experiencing some false positives.

Next up, we will assess a Mars Model

set.seed(123) # for reproducibility
mars_split <- initial_split(retention, prop = .7, strata = Status)
mars_train <- training(mars_split)
mars_test <- testing(mars_split)
mars_recipe <- recipe(Status ~ ., data = mars_train) %>%
  step_YeoJohnson(all_numeric_predictors()) %>%
  step_normalize(all_numeric_predictors())
set.seed(123)
kfolds <- vfold_cv(mars_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 = 40)
mars_wf <- workflow() %>%
  add_recipe(mars_recipe) %>%
  add_model(mars_mod)
tuning_results <- mars_wf %>%
  tune_grid(resamples = kfolds, grid = mars_grid)
tuning_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        15           1 roc_auc binary     0.850     5 0.00509 Preprocessor1_M…
##  2        16           1 roc_auc binary     0.849     5 0.00502 Preprocessor1_M…
##  3        18           1 roc_auc binary     0.849     5 0.00503 Preprocessor1_M…
##  4        19           1 roc_auc binary     0.849     5 0.00503 Preprocessor1_M…
##  5        20           1 roc_auc binary     0.849     5 0.00503 Preprocessor1_M…
##  6        21           1 roc_auc binary     0.849     5 0.00503 Preprocessor1_M…
##  7        22           1 roc_auc binary     0.849     5 0.00503 Preprocessor1_M…
##  8        23           1 roc_auc binary     0.849     5 0.00503 Preprocessor1_M…
##  9        24           1 roc_auc binary     0.849     5 0.00503 Preprocessor1_M…
## 10        25           1 roc_auc binary     0.849     5 0.00503 Preprocessor1_M…
## # ℹ 50 more rows

Here we find a max roc_auc of .8499, but we hope that we can still do better. We achieve this value with 15 terms and a degree of 1, making it a lasso model

autoplot(tuning_results)

Next we will finalize our workflow object with the optimal hyperparameter values,

best_hyperparameters <- select_best(tuning_results, metric = "roc_auc")
final_wf <- workflow() %>%
  add_recipe(mars_recipe) %>%
  add_model(mars_mod) %>%
  finalize_workflow(best_hyperparameters)

Then we fit our final workflow object across the full training set data

final_fit <- final_wf %>%
  fit(data = mars_train)
final_fit %>%
  extract_fit_parsnip() %>%
  vip()

This model shows our most important features to be tenure, total charges, and payment method.

final_fit %>%
  predict(mars_test) %>%
  bind_cols(mars_test %>% select(Status)) %>%
  conf_mat(truth = Status, estimate = .pred_class)
##           Truth
## Prediction Current Left
##    Current    1382  245
##    Left        158  312

Decision Trees

dt_mod <- decision_tree(mode = 'classification') %>%
  set_engine("rpart")
model_recipe <- recipe(Status ~ ., data = log_train)
dt_fit <- workflow() %>%
  add_recipe(model_recipe) %>%
  add_model(dt_mod) %>%
  fit(data = log_train)
rpart.plot::rpart.plot(dt_fit$fit$fit$fit)

The root node of the decision tree is contract = 1 year or 2 years

set.seed(123)
kfold <- vfold_cv(log_train, v = 5)
dt_results <- fit_resamples(dt_mod, model_recipe, kfold)
collect_metrics(dt_results)
## # A tibble: 2 × 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy binary     0.785     5 0.00552 Preprocessor1_Model1
## 2 roc_auc  binary     0.803     5 0.00777 Preprocessor1_Model1

The decision tree gave us a roc_auc of .8033, which is not an optimal value. With further tuning, we can hope to create an even better model.

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         15    40 roc_auc binary     0.814     5 0.00780
## 2    0.0000000178         15    40 roc_auc binary     0.814     5 0.00780
## 3    0.00000316           15    40 roc_auc binary     0.814     5 0.00780
## 4    0.0000000001         11    40 roc_auc binary     0.814     5 0.00827
## 5    0.0000000178         11    40 roc_auc binary     0.814     5 0.00827
## # ℹ 1 more variable: .config <chr>

We now got the roc_auc up to .8142, which is better, but still not the optimal model.

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 = log_train)
dt_final_fit %>%
  extract_fit_parsnip() %>%
  vip(20)

We find that with decision trees, the most important features are tenure, contract, and total charges.

dt_final_fit %>%
   predict(log_test) %>%
   bind_cols(log_test %>% select(Status)) %>%
   conf_mat(truth = Status, estimate = .pred_class)
##           Truth
## Prediction Current Left
##    Current    1346  275
##    Left        194  282

Finally, we will look at a Random Forest

set.seed(123)
rf_split <- initial_split(retention, prop = 0.7, strata = Status)
rf_train <- training(rf_split)
rf_test <- testing(rf_split)

rf_recipe <- recipe(Status ~ ., data = rf_train)

rf_mod <- rand_forest(mode = "classification") %>%
  set_engine("ranger")

set.seed(123)
rf_kfold <- vfold_cv(rf_train, v = 5)

rf_results <- fit_resamples(rf_mod, rf_recipe, rf_kfold)

collect_metrics(rf_results)
## # A tibble: 2 × 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy binary     0.798     5 0.00509 Preprocessor1_Model1
## 2 roc_auc  binary     0.837     5 0.0103  Preprocessor1_Model1

This model gives us a roc_auc of .8368, but with hyperparameter tuning we hope to do even better.

rf_mod <- rand_forest(
  mode = "classification",
  trees = tune(),
  mtry = tune(),
  min_n = tune()) %>%
  set_engine("ranger", importance = "impurity")

rf_hyper_grid <- grid_regular(
  trees(range = c(50, 800)),
  mtry(range = c(2, 50)),
  min_n(range = c(1, 20)),
  levels = 5)

set.seed(123)
rf_results <- tune_grid(rf_mod, rf_recipe, resamples = rf_kfold, grid = rf_hyper_grid)

show_best(rf_results, metric = "roc_auc")
## # A tibble: 5 × 9
##    mtry trees min_n .metric .estimator  mean     n std_err .config              
##   <int> <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
## 1     2   612    20 roc_auc binary     0.844     5 0.0103  Preprocessor1_Model1…
## 2     2   800    20 roc_auc binary     0.843     5 0.0103  Preprocessor1_Model1…
## 3     2   425    15 roc_auc binary     0.843     5 0.00984 Preprocessor1_Model0…
## 4     2   425    20 roc_auc binary     0.843     5 0.0103  Preprocessor1_Model1…
## 5     2   237    20 roc_auc binary     0.843     5 0.0107  Preprocessor1_Model1…

With hyperparameter tuning, we get the roc_auc to .8435

rf_best_hyperparameters <- select_best(rf_results, metric = "roc_auc")
final_rf_wf <- workflow() %>%
  add_recipe(rf_recipe) %>%
  add_model(rf_mod) %>%
  finalize_workflow(rf_best_hyperparameters)
# fit final workflow object
rf_final_fit <- final_rf_wf %>%
  fit(data = rf_train)
# plot feature importance
rf_final_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 20)

With our random forest model, we find the most important features to be tenure, total charges, and contract

rf_final_fit %>%
   predict(rf_test) %>%
   bind_cols(rf_test %>% select(Status)) %>%
   conf_mat(truth = Status, estimate = .pred_class)
##           Truth
## Prediction Current Left
##    Current    1393  276
##    Left        147  281

Conclusion:

After testing four models, we have found the most accurate model to be the MARS model. It had the highest AUC of .850 when it had 15 terms and a degree of 1. With the MARS model, we found the most important features to be tenure, total charges, and payment method electronic check. Tenure and total charges were also 2 of the top features in our random forest, decision tree, and logistic regression model. The MARS model was also performed well with its confusion matrix. It performed with an accuracy of about 81%.

In order to better retain customers, Regork should focus on on the most influential factors in whether a customer stayed or left, which were tenure, total charges, payment method, and contract. To do this, they could create an incentive program where customers are rewarded for signing longer phone contracts, and to make sure that the contracts do not have many slight up-charges or random monthly charges. A long term contract where the fees are clearly listed and built into the monthly cost should improve customer satisfaction and lead to them staying longer. This contract should accept all payment methods, which again will increase customer satisfaction. Because the customers are signing longer contracts to begin with, their tenure with Regork will increase, and brand loyalty may be produced, and as tenure is the most important factor, getting customers to sign long contracts should be the most focused on factor.

To reward customers for their loyalty, Regork should offer customers incentives such as a discount on their phone plan if they sign a contract over 1 to 2 years in length, offering major discounts or even a free device if they reach a certain amount of time with Regork, or offering discounts when they sign up for multiple phones at a time.

One limitation to our analysis was having a dataset of only ~ 7000 variables. Having a larger dataset to sample from would have allowed for a more accurate model. Our data also had no information on time or the date it is from, so we cannot access trends in the data over time. Another limitation is in our machine learning, in that we did not test every possible model, and we could have tested more of a range in values. Unfortunately with machine learning, there are so many models and so many value ranges to test, so it is very difficult to find the definitive number 1 model. Despite this, we are still very happy with our model and believe that it is very accurate.