Week 06 - BankFlip
Alison
2024-02-28
rm(list=ls())
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(rpart)
library(rpart.plot)
library(randomForest)## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
##
## The following object is masked from 'package:dplyr':
##
## combine
##
## The following object is masked from 'package:ggplot2':
##
## marginPart 01: The data
https://www.kaggle.com/shrutimechlearn/churn-modelling
This data set contains details of a bank’s customers and the target variable is a binary variable reflecting the fact whether the customer left the bank (closed their account) or they continue to be a customer.
- RowNumber Row Numbers from 1 to 10000
- CustomerId Unique Ids for bank customer identification
- Surname Customer’s last name
- CreditScore Credit score of the customer
- Geography The country from which the customer belongs
- Gender Male or Female
- Age Age of the customer
- Tenure Number of years for which the customer has been with the bank
- Balance Bank balance of the customer
- NumOfProducts Number of bank products the customer is utilising
- HasCrCard Binary Flag for whether the customer holds a credit card with the bank or not
- IsActiveMember Binary Flag for whether the customer is an active member with the bank or not
- EstimatedSalary Estimated salary of the customer in Dollars
- Exited Binary flag 1 if the customer closed account with bank and 0 if the customer is retained
churn <- read_csv('./Churn_Modelling.csv')## Rows: 10000 Columns: 14
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): Surname, Geography, Gender
## dbl (11): RowNumber, CustomerId, CreditScore, Age, Tenure, Balance, NumOfPro...
##
## ℹ 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.glimpse(churn)## Rows: 10,000
## Columns: 14
## $ RowNumber <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,…
## $ CustomerId <dbl> 15634602, 15647311, 15619304, 15701354, 15737888, 1557…
## $ Surname <chr> "Hargrave", "Hill", "Onio", "Boni", "Mitchell", "Chu",…
## $ CreditScore <dbl> 619, 608, 502, 699, 850, 645, 822, 376, 501, 684, 528,…
## $ Geography <chr> "France", "Spain", "France", "France", "Spain", "Spain…
## $ Gender <chr> "Female", "Female", "Female", "Female", "Female", "Mal…
## $ Age <dbl> 42, 41, 42, 39, 43, 44, 50, 29, 44, 27, 31, 24, 34, 25…
## $ Tenure <dbl> 2, 1, 8, 1, 2, 8, 7, 4, 4, 2, 6, 3, 10, 5, 7, 3, 1, 9,…
## $ Balance <dbl> 0.00, 83807.86, 159660.80, 0.00, 125510.82, 113755.78,…
## $ NumOfProducts <dbl> 1, 1, 3, 2, 1, 2, 2, 4, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, …
## $ HasCrCard <dbl> 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, …
## $ IsActiveMember <dbl> 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, …
## $ EstimatedSalary <dbl> 101348.88, 112542.58, 113931.57, 93826.63, 79084.10, 1…
## $ Exited <dbl> 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, …Create a new outcome variable with labels.
churn <- churn %>%
mutate(exit = factor(Exited, labels = c('Stayed', 'Left'))) %>%
mutate(Gender = as.factor(Gender))Test it real quick to make sure it makes sense.
table(churn$exit, churn$Exited)##
## 0 1
## Stayed 7963 0
## Left 0 2037The outcome variable is unbalanced. But not to a huge amount.
ggplot(churn, aes(x = exit)) +
geom_bar() +
scale_x_discrete("Customer Status") +
theme_minimal()
We want to predict if customers leave the bank. What if we just create a model that just naively states that everyone stays and never leaves? How accurate would we be?
library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftthe_model_results <- factor(rep(0, nrow(churn)),
levels = c(0,1),
labels = c('Stayed', 'Left'))
the_actual_data <- churn$exit
confusionMatrix(the_model_results, the_actual_data)## Confusion Matrix and Statistics
##
## Reference
## Prediction Stayed Left
## Stayed 7963 2037
## Left 0 0
##
## Accuracy : 0.7963
## 95% CI : (0.7883, 0.8042)
## No Information Rate : 0.7963
## P-Value [Acc > NIR] : 0.5059
##
## Kappa : 0
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 1.0000
## Specificity : 0.0000
## Pos Pred Value : 0.7963
## Neg Pred Value : NaN
## Prevalence : 0.7963
## Detection Rate : 0.7963
## Detection Prevalence : 1.0000
## Balanced Accuracy : 0.5000
##
## 'Positive' Class : Stayed
## Look at the accuracy value - why do you think the accuracy is so high? Any prediction model we build needs to be better than this!
Look at the specificity and sensitivity.
Specificity measures how well the model correctly identifies negative cases (the majority class, “Stayed”). As our naive model predicts that everyone stays, the model cannot correctly identify the minority class at all - it is not learning “Left” at all.
Sensitivity is the true positive rate. Here it is 1, as the model
hasn’t seen any minority cases (“Left”), hence is treating
stayed as the positive class.
Fit a basic decision tree to the churn data.
churn_tree1 <- rpart(exit ~ Age + Gender + Tenure + Balance +
NumOfProducts + HasCrCard + IsActiveMember +
EstimatedSalary + CreditScore,
data = churn)
rpart.plot(churn_tree1)
This is better, but because leavers are relatively rare, the accuracy really doesn’t go up too much. The specificity got way better.
p <- predict(churn_tree1)[,'Left']
the_model_results <- factor(ifelse(p > 0.5, 1, 0),
levels = c(0, 1),
labels = c('Stayed', 'Left'))
the_actual_data <- churn$exit
confusionMatrix(the_model_results, the_actual_data)## Confusion Matrix and Statistics
##
## Reference
## Prediction Stayed Left
## Stayed 7641 1132
## Left 322 905
##
## Accuracy : 0.8546
## 95% CI : (0.8475, 0.8615)
## No Information Rate : 0.7963
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.474
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9596
## Specificity : 0.4443
## Pos Pred Value : 0.8710
## Neg Pred Value : 0.7376
## Prevalence : 0.7963
## Detection Rate : 0.7641
## Detection Prevalence : 0.8773
## Balanced Accuracy : 0.7019
##
## 'Positive' Class : Stayed
## Let’s make a really complex model.
rpart.control()
is used to control the complexity parameter (cp) when building a decision tree using therpart package`.
What does
cp = 0do? The cp parameter stands forcomplexity parameterin therpart.control()function. It prevents pruning by setting cp = 0, which means the tree continues growing until all possible splits are made. In most cases, this leads to an overly complex tree that is overfitted to the data.
churn_tree2 <- rpart(exit ~ Age + Gender + Tenure + Balance +
NumOfProducts + HasCrCard + IsActiveMember +
EstimatedSalary + CreditScore,
data = churn,
control = rpart.control(cp = 0))
rpart.plot(churn_tree2)## Warning: labs do not fit even at cex 0.15, there may be some overplotting
Did we gain anything from increasing the complexity?
p <- predict(churn_tree2)[,'Left']
the_model_results <- factor(ifelse(p > 0.5, 1, 0),
levels = c(0, 1),
labels = c('Stayed', 'Left'))
the_actual_data <- churn$exit
confusionMatrix(the_model_results, the_actual_data)## Confusion Matrix and Statistics
##
## Reference
## Prediction Stayed Left
## Stayed 7631 773
## Left 332 1264
##
## Accuracy : 0.8895
## 95% CI : (0.8832, 0.8956)
## No Information Rate : 0.7963
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6295
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9583
## Specificity : 0.6205
## Pos Pred Value : 0.9080
## Neg Pred Value : 0.7920
## Prevalence : 0.7963
## Detection Rate : 0.7631
## Detection Prevalence : 0.8404
## Balanced Accuracy : 0.7894
##
## 'Positive' Class : Stayed
## Part 02: ROC Curves and AUC
The output of the decision tree is a probability, not a 1 or 0. Using 0.50 is a natural threshold to try, but if you try different thresholds can you get a better output?
Let’s try it at 0.6:
p <- predict(churn_tree2)[,'Left']
the_model_results <- factor(ifelse(p > 0.6, 1, 0),
levels = c(0, 1),
labels = c('Stayed', 'Left'))
the_actual_data <- churn$exit
cm <- confusionMatrix(the_model_results, the_actual_data)We can plot and visualise - you can see the trade-off between the True Positive rate (sensitivity) and the True Negative rate (specificity) as you change the threshold:
library(ggtext)
## threshold for predicting positive outcome
p_threshold <- 0.5
## model predictions
churn_pred <- churn %>%
mutate(p_exit = predict(churn_tree2)[,'Left']) %>%
mutate(pred_exit = factor(ifelse(p_exit > p_threshold, 1, 0),
levels = c(0, 1),
labels = c('Stayed', 'Left')))
## confusion matrix
cm <- confusionMatrix(churn_pred$pred_exit,
churn_pred$exit)
## create some annotation messages
msg_tp <- sprintf("TP:<br> %d (%1.2f)", cm$table['Left', 'Left'], cm$table['Left', 'Left'] / nrow(churn))
msg_fp <- sprintf("FP:<br> %d (%1.2f)", cm$table['Left', 'Stayed'], cm$table['Left', 'Stayed'] / nrow(churn))
msg_tn <- sprintf("TN:<br> %d (%1.2f)", cm$table['Stayed', 'Stayed'], cm$table['Stayed', 'Stayed'] / nrow(churn))
msg_fn <- sprintf("FN:<br> %d (%1.2f)", cm$table['Stayed', 'Left'], cm$table['Stayed', 'Left'] / nrow(churn))
## create the plot
pal <- c("Left" = "#FF4500", "Stayed" = "#008B45")
ggplot(churn_pred, aes(x = p_exit, y = exit, color = pred_exit)) +
geom_point(position = position_jitter(height = 0.3),
alpha = 0.5) +
scale_color_manual("Predicted\nOutcome",
values = pal,
labels = c(Left = "<span style='color:#FF4500'>Left</span>",
Stayed = "<span style='color:#008B45'>Stayed</span>")) +
annotate(geom = "richtext", fill = NA, label.color = NA,
x = p_threshold, y = 1.4, hjust = "left",
label = sprintf("Threshold:<br>p = <b>%1.2f</b>", p_threshold)) +
annotate(geom = 'segment', x = p_threshold, xend = p_threshold,
y = 2.4, yend = 0.6, linetype = 2) +
annotate(geom = "richtext", fill = NA, label.color = NA,
x = p_threshold - 0.02, y = 2.4, hjust = "right",
label = "Predicted:<br><b style='color:#008B45'>Stayed</b>") +
annotate(geom = "richtext", fill = NA, label.color = NA,
x = p_threshold + 0.02, y = 2.4, hjust = "left",
label = "Predicted:<br><b style='color:#FF4500'>Left</b>") +
annotate(geom = "richtext", fill = NA, label.color = NA,
hjust = "left",
x = p_threshold + (1 - p_threshold)/2,
y = 1.6, label = msg_tp) +
annotate(geom = "richtext", fill = NA, label.color = NA,
hjust = "left",
x = p_threshold/2,
y = 0.6, label = msg_tn) +
annotate(geom = "richtext", fill = NA, label.color = NA,
hjust = "left",
x = p_threshold + (1 - p_threshold)/2,
y = 0.6, label = msg_fp) +
annotate(geom = "richtext", fill = NA, label.color = NA,
hjust = "left",
x = p_threshold/2,
y = 1.6, label = msg_fn) +
annotate(geom = "richtext", fill = NA, label.color = NA,
hjust = "right",
x = 0,
y = 2, label = "Observed:<br><b>Left</b>") +
annotate(geom = "richtext", fill = NA, label.color = NA,
hjust = "right",
x = 0,
y = 1, label = "Observed:<br><b>Stayed</b>") +
scale_x_continuous(position = "top", breaks = seq(0, 1, length.out = 5)) +
xlab("Predicted Probability of Leaving") +
theme_minimal() +
theme(legend.position = 'none', axis.text = element_markdown(),
axis.text.y = element_blank(), axis.title.y = element_blank()) +
coord_cartesian(xlim = c(-0.2, 1.2))
Let’s try all the thresholds and see what value we get for them all.
p <- predict(churn_tree2)[,'Left'] #predicted probabilities of 'left'
ks <- seq(0,1,length.out = 100)
all_threshold <- lapply(ks, function(k) {
the_model_results <- factor(ifelse(p > k, 1, 0),
levels = c(0, 1),
labels = c('Stayed', 'Left'))
the_actual_data <- churn$exit
cm <- confusionMatrix(the_model_results, the_actual_data)
tibble(accuracy = cm$overall['Accuracy'],
sens = cm$byClass['Sensitivity'],
spec = cm$byClass['Specificity'])
}) %>% bind_rows() %>%
mutate(threshold = ks)Plot the accuracy, sensitivity, and specificity for each threshold.
ggplot(gather(all_threshold, stat, value, -threshold),
aes(x = threshold, y = value, group = stat, color = stat)) +
geom_line() +
scale_x_continuous(breaks = seq(0,1,by =0.1)) +
theme_classic()## Warning: attributes are not identical across measure variables; they will be
## dropped
Another way of plotting this data is with the x-axis as the Specificity and the y-axis as the Sensitivity. Actually, we reverse the x-axis, as a ROC plot actually plots sensitivity (True Positive rate) against 1-Specificity (False Positive rate)
ggplot(data = all_threshold, aes(x = spec, y = sens, color = threshold)) +
geom_point() +
scale_color_viridis_c('Threshold') +
geom_abline(intercept = 1, slope = 1, linetype = 2) +
scale_x_reverse() +
coord_equal() +
labs(y = 'Sensitivity', x = 'Specificity') +
theme_minimal()
pROC is a package just for ROC curves like this. Use
asp=NA to control scaling
library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, varpred <- tibble(p_exit = as.numeric(predict(churn_tree2)[,'Left']),
exit = churn$Exited)
churn_roc1 <- roc(data = pred, response = exit, predictor = p_exit)## Setting levels: control = 0, case = 1## Setting direction: controls < casesplot(churn_roc1, asp=NA)
The AUC or Area Under the Curve is the measure of predictive performance for the model.
Interpretation:
AUC = 1.0 Perfect classification.
AUC = 0.5 Random guessing (bad model).
AUC < 0.5 Worse than random (something is wrong).
auc(churn_roc1)## Area under the curve: 0.8946ggroc() creates a ROC curve with more control over
aesthetics:
Adds a diagonal reference line (geom_abline()) to compare against random guessing.
Annotates the AUC score on the plot.
Uses coord_equal() to ensure correct scaling.
ggroc(churn_roc1) +
geom_abline(intercept = 1, slope = 1, linetype = 2) +
annotate('text', x = 0.75, y = 0.5,
label = sprintf('AUC: %1.2f', churn_roc1$auc)) +
labs(y = 'Sensitivity', x = 'Specificity') +
coord_equal() +
theme_minimal()
We can ask the ROC object what the best threshold would be to maximize both sensitivity and specificity. It’s considerably different from 0.5
coords(churn_roc1, 'best')## threshold specificity sensitivity
## 1 0.2034483 0.8867261 0.7506136We can add it to the plot:
best_coords <- coords(churn_roc1, 'best')
ggroc(churn_roc1) +
geom_abline(intercept = 1, slope = 1, linetype = 2) +
annotate('text', x = 0.75, y = 0.5,
label = sprintf('AUC: %1.2f', churn_roc1$auc)) +
annotate('point', x = best_coords$specificity,
y = best_coords$sensitivity) +
labs(y = 'Sensitivity', x = 'Specificity') +
coord_equal() +
theme_minimal()
This code is computing a confusion matrix using a custom threshold (best_coords$threshold), which is obtained from the ROC analysis. The sensitivity and specificity are maximised.
p <- predict(churn_tree2)[,'Left']
the_model_results <- factor(ifelse(p > best_coords$threshold, 1, 0),
levels = c(0, 1),
labels = c('Stayed', 'Left'))
the_actual_data <- churn$exit
confusionMatrix(the_model_results, the_actual_data)## Confusion Matrix and Statistics
##
## Reference
## Prediction Stayed Left
## Stayed 7061 508
## Left 902 1529
##
## Accuracy : 0.859
## 95% CI : (0.852, 0.8658)
## No Information Rate : 0.7963
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.5945
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.8867
## Specificity : 0.7506
## Pos Pred Value : 0.9329
## Neg Pred Value : 0.6290
## Prevalence : 0.7963
## Detection Rate : 0.7061
## Detection Prevalence : 0.7569
## Balanced Accuracy : 0.8187
##
## 'Positive' Class : Stayed
## We can present our outcomes in different ways. Why might you want to do this? Here are a couple of examples:
Lift curve
A Lift Curve is a powerful tool for evaluating how well a predictive model identifies churners (exits) compared to random guessing. It provides actionable insights into model effectiveness.
X-axis (perc_pop): The proportion of the population examined (from highest to lowest predicted exit probability). Y-axis (Lift): How much better the model is at identifying churners than a random selection of customers.
A lift of 1.0 means the model is performing no better than random guessing. A lift > 1.0 means the model is identifying exits better than random. The higher the lift, the better the model’s predictive power.
# Define population cutoffs
top_x_percentages <- c(0.05, 0.10, 0.20, 0.50) # Top 5%, 10%, 20%, 50%
# Compute ranked lift curve
nr <- nrow(pred)
pred2 <- pred %>%
arrange(desc(p_exit)) %>% # Sort by predicted churn probability
mutate(rn = row_number()) %>% # Assign rank
mutate(perc_pop = rn / nr) %>% # Compute % of the population seen
mutate(exit_sum = cumsum(exit)) %>% # Cumulative sum of actual exits
mutate(exit_p = exit_sum / sum(exit)) %>% # Percentage of total exits identified
mutate(lift = exit_p / perc_pop) # Compute lift
# Compute lift at different top X% cutoffs
top_x_lift <- pred2 %>%
filter(perc_pop %in% top_x_percentages) %>%
select(perc_pop, lift)
# Print results
print(top_x_lift)## # A tibble: 4 × 2
## perc_pop lift
## <dbl> <dbl>
## 1 0.05 4.63
## 2 0.1 4.36
## 3 0.2 3.49
## 4 0.5 1.83e.g. A Lift of 3.5 for the Top 5% means our model is identifying 3.5× more churners than random selection at that level. As we increase the population size (examining more customers), Lift decreases because we’re including lower-risk customers. Higher Lift in early cutoffs means a strong predictive model (good at ranking high-risk churners early).
Let’s visualise it:
ggplot(pred2, aes(x = perc_pop, y = lift)) +
geom_line(color = "blue", size = 1.2) +
geom_point(data = top_x_lift, aes(x = perc_pop, y = lift), color = "red", size = 3) +
scale_x_continuous(labels = scales::percent) + # Convert X-axis to percentage
labs(title = "Lift Across Different Population Cutoffs",
x = "Proportion of Population Examined",
y = "Lift Value") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Density plots
The plot below shows how predicted churn probabilities are distributed for customers who actually stayed (blue) vs. customers who actually left (red).
Separation Between Stayed vs. Left Customers:
If red (Left) and blue (Stayed) are well separated, the model is good at distinguishing churners. If they overlap significantly, the model struggles to differentiate them.
Impact of the Decision Threshold (Vertical Line):
The dashed black threshold line represents the cutoff for
classification. Customers to the right of this line are predicted as
“Left” (churners), while those to the left are predicted as
“Stayed.”
Misclassification Issues:
Blue customers are to the right of the threshold are False Positives (customers predicted to leave but actually stayed). Red customers to the left of the threshold are False Negatives (missed churners who should have been classified as leaving).
Implications for Decision-Making:
Lowering the threshold (moving the line left) catches more churners (higher recall) but increases false alarms (lower precision). Raising the threshold (moving the line right) reduces false alarms but misses more actual churners (lower recall). A well-placed threshold balances false positives and false negatives, depending on business needs.
When to Adjust the Threshold?
If the cost of losing a customer is high, use a lower threshold (higher recall, even if it causes more false positives).
If false churn alerts are costly (e.g., retention campaigns are expensive), use a higher threshold (reducing false positives).
threshold_value <- best_coords[1,1]
ggplot(pred2, aes(x = perc_pop, fill = factor(exit))) +
geom_density(alpha = 0.4) +
scale_fill_manual(values = c("red", "blue"), labels = c("Left", "Stayed")) + # Ensure correct mapping
geom_vline(xintercept = threshold_value, linetype = "dashed", color = "black", linewidth = 1) +
annotate("text", x = threshold_value, y = 2, label = sprintf("Threshold: %.2f", threshold_value),
color = "black", vjust = -0.5, hjust = -0.1) +
labs(title = "Predicted Churn Probability Distribution",
x = "Predicted Churn Probability", y = "Density", fill = "Actual Outcome") +
theme_minimal()
Part 03: Cross Validation
Create cross validation samples. We will actually just use one as a hold-out set here
churn$set <- sample(1:10, nrow(churn), replace = T)
churn_train <- filter(churn, set != 1)
churn_test <- filter(churn, set == 1)
dim(churn)## [1] 10000 16dim(churn_train)## [1] 9056 16dim(churn_test)## [1] 944 16Decision trees
churn_tree3 <- rpart(exit ~ Age + Gender + Tenure + Balance +
NumOfProducts + HasCrCard + IsActiveMember +
EstimatedSalary + CreditScore,
data = churn_train,
control = rpart.control(cp = 0))
rpart.plot(churn_tree3)## Warning: labs do not fit even at cex 0.15, there may be some overplotting
Use this model (churn_tree3), that was created using the
training data (churn_test) and test it on the “holdout
data” (churn_test).
predict(churn_tree3, newdata = churn_test) %>% head## Stayed Left
## 1 0.9795699 0.02043011
## 2 0.8214286 0.17857143
## 3 0.9618768 0.03812317
## 4 0.2083333 0.79166667
## 5 0.2500000 0.75000000
## 6 0.1538462 0.84615385Checking out the ROC and AUC for the training set.
pred <- tibble(p_exit = predict(churn_tree3, newdata = churn_train)[,'Left'],
exit = churn_train$Exited)
churn_roc <- roc(data = pred, response = exit, predictor = p_exit)## Setting levels: control = 0, case = 1## Setting direction: controls < caseschurn_roc$auc## Area under the curve: 0.889plot(churn_roc)
Checking out the ROC and AUC for the test set.
pred <- tibble(p_exit = predict(churn_tree3, newdata = churn_test)[,'Left'],
exit = churn_test$Exited)
churn_roc <- roc(data = pred, response = exit, predictor = p_exit)## Setting levels: control = 0, case = 1## Setting direction: controls < caseschurn_roc$auc## Area under the curve: 0.7999plot(churn_roc, asp=NA)
The complexity parameter (cp) controls the
minimum improvement necessary in order to create a new split of the
data. The smaller the value, the smaller the improvement needs to be,
and as a result you see more splits of the data. It turns into lots of
splits to get tiny improvements. But this can result in overfitting to
outlies or idiosyncracies of the training set and won’t generalize to
new data.
Let’s test how different values of the complexity parameter (cp) impact training and test accuracy in a decision tree model (rpart). This helps identify the best trade-off between model complexity and generalization.
What do you see? What would be the best cp?
churn_form <- exit ~ Age + Gender + Tenure + Balance +
NumOfProducts + HasCrCard + IsActiveMember +
EstimatedSalary + CreditScore
cp_seq <- c(0.05, 0.02, 0.01, 0.005, 0.001, 0.0005, 0.0001, 0.00005, 0.00001)
train_test_accuracy <- function(k) {
churn_tree4 <- rpart(churn_form, data = churn_train,
control = rpart.control(cp = k))
accuracy_train <- predict(churn_tree4, newdata = churn_train)[,'Left'] %>%
{ ifelse(. > 0.5, 1, 0)} %>%
{ . == churn_train$Exited} %>%
{ sum(.) / length(.)}
accuracy_test <- predict(churn_tree4, newdata = churn_test)[,'Left'] %>%
{ ifelse(. > 0.5, 1, 0)} %>%
{ . == churn_test$Exited} %>%
{ sum(.) / length(.)}
tibble(accuracy_train = accuracy_train, accuracy_test = accuracy_test)
}
d <- lapply(cp_seq, train_test_accuracy) %>%
bind_rows %>%
mutate(cp = cp_seq)
ggplot(gather(d, set, value, -cp),
aes(x = cp, y = value, group = set, linetype = set)) +
geom_line() +
scale_x_log10(breaks = cp_seq, labels = sprintf('%1.5f', cp_seq)) +
# scale_y_continuous(limits = c(0,1)) + # note the y-scale
labs(x = 'Complexity Parameter\n<--- more complex less complex --->',
y = 'Accuracy') +
theme_classic()
We could experiment with “pruning” our tree.
CP Table: Each row shows a potential pruning point. xerror
(Cross-validation error): Helps select the optimal subtree.
printcp(churn_tree3)##
## Classification tree:
## rpart(formula = exit ~ Age + Gender + Tenure + Balance + NumOfProducts +
## HasCrCard + IsActiveMember + EstimatedSalary + CreditScore,
## data = churn_train, control = rpart.control(cp = 0))
##
## Variables actually used in tree construction:
## [1] Age Balance CreditScore EstimatedSalary
## [5] Gender HasCrCard IsActiveMember NumOfProducts
## [9] Tenure
##
## Root node error: 1838/9056 = 0.20296
##
## n= 9056
##
## CP nsplit rel error xerror xstd
## 1 0.06610446 0 1.00000 1.00000 0.020824
## 2 0.03536453 2 0.86779 0.89010 0.019920
## 3 0.03482046 4 0.79706 0.84603 0.019526
## 4 0.03046790 5 0.76224 0.78618 0.018960
## 5 0.02557127 6 0.73177 0.73993 0.018496
## 6 0.00380849 7 0.70620 0.71763 0.018264
## 7 0.00353645 10 0.69314 0.72579 0.018350
## 8 0.00272035 12 0.68607 0.71545 0.018241
## 9 0.00258433 13 0.68335 0.71872 0.018275
## 10 0.00163221 17 0.67301 0.73286 0.018423
## 11 0.00136017 22 0.66485 0.73885 0.018485
## 12 0.00126950 34 0.63874 0.73776 0.018474
## 13 0.00119695 39 0.63112 0.73830 0.018480
## 14 0.00108814 50 0.61589 0.73830 0.018480
## 15 0.00072543 65 0.59848 0.75245 0.018624
## 16 0.00054407 90 0.57889 0.75571 0.018657
## 17 0.00043526 108 0.56746 0.77911 0.018891
## 18 0.00040805 119 0.56257 0.78509 0.018949
## 19 0.00036271 125 0.55930 0.79489 0.019045
## 20 0.00027203 136 0.55441 0.80359 0.019129
## 21 0.00020403 140 0.55332 0.82318 0.019314
## 22 0.00018136 149 0.55114 0.83351 0.019411
## 23 0.00013602 159 0.54897 0.83841 0.019456
## 24 0.00000000 175 0.54570 0.85092 0.019571Find the best pruning point (cp) with the lowest cross-validation error (xerror). If multiple cp values are close, chooses the simpler tree.
best_cp <- churn_tree3$cptable[which.min(churn_tree3$cptable[,"xerror"]), "CP"]
print(best_cp)## [1] 0.002720348Prune the tree! Remove splits that don’t significantly improve predictive power. Keep the best balance between complexity and generalization.
churn_tree_pruned <- prune(churn_tree3, cp = best_cp)Visualise the pruned tree
rpart.plot(churn_tree_pruned)
Random Forest
churn_rf <- randomForest(exit ~ Age + Gender + Tenure + Balance +
NumOfProducts + HasCrCard + IsActiveMember +
EstimatedSalary + CreditScore,
data = churn_train)p_train <- predict(churn_rf, type = 'prob')[,'Left']
the_model_results <- factor(ifelse(p_train > 0.45, 1, 0),
levels = c(0, 1),
labels = c('Stayed', 'Left'))
the_actual_data <- churn_train$exit
confusionMatrix(the_model_results, the_actual_data)## Confusion Matrix and Statistics
##
## Reference
## Prediction Stayed Left
## Stayed 6860 976
## Left 358 862
##
## Accuracy : 0.8527
## 95% CI : (0.8452, 0.8599)
## No Information Rate : 0.797
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4795
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9504
## Specificity : 0.4690
## Pos Pred Value : 0.8754
## Neg Pred Value : 0.7066
## Prevalence : 0.7970
## Detection Rate : 0.7575
## Detection Prevalence : 0.8653
## Balanced Accuracy : 0.7097
##
## 'Positive' Class : Stayed
## pred <- tibble(p_exit = p_train, exit = churn_train$Exited)
churn_roc_train <- roc(data = pred, response = exit, predictor = p_exit)## Setting levels: control = 0, case = 1## Setting direction: controls < casesplot(churn_roc, print.auc = T, auc.polygon = T, print.thres = 'best')
ggplot(pred) +
geom_jitter(aes(x = p_exit, y = exit))
coords(churn_roc, 'best', transpose = FALSE)## threshold specificity sensitivity
## 1 0.1197723 0.7395973 0.7386935p_test <- predict(churn_rf, newdata = churn_test, type = 'prob')[,'Left']
pred_test <- tibble(p_exit = p_test,exit = churn_test$Exited)
churn_roc_test <- roc(data = pred_test, response = exit, predictor = p_exit)## Setting levels: control = 0, case = 1## Setting direction: controls < casesplot(churn_roc_test, print.auc = T, auc.polygon = T, print.thres = 'best')
Plot a comparison of test and train.
# Plot Train ROC first (sets the correct scale)
plot(churn_roc_train, col = "black", lwd = 2,
xlim = c(1, 0), ylim = c(0, 1), # Ensure correct axis scaling
asp = NA,
xlab = "Specificity (True Negative Rate)",
ylab = "Sensitivity (True Positive Rate)",
main = "ROC Curve: Train vs Test")
# Add Test ROC (ensure it aligns correctly)
plot(churn_roc_test, add = TRUE, col = "red", lwd = 2)
# Add a legend with reduced font size
legend("bottomright", # Automatic positioning
legend = c(sprintf("Train (AUC=%1.2f)", auc(churn_roc_train)),
sprintf("Test (AUC=%1.2f)", auc(churn_roc_test))),
col = c("black", "red"),
lty = 1, lwd = 2,
cex = 0.8, # Reduce font size (default is 1.0)
pt.cex = 0.5, # Reduce legend size
bty = "n") # Remove box around legend
Part 04: Trying CV with Caret
library(caret)Use cross validation and a range of mtry parameters to define training control parameters.
mtry controls how many features are randomly chosen at
each split in the Random Forest.
We can create a manual grid to test, e.g:
> # Define a tuning grid - > mtry_values <- seq(2, 4) >
tuneGrid = expand.grid(.mtry = mtry_values)
Then set
tuneGrid = tuneGrid.
Alternatively, we have set tuneLength = 5 which
automatically tests 5 different mtry parameters.
train_control <- trainControl(method = 'cv',
number = 10,
savePredictions = 'final',
classProbs = TRUE,
verboseIter = TRUE)
churn_cv <- train(exit ~ Age + Gender + Tenure + Balance +
NumOfProducts + HasCrCard + IsActiveMember +
EstimatedSalary + CreditScore,
data = churn_train,
trControl = train_control,
preProcess = c('center', 'scale'),
tunelength = 5,
method = 'rf')## + Fold01: mtry=2
## - Fold01: mtry=2
## + Fold01: mtry=5
## - Fold01: mtry=5
## + Fold01: mtry=9
## - Fold01: mtry=9
## + Fold02: mtry=2
## - Fold02: mtry=2
## + Fold02: mtry=5
## - Fold02: mtry=5
## + Fold02: mtry=9
## - Fold02: mtry=9
## + Fold03: mtry=2
## - Fold03: mtry=2
## + Fold03: mtry=5
## - Fold03: mtry=5
## + Fold03: mtry=9
## - Fold03: mtry=9
## + Fold04: mtry=2
## - Fold04: mtry=2
## + Fold04: mtry=5
## - Fold04: mtry=5
## + Fold04: mtry=9
## - Fold04: mtry=9
## + Fold05: mtry=2
## - Fold05: mtry=2
## + Fold05: mtry=5
## - Fold05: mtry=5
## + Fold05: mtry=9
## - Fold05: mtry=9
## + Fold06: mtry=2
## - Fold06: mtry=2
## + Fold06: mtry=5
## - Fold06: mtry=5
## + Fold06: mtry=9
## - Fold06: mtry=9
## + Fold07: mtry=2
## - Fold07: mtry=2
## + Fold07: mtry=5
## - Fold07: mtry=5
## + Fold07: mtry=9
## - Fold07: mtry=9
## + Fold08: mtry=2
## - Fold08: mtry=2
## + Fold08: mtry=5
## - Fold08: mtry=5
## + Fold08: mtry=9
## - Fold08: mtry=9
## + Fold09: mtry=2
## - Fold09: mtry=2
## + Fold09: mtry=5
## - Fold09: mtry=5
## + Fold09: mtry=9
## - Fold09: mtry=9
## + Fold10: mtry=2
## - Fold10: mtry=2
## + Fold10: mtry=5
## - Fold10: mtry=5
## + Fold10: mtry=9
## - Fold10: mtry=9
## Aggregating results
## Selecting tuning parameters
## Fitting mtry = 2 on full training set#churn_cv_eval <- evalm(churn_cv)
confusionMatrix(churn_cv$pred$pred, churn_cv$pred$obs)## Confusion Matrix and Statistics
##
## Reference
## Prediction Stayed Left
## Stayed 7030 1107
## Left 188 731
##
## Accuracy : 0.857
## 95% CI : (0.8496, 0.8642)
## No Information Rate : 0.797
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4568
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9740
## Specificity : 0.3977
## Pos Pred Value : 0.8640
## Neg Pred Value : 0.7954
## Prevalence : 0.7970
## Detection Rate : 0.7763
## Detection Prevalence : 0.8985
## Balanced Accuracy : 0.6858
##
## 'Positive' Class : Stayed
## ggplot(churn_cv) +
labs(title = "RF cross-validation results")
The cross-validation chooses the best mtry. We can plot the
cross-validation results and see how accuracy changes across differnt
mtry
rf_cv <- churn_cv$finalModelAfter training with caret::train(), the predictions from all
cross-validation folds are stored in chrun_cv$pred.
This contains:
.mtry: The hyperparameter value for that model.
obs: The actual class labels (Stayed or Left).
pred: The predicted class (Stayed or Left).
Class probabilities for each class (e.g., Stayed, Left).
We can plot all cross-validation predictions together, treating them as one single dataset. It does not show separate ROC curves for each fold. The resulting ROC curve is an aggregate of all folds.
roc_obj <- roc(response = churn_cv$pred$obs,
predictor = churn_cv$pred$Left,
levels = rev(levels(churn_cv$pred$obs)))## Setting direction: controls > casesplot(roc_obj, asp=NA)
auc(roc_obj)## Area under the curve: 0.8432Visualise feature importance (which features contribute most to predictions). Helps interpretability by identifying key factors influencing churn.
varImpPlot(rf_cv)