Bank Customer Churn Prediction

Introduction

Business Context

Customer churn, also known as customer attrition, occurs when customers stop doing business with a company. In the banking industry, this means customers closing their accounts or switching to another bank. Addressing churn is critical as acquiring new customers typically costs 5-25 times more than retaining existing ones.

This analysis focuses on predicting which bank customers are likely to leave in the near future, enabling proactive retention strategies. By identifying at-risk customers before they leave, banks can implement targeted interventions that:

1. Reduce revenue loss from departing customers
2. Lower customer acquisition costs by maintaining a stable customer base
3. Preserve the investment already made in existing customer relationships
4. Gain competitive advantage through improved customer satisfaction and loyalty

Dataset Description

The dataset contains information about 10,000 bank customers and whether they exited (churned) or remained with the bank. It includes:

• Demographic information: Age, gender, geography
• Financial attributes: Credit score, account balance, estimated salary
• Relationship metrics: Tenure (years as customer), number of products, activity level
• Target variable: Whether the customer exited the bank (1 = yes, 0 = no)

Analysis Approach

This analysis follows a comprehensive machine learning workflow:

1. Exploratory Data Analysis: Understanding data characteristics, distributions, and relationships
2. Feature Engineering: Creating meaningful new features for improved model performance
3. Model Development: Building and comparing multiple predictive models:
   • Logistic Regression (baseline)
   • Random Forest
   • XGBoost
   • Gradient Boosting Machine (GBM)
4. Model Evaluation: Assessing performance using appropriate metrics (AUC-ROC, precision, recall)
5. Customer Segmentation: Identifying high-risk customer segments
6. Business Recommendations: Translating findings into actionable strategies

The ultimate goal is to develop a reliable predictive model that helps identify customers at risk of churning, understand the factors driving that risk, and inform effective retention strategies.

Executive Summary

This analysis develops predictive models to identify bank customers most likely to churn, allowing for proactive retention measures. Using machine learning algorithms applied to a dataset of 10,000 customers with demographic, financial, and relationship attributes, we achieved the following key results:

• Churn Rate: 20.4% of customers in the dataset left the bank
• Model Performance: The XGBoost model achieved the highest AUC of 0.866, indicating strong predictive capability
• Key Churn Drivers: Age, number of products, and account activity status emerged as the most influential predictors
• Customer Profiles: Specific high-risk segments were identified, with “Senior Inactive French” customers showing an 81% churn rate

Strategic recommendations include:

• Targeted retention campaigns for high-risk segments

• Product portfolio optimization, focusing on the ideal 2-product offering

• Special attention to inactive customers with high balances (30.3% churn rate)

• Enhanced focus on the German market, which shows higher churn rates than initially expected

By implementing these model-driven retention strategies, the bank could potentially reduce churn by 25-30%, translating to estimated revenue preservation of €3-5 million.

1. Data Exploration

# Set seed for reproducibility
set.seed(123)

# Load data
churn_data <- read.csv("Churn_Modelling.csv")

# Examine data structure
str(churn_data)

## 'data.frame':    10000 obs. of  14 variables:
##  $ RowNumber      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ CustomerId     : int  15634602 15647311 15619304 15701354 15737888 15574012 15592531 15656148 15792365 15592389 ...
##  $ Surname        : chr  "Hargrave" "Hill" "Onio" "Boni" ...
##  $ CreditScore    : int  619 608 502 699 850 645 822 376 501 684 ...
##  $ Geography      : chr  "France" "Spain" "France" "France" ...
##  $ Gender         : chr  "Female" "Female" "Female" "Female" ...
##  $ Age            : int  42 41 42 39 43 44 50 29 44 27 ...
##  $ Tenure         : int  2 1 8 1 2 8 7 4 4 2 ...
##  $ Balance        : num  0 83808 159661 0 125511 ...
##  $ NumOfProducts  : int  1 1 3 2 1 2 2 4 2 1 ...
##  $ HasCrCard      : int  1 0 1 0 1 1 1 1 0 1 ...
##  $ IsActiveMember : int  1 1 0 0 1 0 1 0 1 1 ...
##  $ EstimatedSalary: num  101349 112543 113932 93827 79084 ...
##  $ Exited         : int  1 0 1 0 0 1 0 1 0 0 ...

# Basic summary statistics
summary(churn_data)

##    RowNumber       CustomerId         Surname           CreditScore   
##  Min.   :    1   Min.   :15565701   Length:10000       Min.   :350.0  
##  1st Qu.: 2501   1st Qu.:15628528   Class :character   1st Qu.:584.0  
##  Median : 5000   Median :15690738   Mode  :character   Median :652.0  
##  Mean   : 5000   Mean   :15690941                      Mean   :650.5  
##  3rd Qu.: 7500   3rd Qu.:15753234                      3rd Qu.:718.0  
##  Max.   :10000   Max.   :15815690                      Max.   :850.0  
##   Geography            Gender               Age            Tenure      
##  Length:10000       Length:10000       Min.   :18.00   Min.   : 0.000  
##  Class :character   Class :character   1st Qu.:32.00   1st Qu.: 3.000  
##  Mode  :character   Mode  :character   Median :37.00   Median : 5.000  
##                                        Mean   :38.92   Mean   : 5.013  
##                                        3rd Qu.:44.00   3rd Qu.: 7.000  
##                                        Max.   :92.00   Max.   :10.000  
##     Balance       NumOfProducts    HasCrCard      IsActiveMember  
##  Min.   :     0   Min.   :1.00   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:     0   1st Qu.:1.00   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median : 97199   Median :1.00   Median :1.0000   Median :1.0000  
##  Mean   : 76486   Mean   :1.53   Mean   :0.7055   Mean   :0.5151  
##  3rd Qu.:127644   3rd Qu.:2.00   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :250898   Max.   :4.00   Max.   :1.0000   Max.   :1.0000  
##  EstimatedSalary         Exited      
##  Min.   :    11.58   Min.   :0.0000  
##  1st Qu.: 51002.11   1st Qu.:0.0000  
##  Median :100193.91   Median :0.0000  
##  Mean   :100090.24   Mean   :0.2037  
##  3rd Qu.:149388.25   3rd Qu.:0.0000  
##  Max.   :199992.48   Max.   :1.0000

# More comprehensive data summary
skim(churn_data)

Data summary

Name	churn_data
Number of rows	10000
Number of columns	14
_______________________
Column type frequency:
character	3
numeric	11
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
Surname	0	1	2	23	0	2932	0
Geography	0	1	5	7	0	3	0
Gender	0	1	4	6	0	2	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
RowNumber	0	1	5000.50	2886.90	1.00	2500.75	5000.50	7500.25	10000.0	▇▇▇▇▇
CustomerId	0	1	15690940.57	71936.19	15565701.00	15628528.25	15690738.00	15753233.75	15815690.0	▇▇▇▇▇
CreditScore	0	1	650.53	96.65	350.00	584.00	652.00	718.00	850.0	▁▃▇▇▃
Age	0	1	38.92	10.49	18.00	32.00	37.00	44.00	92.0	▅▇▂▁▁
Tenure	0	1	5.01	2.89	0.00	3.00	5.00	7.00	10.0	▇▆▆▆▅
Balance	0	1	76485.89	62397.41	0.00	0.00	97198.54	127644.24	250898.1	▇▃▇▂▁
NumOfProducts	0	1	1.53	0.58	1.00	1.00	1.00	2.00	4.0	▇▇▁▁▁
HasCrCard	0	1	0.71	0.46	0.00	0.00	1.00	1.00	1.0	▃▁▁▁▇
IsActiveMember	0	1	0.52	0.50	0.00	0.00	1.00	1.00	1.0	▇▁▁▁▇
EstimatedSalary	0	1	100090.24	57510.49	11.58	51002.11	100193.91	149388.25	199992.5	▇▇▇▇▇
Exited	0	1	0.20	0.40	0.00	0.00	0.00	0.00	1.0	▇▁▁▁▂

# Check for missing values
missing_values <- colSums(is.na(churn_data))
print(missing_values)

##       RowNumber      CustomerId         Surname     CreditScore       Geography 
##               0               0               0               0               0 
##          Gender             Age          Tenure         Balance   NumOfProducts 
##               0               0               0               0               0 
##       HasCrCard  IsActiveMember EstimatedSalary          Exited 
##               0               0               0               0

# Check class imbalance
table(churn_data$Exited) %>% 
  prop.table() %>% 
  round(4) * 100

## 
##     0     1 
## 79.63 20.37

The dataset contains information about 10,000 bank customers with no missing values. Key characteristics include:

• Customer Demographics: Ages range from 18-92 years (median: 37), with balanced gender distribution and three countries represented (France, Spain, Germany)
• Financial Attributes: Credit scores average 651 (range: 350-850), account balances average €76,486 with 36% having zero balance
• Banking Relationship: Most customers (65%) have only one product, 71% have credit cards, and 52% are active members
• Churn Rate: 20.4% of customers exited the bank, showing moderate class imbalance

This clean, comprehensive dataset provides a solid foundation for predictive modeling, with the reasonable class imbalance being manageable for our analysis.

1.1 Data Cleaning and Transformation

# Remove unnecessary columns
churn_clean <- churn_data %>%
  select(-RowNumber, -CustomerId, -Surname)

# Convert categorical variables to factors
churn_clean$Geography <- as.factor(churn_clean$Geography)
churn_clean$Gender <- as.factor(churn_clean$Gender)
churn_clean$Exited <- as.factor(churn_clean$Exited)
churn_clean$HasCrCard <- as.factor(churn_clean$HasCrCard)
churn_clean$IsActiveMember <- as.factor(churn_clean$IsActiveMember)

# Check the structure after transformation
str(churn_clean)

## 'data.frame':    10000 obs. of  11 variables:
##  $ CreditScore    : int  619 608 502 699 850 645 822 376 501 684 ...
##  $ Geography      : Factor w/ 3 levels "France","Germany",..: 1 3 1 1 3 3 1 2 1 1 ...
##  $ Gender         : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 2 2 1 2 2 ...
##  $ Age            : int  42 41 42 39 43 44 50 29 44 27 ...
##  $ Tenure         : int  2 1 8 1 2 8 7 4 4 2 ...
##  $ Balance        : num  0 83808 159661 0 125511 ...
##  $ NumOfProducts  : int  1 1 3 2 1 2 2 4 2 1 ...
##  $ HasCrCard      : Factor w/ 2 levels "0","1": 2 1 2 1 2 2 2 2 1 2 ...
##  $ IsActiveMember : Factor w/ 2 levels "0","1": 2 2 1 1 2 1 2 1 2 2 ...
##  $ EstimatedSalary: num  101349 112543 113932 93827 79084 ...
##  $ Exited         : Factor w/ 2 levels "0","1": 2 1 2 1 1 2 1 2 1 1 ...

1.2 Exploratory Data Analysis

1.2.1 Univariate Analysis

# Function to create standardized histogram
create_histogram <- function(data, variable, title, bins = 30, fill_color = "#3498db") {
  ggplot(data, aes_string(x = variable)) +
    geom_histogram(bins = bins, fill = fill_color, color = "white", alpha = 0.8) +
    labs(title = title, x = variable, y = "Count") +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(hjust = 0.5, face = "bold"))
}

# Numeric variables
p1 <- create_histogram(churn_clean, "CreditScore", "Credit Score Distribution")
p2 <- create_histogram(churn_clean, "Age", "Age Distribution", fill_color = "#2ecc71")
p3 <- create_histogram(churn_clean, "Tenure", "Tenure Distribution", fill_color = "#9b59b6")
p4 <- create_histogram(churn_clean, "Balance", "Balance Distribution", fill_color = "#e74c3c")

# Arrange numeric plots
grid.arrange(p1, p2, p3, p4, ncol = 2)

# Function to create standardized bar plot
create_barplot <- function(data, variable, title, fill_color = "#3498db") {
  ggplot(data, aes_string(x = variable)) +
    geom_bar(fill = fill_color, color = "white", alpha = 0.8) +
    labs(title = title, x = variable, y = "Count") +
    theme_minimal(base_size = 12) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.title = element_text(hjust = 0.5, face = "bold")
    )
}

# Categorical variables
p5 <- create_barplot(churn_clean, "Geography", "Customer Geography", fill_color = "#f39c12")
p6 <- create_barplot(churn_clean, "Gender", "Customer Gender", fill_color = "#d35400")
p7 <- create_barplot(churn_clean, "HasCrCard", "Has Credit Card", fill_color = "#27ae60")
p8 <- create_barplot(churn_clean, "IsActiveMember", "Is Active Member", fill_color = "#8e44ad")

# Arrange categorical plots
grid.arrange(p5, p6, p7, p8, ncol = 2)

# NumOfProducts distribution
ggplot(churn_clean, aes(x = factor(NumOfProducts))) +
  geom_bar(fill = "#1abc9c", color = "white", alpha = 0.8) +
  labs(title = "Number of Products Distribution", x = "Number of Products", y = "Count") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Exited (Churn) distribution
ggplot(churn_clean, aes(x = Exited, fill = Exited)) +
  geom_bar(color = "white", alpha = 0.8) +
  scale_fill_manual(values = c("#3498db", "#e74c3c")) +
  labs(title = "Churn Distribution", x = "Exited", y = "Count") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

1.2.2 Bivariate Analysis

# Function to create comparison plots
create_comparison_plot <- function(data, var_name, title) {
  ggplot(data, aes_string(x = "Exited", y = var_name, fill = "Exited")) +
    geom_boxplot(alpha = 0.8) +
    scale_fill_manual(values = c("#3498db", "#e74c3c")) +
    labs(title = title, y = var_name) +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(hjust = 0.5, face = "bold"))
}

# Create comparison boxplots
p1 <- create_comparison_plot(churn_clean, "CreditScore", "Credit Score by Churn Status")
p2 <- create_comparison_plot(churn_clean, "Age", "Age by Churn Status")
p3 <- create_comparison_plot(churn_clean, "Tenure", "Tenure by Churn Status")
p4 <- create_comparison_plot(churn_clean, "Balance", "Balance by Churn Status")

# Arrange comparison plots
grid.arrange(p1, p2, p3, p4, ncol = 2)

# Function to create barplots for categorical comparisons
create_cat_comparison <- function(data, var_name, title) {
  ggplot(data, aes_string(x = var_name, fill = "Exited")) +
    geom_bar(position = "fill") +
    scale_fill_manual(values = c("#3498db", "#e74c3c"), name = "Churn") +
    labs(title = title, y = "Proportion", x = var_name) +
    theme_minimal(base_size = 12) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.title = element_text(hjust = 0.5, face = "bold")
    )
}

# Create categorical comparison plots
p5 <- create_cat_comparison(churn_clean, "Geography", "Churn Rate by Geography")
p6 <- create_cat_comparison(churn_clean, "Gender", "Churn Rate by Gender")
p7 <- create_cat_comparison(churn_clean, "HasCrCard", "Churn Rate by Credit Card Status")
p8 <- create_cat_comparison(churn_clean, "IsActiveMember", "Churn Rate by Activity Status")

# Arrange categorical comparison plots
grid.arrange(p5, p6, p7, p8, ncol = 2)

# NumOfProducts vs Churn
ggplot(churn_clean, aes(x = factor(NumOfProducts), fill = Exited)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("#3498db", "#e74c3c"), name = "Churn") +
  labs(title = "Churn Rate by Number of Products", y = "Proportion", x = "Number of Products") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Age groups vs Churn
churn_clean <- churn_clean %>%
  mutate(age_group = cut(Age, breaks = c(0, 30, 40, 50, 60, 100), 
                          labels = c("Under 30", "30-40", "40-50", "50-60", "60+")))

ggplot(churn_clean, aes(x = age_group, fill = Exited)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("#3498db", "#e74c3c"), name = "Churn") +
  labs(title = "Churn Rate by Age Group", y = "Proportion", x = "Age Group") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Balance groups vs Churn
churn_clean <- churn_clean %>%
  mutate(balance_group = cut(Balance, 
                             breaks = c(-1, 0.01, 50000, 100000, 150000, 250000), 
                             labels = c("0", "0-50K", "50K-100K", "100K-150K", "150K+")))

ggplot(churn_clean, aes(x = balance_group, fill = Exited)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("#3498db", "#e74c3c"), name = "Churn") +
  labs(title = "Churn Rate by Balance Group", y = "Proportion", x = "Balance Group") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

The bivariate analysis reveals key relationships between customer characteristics and churn behavior:

Demographic Factors:

• Age Impact: A striking relationship exists, with the 50-60 age group showing the highest churn rate (~55%) while younger customers (under 30) have the lowest (~5%)

• Geographic Variation: German customers exhibit higher churn rates (~33%) compared to France and Spain (both ~15%)

• Gender Differences: Female customers show moderately higher churn propensity (~25%) compared to males (~15%)

Banking Relationship Indicators:

• Product Ownership: A clear U-shaped pattern emerges where customers with 2 products show the lowest churn rate (~8%), while those with 3-4 products demonstrate dramatically higher churn (80-100%)

• Account Activity: Inactive members are more than twice as likely to churn (~25%) compared to active members (~12%)

• Balance Relationships: Both very low and very high balance groups show elevated churn rates, suggesting a non-linear relationship

These patterns highlight the importance of product portfolio management, engagement strategies, and age-appropriate offerings in customer retention.

1.2.3 Correlation Analysis

# Select numerical variables for correlation
num_vars <- churn_clean %>%
  select(CreditScore, Age, Tenure, Balance, NumOfProducts, EstimatedSalary) %>%
  mutate(Exited = as.numeric(as.character(churn_clean$Exited)))

# Calculate correlation matrix
correlation <- cor(num_vars)

# Plot correlation matrix
corrplot(correlation, 
         method = "circle", 
         type = "upper", 
         order = "hclust",
         addCoef.col = "black", 
         tl.col = "black", 
         tl.srt = 45,
         diag = FALSE,
         title = "Correlation Matrix of Numerical Variables",
         mar = c(0,0,2,0))

The correlation analysis reveals:

• Age and Churn: Moderate positive correlation (0.29), confirming the importance of age in predicting churn
• Balance and Number of Products: Notable negative correlation (-0.30), indicating customers with higher balances tend to have fewer products
• Balance and Churn: Weak positive correlation (0.12), suggesting a modest relationship
• Limited Linear Relationships: Most variables show minimal correlation with each other, indicating independence of predictors
• Credit Score and Churn: Essentially no correlation (0.00), suggesting credit score is not a key driver of churn

These findings suggest that many relationships may be non-linear, supporting our decision to use tree-based models alongside logistic regression in our modeling approach.

1.3 Key Insights from EDA

Our exploratory analysis reveals several crucial patterns and relationships:

1. Churn Rate (20.4%): Represents a significant retention challenge requiring targeted intervention

2. Age as Key Driver: Customers over 50 show dramatically higher churn rates, with the 50-60 age group reaching 55%

3. Product Portfolio Effect: A striking U-shaped relationship where customers with 2 products have optimal retention (8% churn), while those with 3-4 products show extremely high churn rates (80-100%)

4. Geographic Variation: German customers exhibit substantially higher churn rates (~33%) compared to other markets

5. Activity Status Impact: Inactive members have more than double the churn rate of active members

6. Balance-Churn Relationship: Both very low and high-balance customers show elevated churn risk, suggesting different strategies may be needed for various wealth segments

7. Interaction Effects: The combination of factors (age, product count, activity status, geography) amplifies churn risk far beyond individual effects

8. Weak Linear Correlations: Most numerical variables show limited linear correlation, indicating the need for models that can capture non-linear relationships

These insights inform our feature engineering approach and model selection in subsequent analysis.

2. Feature Engineering and Data Preparation

2.1 Feature Engineering

# Create new features
churn_featured <- churn_clean %>%
  mutate(
    # Balance to Salary Ratio (indicates what proportion of annual income is kept in the bank)
    BalanceToSalary = Balance / EstimatedSalary,
    
    # Zero Balance Flag (many customers have exactly zero balance)
    ZeroBalance = ifelse(Balance == 0, 1, 0),
    
    # Products per Year (product adoption rate)
    ProductsPerYear = ifelse(Tenure == 0, NumOfProducts, NumOfProducts / Tenure),
    
    # Age Groups (categorical version of age)
    AgeGroup = cut(Age, breaks = c(0, 30, 40, 50, 60, 100), 
                   labels = c("Under 30", "30-40", "40-50", "50-60", "60+")),
    
    # Tenure Groups
    TenureGroup = cut(Tenure, breaks = c(-1, 2, 5, 8, 10), 
                      labels = c("New (<2)", "Developing (2-5)", "Established (5-8)", "Mature (8+)"))
  )

# Create interaction features using numeric format to avoid factor level issues
churn_featured <- churn_featured %>%
  mutate(
    # Interaction between activity and balance
    InactiveHighBalance = ifelse(as.numeric(as.character(IsActiveMember)) == 0 & Balance > 100000, 1, 0),
    
    # Interaction between age and product count
    SeniorSingleProduct = ifelse(Age > 50 & NumOfProducts == 1, 1, 0),
    
    # Interaction between geography and gender
    FranceFemale = ifelse(Geography == "France" & Gender == "Female", 1, 0)
  )

# Convert new binary features to factors
churn_featured$ZeroBalance <- as.factor(churn_featured$ZeroBalance)
churn_featured$InactiveHighBalance <- as.factor(churn_featured$InactiveHighBalance)
churn_featured$SeniorSingleProduct <- as.factor(churn_featured$SeniorSingleProduct)
churn_featured$FranceFemale <- as.factor(churn_featured$FranceFemale)

# Examine new features
summary(churn_featured[, c("BalanceToSalary", "ZeroBalance", "ProductsPerYear", 
                           "InactiveHighBalance", "SeniorSingleProduct", "FranceFemale")])

##  BalanceToSalary     ZeroBalance ProductsPerYear  InactiveHighBalance
##  Min.   :    0.000   0:6383      Min.   :0.1000   0:7644             
##  1st Qu.:    0.000   1:3617      1st Qu.:0.2000   1:2356             
##  Median :    0.747               Median :0.3333                      
##  Mean   :    3.879               Mean   :0.5085                      
##  3rd Qu.:    1.514               3rd Qu.:0.6667                      
##  Max.   :10614.655               Max.   :4.0000                      
##  SeniorSingleProduct FranceFemale
##  0:9251              0:7739      
##  1: 749              1:2261      
##                                  
##                                  
##                                  
## 

Based on our exploratory findings, we’ve created several targeted features to enhance model performance:

1. Financial Relationship Features:
   • Balance-to-Salary Ratio: Measures deposit concentration relative to income (mean: 3.88)
   • Zero Balance Flag: Identifies the 36.2% of customers with no deposits
2. Engagement Metrics:
   • Products-per-Year: Measures product adoption rate over time (mean: 0.51)
   • Inactive High-Balance: Flags the high-risk segment of inactive customers with substantial balances (23.6%)
3. Demographic Interaction Terms:
   • Senior Single-Product: Identifies older customers with limited product engagement (7.5%)
   • France Female: Captures the gender-geography interaction seen in our EDA (22.6%)
4. Categorical Transformations:
   • Age Groups: Converts continuous age to meaningful life-stage categories
   • Tenure Groups: Transforms customer tenure into relationship maturity segments

These engineered features incorporate domain knowledge and EDA insights to capture relevant patterns and interactions that may improve predictive power.

2.2 Data Preprocessing for Modeling

# One-hot encode categorical variables
cat_vars <- c("Geography", "Gender", "AgeGroup", "TenureGroup")
dummies <- model.matrix(~ Geography + Gender + AgeGroup + TenureGroup - 1, data = churn_featured)

# Prepare numerical data
numeric_data <- churn_featured %>%
  select(CreditScore, Age, Tenure, Balance, NumOfProducts, HasCrCard, 
         IsActiveMember, EstimatedSalary, BalanceToSalary, ZeroBalance, 
         ProductsPerYear, InactiveHighBalance, SeniorSingleProduct, 
         FranceFemale, Exited)

# Convert factor columns to numeric for ML compatibility
numeric_data$HasCrCard <- as.numeric(as.character(numeric_data$HasCrCard))
numeric_data$IsActiveMember <- as.numeric(as.character(numeric_data$IsActiveMember))
numeric_data$ZeroBalance <- as.numeric(as.character(numeric_data$ZeroBalance))
numeric_data$InactiveHighBalance <- as.numeric(as.character(numeric_data$InactiveHighBalance))
numeric_data$SeniorSingleProduct <- as.numeric(as.character(numeric_data$SeniorSingleProduct))
numeric_data$FranceFemale <- as.numeric(as.character(numeric_data$FranceFemale))

# Convert target variable to factor with meaningful labels
numeric_data$Exited <- factor(numeric_data$Exited, levels = c("0", "1"), labels = c("No", "Yes"))

# Combine numerical data with dummy variables
model_data <- cbind(numeric_data, dummies)

# Scale numeric features (excluding binary features)
numeric_cols <- c("CreditScore", "Age", "Tenure", "Balance", "NumOfProducts", 
                  "EstimatedSalary", "BalanceToSalary", "ProductsPerYear")
model_data[numeric_cols] <- scale(model_data[numeric_cols])

# Split data into training (70%) and testing (30%) sets
set.seed(42)
train_index <- createDataPartition(model_data$Exited, p = 0.7, list = FALSE)
train_data <- model_data[train_index, ]
test_data <- model_data[-train_index, ]

# Check for NA values
na_counts_train <- colSums(is.na(train_data))
na_counts_test <- colSums(is.na(test_data))
print("NA values in training data:")

## [1] "NA values in training data:"

print(sum(na_counts_train))

## [1] 0

print("NA values in test data:")

## [1] "NA values in test data:"

print(sum(na_counts_test))

## [1] 0

# Handle any NA values if present
if(sum(na_counts_train) > 0 || sum(na_counts_test) > 0) {
  # Identify columns with NA values
  na_cols_train <- names(which(na_counts_train > 0))
  na_cols_test <- names(which(na_counts_test > 0))
  na_cols <- unique(c(na_cols_train, na_cols_test))
  
  # Remove these columns from both datasets
  if(length(na_cols) > 0) {
    train_data <- train_data[, !names(train_data) %in% na_cols]
    test_data <- test_data[, !names(test_data) %in% na_cols]
    print(paste("Removed columns with NA values:", paste(na_cols, collapse=", ")))
  }
}

# Verify class distribution is maintained
print("Class distribution in training set:")

## [1] "Class distribution in training set:"

table(train_data$Exited) %>% prop.table() %>% round(4) * 100

## 
##    No   Yes 
## 79.63 20.37

print("Class distribution in test set:")

## [1] "Class distribution in test set:"

table(test_data$Exited) %>% prop.table() %>% round(4) * 100

## 
##    No   Yes 
## 79.63 20.37

Our data preparation process creates a robust foundation for predictive modeling through:

1. Categorical Encoding: One-hot encoding of categorical variables to ensure proper model interpretation
2. Feature Scaling: Standardization of numerical features to prevent dominance of variables with large ranges
3. Data Partitioning: 70/30 train-test split while preserving the original class distribution
4. Data Quality Assurance: Verification of data integrity and handling of any NA values
5. Target Variable Preparation: Conversion to factor with meaningful labels for model compatibility

The final preprocessed dataset maintains the original churn ratio (~20.4%) in both training and test sets, ensuring reliable model evaluation.

3. Predictive Modeling

3.1 Model Training Setup

# Define cross-validation settings
train_control <- trainControl(
  method = "cv",
  number = 5,
  classProbs = TRUE,
  summaryFunction = twoClassSummary,
  verboseIter = FALSE
)

# Define evaluation metric
metric <- "ROC"

# Prepare formula (use dot notation to include all variables)
formula <- Exited ~ .

3.2 Logistic Regression

# Train logistic regression model
set.seed(123)
logit_model <- train(
  formula,
  data = train_data,
  method = "glm",
  family = "binomial",
  metric = metric,
  trControl = train_control
)

# Model summary
summary(logit_model$finalModel)

## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1898  -0.6283  -0.4216  -0.2513   2.9938  
## 
## Coefficients: (1 not defined because of singularities)
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -1.48314    0.27336  -5.426 5.77e-08 ***
## CreditScore                          -0.06095    0.03315  -1.838 0.066014 .  
## Age                                   0.62233    0.11278   5.518 3.43e-08 ***
## Tenure                               -0.16941    0.14264  -1.188 0.234974    
## Balance                               0.07955    0.09368   0.849 0.395775    
## NumOfProducts                        -0.00425    0.04575  -0.093 0.925991    
## HasCrCard                            -0.03694    0.07260  -0.509 0.610910    
## IsActiveMember                       -0.93598    0.09456  -9.898  < 2e-16 ***
## EstimatedSalary                       0.05065    0.03370   1.503 0.132787    
## BalanceToSalary                       0.06454    0.10654   0.606 0.544645    
## ZeroBalance                          -0.07057    0.18357  -0.384 0.700657    
## ProductsPerYear                      -0.05067    0.07567  -0.670 0.503106    
## InactiveHighBalance                   0.17321    0.11738   1.476 0.140045    
## SeniorSingleProduct                   0.59479    0.17386   3.421 0.000624 ***
## FranceFemale                          0.27874    0.13443   2.073 0.038134 *  
## GeographyFrance                      -0.19428    0.11159  -1.741 0.081669 .  
## GeographyGermany                      0.67232    0.09857   6.820 9.08e-12 ***
## GeographySpain                             NA         NA      NA       NA    
## GenderMale                           -0.39144    0.08902  -4.397 1.10e-05 ***
## `\\`AgeGroup30-40\\``                -0.03261    0.15143  -0.215 0.829513    
## `\\`AgeGroup40-50\\``                 0.80083    0.22254   3.599 0.000320 ***
## `\\`AgeGroup50-60\\``                 0.76031    0.35095   2.166 0.030278 *  
## `\\`AgeGroup60+\\``                  -0.97758    0.48895  -1.999 0.045570 *  
## `\\`TenureGroupDeveloping (2-5)\\``   0.07472    0.16146   0.463 0.643523    
## `\\`TenureGroupEstablished (5-8)\\``  0.13494    0.27126   0.497 0.618854    
## `\\`TenureGroupMature (8+)\\``        0.30610    0.37499   0.816 0.414345    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7077.6  on 7000  degrees of freedom
## Residual deviance: 5760.6  on 6976  degrees of freedom
## AIC: 5810.6
## 
## Number of Fisher Scoring iterations: 6

# Make predictions on test data
logit_pred <- predict(logit_model, newdata = test_data, type = "prob")[, "Yes"]

# Calculate ROC
logit_roc <- roc(test_data$Exited, logit_pred)
logit_auc <- auc(logit_roc)

# Plot ROC curve
plot(logit_roc, main = "ROC Curve for Logistic Regression", 
     col = "#3498db", lwd = 2)
text(0.7, 0.3, paste("AUC =", round(logit_auc, 3)), col = "#3498db")

# Feature importance
logit_importance <- varImp(logit_model)
plot(logit_importance, top = 20, main = "Logistic Regression - Feature Importance")

The logistic regression model achieves an AUC of 0.801, providing a solid baseline and valuable insights into key churn drivers:

1. Activity Status (100% importance): Highly significant negative association (-0.94, p < 0.001)
2. Geography: Germany (74%): Strong positive association (0.67, p < 0.001)
3. Age (56%): Significant positive relationship (0.62, p < 0.001)
4. Gender: Male (45%): Negative association (-0.39, p < 0.001)
5. Senior Single-Product (38%): Significant positive effect (0.59, p < 0.001)

This model highlights the critical role of customer engagement (activity status), geographic factors, and age-related variables in determining churn risk. The model’s performance provides a strong foundation for comparison with more complex approaches.

3.3 Random Forest

# Train Random Forest model
set.seed(123)
rf_params <- expand.grid(
  mtry = c(3, 5, 7)
)

rf_model <- train(
  formula,
  data = train_data,
  method = "rf",
  metric = metric,
  tuneGrid = rf_params,
  trControl = train_control,
  importance = TRUE
)

# Print best tuning parameters
print(rf_model$bestTune)

##   mtry
## 3    7

# Model performance
print(rf_model)

## Random Forest 
## 
## 7001 samples
##   25 predictor
##    2 classes: 'No', 'Yes' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 5600, 5601, 5601, 5601, 5601 
## Resampling results across tuning parameters:
## 
##   mtry  ROC        Sens       Spec     
##   3     0.8496266  0.9770404  0.3842964
##   5     0.8525927  0.9635874  0.4572273
##   7     0.8529607  0.9580269  0.4698540
## 
## ROC was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 7.

# Make predictions on test data
rf_pred <- predict(rf_model, newdata = test_data, type = "prob")[, "Yes"]

# Calculate ROC
rf_roc <- roc(test_data$Exited, rf_pred)
rf_auc <- auc(rf_roc)

# Plot ROC curve
plot(rf_roc, main = "ROC Curve for Random Forest", 
     col = "#2ecc71", lwd = 2)
text(0.7, 0.3, paste("AUC =", round(rf_auc, 3)), col = "#2ecc71")

# Feature importance
rf_importance <- varImp(rf_model)
plot(rf_importance, top = 20, main = "Random Forest - Feature Importance")

The Random Forest model achieves an AUC of 0.848, a substantial improvement over logistic regression. With optimal parameter mtry = 7, the model reveals a notably different feature importance pattern:

1. Number of Products (100% importance): Emerges as the dominant predictor, capturing the non-linear relationship observed in our EDA
2. Age (48%): Remains critical but with reduced relative importance
3. Activity Status (40%): Still influential but less dominant than in logistic regression
4. Balance-Related Factors (23-24%): Account balance and related metrics show moderate importance
5. Geography: Germany (23%): Regional market dynamics remain significant

The Random Forest model captures complex non-linear relationships and interaction effects more effectively than logistic regression, particularly highlighting the critical role of product portfolio management in customer retention strategies.

3.4 XGBoost

# Train XGBoost model
set.seed(123)
xgb_params <- expand.grid(
  nrounds = c(100, 200),
  max_depth = c(3, 6),
  eta = c(0.1, 0.3),
  gamma = 0,
  colsample_bytree = 0.7,
  min_child_weight = 1,
  subsample = 0.8
)

xgb_model <- train(
  formula,
  data = train_data,
  method = "xgbTree",
  metric = metric,
  tuneGrid = xgb_params,
  trControl = train_control,
  verbose = FALSE
)

## [20:16:09] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:09] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:10] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:10] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:10] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:10] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:11] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:11] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:12] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:12] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:13] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:13] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:13] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:13] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:15] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:15] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:15] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:15] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:16] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:16] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:17] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:17] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:18] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:18] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:18] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:18] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:19] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:19] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:20] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:20] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:21] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:21] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:22] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:22] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:23] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:23] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:23] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:23] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:24] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.
## [20:16:24] WARNING: src/c_api/c_api.cc:935: `ntree_limit` is deprecated, use `iteration_range` instead.

# Print best tuning parameters
print(xgb_model$bestTune)

##   nrounds max_depth eta gamma colsample_bytree min_child_weight subsample
## 1     100         3 0.1     0              0.7                1       0.8

# Model performance
print(xgb_model)

## eXtreme Gradient Boosting 
## 
## 7001 samples
##   25 predictor
##    2 classes: 'No', 'Yes' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 5600, 5601, 5601, 5601, 5601 
## Resampling results across tuning parameters:
## 
##   eta  max_depth  nrounds  ROC        Sens       Spec     
##   0.1  3          100      0.8619327  0.9639462  0.4635358
##   0.1  3          200      0.8611089  0.9608969  0.4810747
##   0.1  6          100      0.8574260  0.9565919  0.4782505
##   0.1  6          200      0.8484684  0.9510314  0.4922807
##   0.3  3          100      0.8533149  0.9540807  0.4908772
##   0.3  3          200      0.8456983  0.9470852  0.4993056
##   0.3  6          100      0.8367451  0.9434978  0.5021151
##   0.3  6          200      0.8301747  0.9395516  0.4979070
## 
## Tuning parameter 'gamma' was held constant at a value of 0
## Tuning
## 
## Tuning parameter 'min_child_weight' was held constant at a value of 1
## 
## Tuning parameter 'subsample' was held constant at a value of 0.8
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were nrounds = 100, max_depth = 3, eta
##  = 0.1, gamma = 0, colsample_bytree = 0.7, min_child_weight = 1 and subsample
##  = 0.8.

# Make predictions on test data
xgb_pred <- predict(xgb_model, newdata = test_data, type = "prob")[, "Yes"]

# Calculate ROC
xgb_roc <- roc(test_data$Exited, xgb_pred)
xgb_auc <- auc(xgb_roc)

# Plot ROC curve
plot(xgb_roc, main = "ROC Curve for XGBoost", 
     col = "#e74c3c", lwd = 2)
text(0.7, 0.3, paste("AUC =", round(xgb_auc, 3)), col = "#e74c3c")

# Feature importance
xgb_importance <- varImp(xgb_model)
plot(xgb_importance, top = 20, main = "XGBoost - Feature Importance")

The XGBoost model demonstrates the highest predictive performance with an AUC of 0.866, outperforming both logistic regression (0.801) and Random Forest (0.848). Through hyperparameter tuning, the optimal configuration was identified with shallower trees (depth 3), a conservative learning rate (0.1), and 100 boosting rounds.

Feature importance analysis reveals:

1. Age (100% importance): The single most influential predictor
2. Number of Products (84%): Nearly as important as age in this model
3. Activity Status (34%): Remains highly significant but with reduced relative importance
4. Balance (26%): Moderate importance consistent across tree-based models
5. Geography: Germany (17%): Still influential but less dominant than in logistic regression

The XGBoost model’s outstanding performance comes from its ability to capture both linear and non-linear relationships while minimizing overfitting through regularization. The preference for shallower trees indicates that simpler decision boundaries generalize better to unseen data in this domain.

3.5 Gradient Boosting Machine

# Train GBM model (as alternative to LightGBM)
set.seed(123)
gbm_params <- expand.grid(
  n.trees = c(100, 200),
  interaction.depth = c(3, 5),
  shrinkage = c(0.01, 0.1),
  n.minobsinnode = 10
)

gbm_model <- train(
  formula,
  data = train_data,
  method = "gbm",
  metric = metric,
  tuneGrid = gbm_params,
  trControl = train_control,
  verbose = FALSE
)

# Print best tuning parameters
print(gbm_model$bestTune)

##   n.trees interaction.depth shrinkage n.minobsinnode
## 5     100                 3       0.1             10

# Model performance
print(gbm_model)

## Stochastic Gradient Boosting 
## 
## 7001 samples
##   25 predictor
##    2 classes: 'No', 'Yes' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 5600, 5601, 5601, 5601, 5601 
## Resampling results across tuning parameters:
## 
##   shrinkage  interaction.depth  n.trees  ROC        Sens       Spec     
##   0.01       3                  100      0.8286604  0.9967713  0.1409471
##   0.01       3                  200      0.8438963  0.9757848  0.3850031
##   0.01       5                  100      0.8379889  0.9890583  0.3078444
##   0.01       5                  200      0.8532030  0.9770404  0.3913238
##   0.10       3                  100      0.8624666  0.9610762  0.4733530
##   0.10       3                  200      0.8609870  0.9587444  0.4915912
##   0.10       5                  100      0.8604469  0.9594619  0.4972028
##   0.10       5                  200      0.8559009  0.9565919  0.4964863
## 
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were n.trees = 100, interaction.depth =
##  3, shrinkage = 0.1 and n.minobsinnode = 10.

# Make predictions on test data
gbm_pred <- predict(gbm_model, newdata = test_data, type = "prob")[, "Yes"]

# Calculate ROC
gbm_roc <- roc(test_data$Exited, gbm_pred)
gbm_auc <- auc(gbm_roc)

# Plot ROC curve
plot(gbm_roc, main = "ROC Curve for GBM", 
     col = "#9b59b6", lwd = 2)
text(0.7, 0.3, paste("AUC =", round(gbm_auc, 3)), col = "#9b59b6")

# Feature importance
gbm_importance <- varImp(gbm_model)
plot(gbm_importance, top = 20, main = "GBM - Feature Importance")

The Gradient Boosting Machine (GBM) model delivers exceptional predictive performance with an AUC of 0.864, nearly matching XGBoost’s 0.866. The optimal configuration includes 100 trees with depth 3 and a learning rate of 0.1, mirroring XGBoost’s parameters.

Feature importance is highly consistent with XGBoost:

1. Age (100% importance): Confirmed as the dominant predictor
2. Number of Products (87%): Nearly equal in importance to age
3. Activity Status (33%): Moderate but substantial contribution
4. Balance (21%) and Geography: Germany (20%): Secondary but meaningful predictors

The remarkable similarity between GBM and XGBoost results confirms the robustness of our findings and suggests that we may be approaching the maximum predictive power available from this feature set. The preference for similar hyperparameters (100 trees, depth 3) across both boosting implementations indicates this represents an appropriate model complexity for the underlying data structure.

3.6 Model Comparison

# Combine all ROC curves
plot(logit_roc, col = "#3498db", lwd = 2, main = "Model Comparison - ROC Curves")
plot(rf_roc, col = "#2ecc71", lwd = 2, add = TRUE)
plot(xgb_roc, col = "#e74c3c", lwd = 2, add = TRUE)
plot(gbm_roc, col = "#9b59b6", lwd = 2, add = TRUE)
legend("bottomright", 
       legend = c(paste("Logistic Regression (AUC:", round(logit_auc, 3), ")"),
                  paste("Random Forest (AUC:", round(rf_auc, 3), ")"),
                  paste("XGBoost (AUC:", round(xgb_auc, 3), ")"),
                  paste("GBM (AUC:", round(gbm_auc, 3), ")")),
       col = c("#3498db", "#2ecc71", "#e74c3c", "#9b59b6"),
       lwd = 2)

# Create a table of results
model_results <- data.frame(
  Model = c("Logistic Regression", "Random Forest", "XGBoost", "GBM"),
  AUC = c(logit_auc, rf_auc, xgb_auc, gbm_auc)
)

# Sort by performance
model_results <- model_results[order(-model_results$AUC), ]
kable(model_results, caption = "Model Performance Comparison", digits = 4)

Model Performance Comparison

	Model	AUC
3	XGBoost	0.8661
4	GBM	0.8640
2	Random Forest	0.8478
1	Logistic Regression	0.8011

Our comprehensive modeling approach included four different algorithms, each offering unique insights into customer churn prediction:

Model	AUC	Key Strengths
XGBoost	0.866	Best overall performance, balanced feature importance
GBM	0.864	Nearly identical to XGBoost, confirms key patterns
Random Forest	0.848	Strong performance, emphasizes product relationships
Logistic Regression	0.801	Baseline model, high interpretability

The ROC curves visualization clearly demonstrates the performance hierarchy, with both gradient boosting implementations (XGBoost and GBM) achieving superior discrimination ability. The convergence of these curves at both extremes indicates that all models perform similarly for the most obvious cases, but differ in their ability to correctly classify borderline customers.

Based on these results, we select XGBoost as our best model for further evaluation and deployment due to its superior performance and well-balanced feature importance profile.

4. Best Model Evaluation

# Determine best model based on AUC
best_model_name <- model_results$Model[1]
if (best_model_name == "Logistic Regression") {
  best_model <- logit_model
  best_pred <- logit_pred
} else if (best_model_name == "Random Forest") {
  best_model <- rf_model
  best_pred <- rf_pred
} else if (best_model_name == "XGBoost") {
  best_model <- xgb_model
  best_pred <- xgb_pred
} else if (best_model_name == "GBM") {
  best_model <- gbm_model
  best_pred <- gbm_pred
}

# Get predicted probabilities from best model
best_pred_probs <- best_pred

# Convert to binary predictions using different thresholds
evaluate_threshold <- function(threshold) {
  predicted_classes <- ifelse(best_pred_probs >= threshold, "Yes", "No")
  conf_matrix <- confusionMatrix(
    factor(predicted_classes, levels = c("No", "Yes")),
    test_data$Exited
  )
  
  precision <- conf_matrix$byClass["Pos Pred Value"]
  recall <- conf_matrix$byClass["Sensitivity"]
  f1 <- 2 * precision * recall / (precision + recall)
  
  return(data.frame(
    Threshold = threshold,
    Accuracy = conf_matrix$overall["Accuracy"],
    Precision = precision,
    Recall = recall,
    F1_Score = f1,
    TP = conf_matrix$table[2, 2],
    FP = conf_matrix$table[2, 1],
    TN = conf_matrix$table[1, 1],
    FN = conf_matrix$table[1, 2]
  ))
}

# Try different thresholds
thresholds <- seq(0.1, 0.9, by = 0.1)
threshold_results <- do.call(rbind, lapply(thresholds, evaluate_threshold))

# Display results
kable(threshold_results, caption = "Performance Metrics at Different Threshold Values", digits = 4)

Performance Metrics at Different Threshold Values

	Threshold	Accuracy	Precision	Recall	F1_Score	TP	FP	TN	FN
Accuracy	0.1	0.6532	0.9610	0.5884	0.7299	554	983	1405	57
Accuracy1	0.2	0.7953	0.9277	0.8057	0.8624	461	464	1924	150
Accuracy2	0.3	0.8409	0.9030	0.8966	0.8998	381	247	2141	230
Accuracy3	0.4	0.8530	0.8874	0.9338	0.9100	328	158	2230	283
Accuracy4	0.5	0.8633	0.8769	0.9636	0.9182	288	87	2301	323
Accuracy5	0.6	0.8633	0.8636	0.9837	0.9197	240	39	2349	371
Accuracy6	0.7	0.8570	0.8527	0.9916	0.9169	202	20	2368	409
Accuracy7	0.8	0.8433	0.8367	0.9979	0.9102	146	5	2383	465
Accuracy8	0.9	0.8216	0.8172	0.9996	0.8992	77	1	2387	534

# Plot threshold impact on metrics
threshold_results_long <- threshold_results %>%
  select(Threshold, Accuracy, Precision, Recall, F1_Score) %>%
  pivot_longer(cols = c("Accuracy", "Precision", "Recall", "F1_Score"),
               names_to = "Metric", values_to = "Value")

ggplot(threshold_results_long, aes(x = Threshold, y = Value, color = Metric, group = Metric)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_color_viridis_d(option = "D") +
  labs(title = "Impact of Threshold on Performance Metrics",
       x = "Threshold",
       y = "Metric Value") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Select optimal threshold based on F1 score
optimal_threshold <- threshold_results %>%
  filter(F1_Score == max(F1_Score)) %>%
  pull(Threshold)

# Create final predictions with optimal threshold
final_predictions <- ifelse(best_pred_probs >= optimal_threshold, "Yes", "No")
final_predictions <- factor(final_predictions, levels = c("No", "Yes"))

# Create final confusion matrix
conf_matrix <- confusionMatrix(final_predictions, test_data$Exited)
print(conf_matrix)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   No  Yes
##        No  2349  371
##        Yes   39  240
##                                           
##                Accuracy : 0.8633          
##                  95% CI : (0.8505, 0.8754)
##     No Information Rate : 0.7963          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.4719          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9837          
##             Specificity : 0.3928          
##          Pos Pred Value : 0.8636          
##          Neg Pred Value : 0.8602          
##              Prevalence : 0.7963          
##          Detection Rate : 0.7833          
##    Detection Prevalence : 0.9070          
##       Balanced Accuracy : 0.6882          
##                                           
##        'Positive' Class : No              
## 

# Visualize confusion matrix
conf_data <- as.data.frame(conf_matrix$table)
colnames(conf_data) <- c("Predicted", "Actual", "Freq")

ggplot(conf_data, aes(x = Predicted, y = Actual, fill = Freq)) +
  geom_tile() +
  geom_text(aes(label = Freq), color = "white", size = 10) +
  scale_fill_gradient(low = "#3498db", high = "#e74c3c") +
  labs(title = paste("Confusion Matrix (Threshold =", optimal_threshold, ")"),
       x = "Predicted Class",
       y = "Actual Class") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

The prediction threshold fundamentally impacts the practical application of our model. Our threshold analysis reveals critical trade-offs:

• Low thresholds (0.1-0.3): Maximize recall (identifying nearly all churning customers) but generate many false positives, reducing precision
• High thresholds (0.7-0.9): Maximize precision (high confidence in predicted churners) but miss many actual churners
• Optimal threshold (0.6): Balances trade-offs with highest F1-score (0.92)

The confusion matrix at our optimal threshold shows:

• True Negatives (2,349): Correctly identified retained customers

• True Positives (240): Correctly identified churning customers

• False Negatives (371): Churners incorrectly classified as staying

• False Positives (39): Stayers incorrectly classified as churning

This configuration provides an optimal balance between identifying customers at risk of churning and minimizing false alarms that would lead to unnecessary retention interventions.

4.1 Feature Importance Analysis

# Extract feature importance from best model
best_importance <- varImp(best_model)

# Plot feature importance
plot(best_importance, top = 15, 
     main = paste(best_model_name, "- Top 15 Features"))

# Get importance values as data frame
imp_df <- best_importance$importance
imp_df <- data.frame(
  Feature = rownames(imp_df),
  Importance = imp_df$Overall
)

# Sort by importance
imp_df <- imp_df[order(-imp_df$Importance), ]

# Display top 15 features
kable(head(imp_df, 15), caption = "Top 15 Most Important Features", digits = 2)

Top 15 Most Important Features

Feature	Importance
Age	100.00
NumOfProducts	84.98
IsActiveMember	27.99
Balance	19.54
GeographyGermany	13.76
BalanceToSalary	9.78
CreditScore	8.89
EstimatedSalary	7.41
GenderMale	6.22
InactiveHighBalance	4.31
AgeGroup50-60	2.52
GeographyFrance	2.07
ProductsPerYear	2.03
ZeroBalance	1.99
Tenure	1.54

Across all models, a clear consensus emerges on the key drivers of customer churn:

1. Age: Consistently the top predictor, with dramatically higher churn rates among older customers
2. Number of Products: Strong predictor with a non-linear relationship showing highest churn at extremes (1 or 3-4 products)
3. Activity Status: Critical engagement indicator, with inactive customers showing substantially higher churn risk
4. Balance: Moderate predictor with complex relationships to churn behavior
5. Geography: Germany: Regional market dynamics significantly impact retention

These five features together capture the essential dimensions of churn risk: customer demographics, product relationships, engagement level, financial position, and market-specific factors. The robust performance across different modeling approaches and the consistent feature importance rankings provide high confidence in these findings.

5. Customer Risk Profiling

# Add prediction probabilities back to the test data
test_profiles <- test_data
test_profiles$PredictedProbability <- best_pred_probs
test_profiles$PredictedChurn <- final_predictions

# Create risk groups based on predicted probabilities
test_profiles$RiskLevel <- cut(
  test_profiles$PredictedProbability,
  breaks = c(0, 0.25, 0.5, 0.75, 1),
  labels = c("Low Risk", "Medium Risk", "High Risk", "Very High Risk"),
  include.lowest = TRUE
)

# Distribution of risk levels
ggplot(test_profiles, aes(x = RiskLevel, fill = RiskLevel)) +
  geom_bar() +
  scale_fill_viridis_d(option = "D") +
  labs(title = "Distribution of Customer Risk Levels",
       x = "Risk Level",
       y = "Count") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Analyze characteristics of high risk customers
high_risk_profile <- test_profiles %>%
  filter(RiskLevel %in% c("High Risk", "Very High Risk")) %>%
  select(Age, Balance, NumOfProducts, IsActiveMember, GeographyFrance, GeographyGermany, GenderMale, Exited) %>%
  mutate(
    Age = Age * sd(churn_clean$Age) + mean(churn_clean$Age),  # Denormalize
    Balance = Balance * sd(churn_clean$Balance) + mean(churn_clean$Balance),  # Denormalize
    IsActiveMember = factor(IsActiveMember, levels = c(0, 1), labels = c("Inactive", "Active")),
    Geography = case_when(
      GeographyFrance == 1 ~ "France",
      GeographyGermany == 1 ~ "Germany",
      TRUE ~ "Spain"
    ),
    Gender = ifelse(GenderMale == 1, "Male", "Female")
  ) %>%
  select(-GeographyFrance, -GeographyGermany, -GenderMale)

# Summarize high risk segment
high_risk_summary <- high_risk_profile %>%
  group_by(Exited) %>%
  summarise(
    Count = n(),
    AvgAge = mean(Age),
    AvgBalance = mean(Balance),
    PercentInactive = mean(IsActiveMember == "Inactive") * 100,
    PercentFrance = mean(Geography == "France") * 100,
    PercentFemale = mean(Gender == "Female") * 100,
    AvgProducts = mean(NumOfProducts)
  )

# Display high risk profile summary
kable(high_risk_summary, caption = "High Risk Customer Profile", digits = 1)

High Risk Customer Profile

Exited	Count	AvgAge	AvgBalance	PercentInactive	PercentFrance	PercentFemale	AvgProducts
No	87	47.7	90890.3	81.6	37.9	52.9	-0.5
Yes	288	49.9	94208.7	76.7	37.2	60.1	0.3

# Create key risk segments
risk_segments <- test_profiles %>%
  mutate(
    Age = Age * sd(churn_clean$Age) + mean(churn_clean$Age),  # Denormalize
    Balance = Balance * sd(churn_clean$Balance) + mean(churn_clean$Balance),  # Denormalize
    IsActiveMember = factor(IsActiveMember, levels = c(0, 1), labels = c("Inactive", "Active")),
    Geography = case_when(
      GeographyFrance == 1 ~ "France",
      GeographyGermany == 1 ~ "Germany",
      TRUE ~ "Spain"
    ),
    Gender = ifelse(GenderMale == 1, "Male", "Female"),
    SeniorCustomer = Age >= 50,
    HighBalance = Balance >= 100000,
    SingleProduct = NumOfProducts == 1
  )

# Define key segments
segments <- risk_segments %>%
  mutate(
    Segment = case_when(
      SeniorCustomer & Geography == "France" & IsActiveMember == "Inactive" ~ "Senior Inactive French",
      SeniorCustomer & SingleProduct ~ "Senior Single-Product",
      HighBalance & IsActiveMember == "Inactive" ~ "Inactive High-Balance",
      Geography == "France" & Gender == "Female" ~ "French Female",
      TRUE ~ "Other"
    )
  )

# Analyze churn rates by segment
segment_analysis <- segments %>%
  group_by(Segment) %>%
  summarise(
    Count = n(),
    ChurnRate = mean(Exited == "Yes") * 100,
    AvgRiskScore = mean(PredictedProbability) * 100
  ) %>%
  filter(Count >= 30) %>%  # Filter for significant segments
  arrange(desc(ChurnRate))

# Display segment analysis
kable(segment_analysis, caption = "Key Customer Segments and Churn Risk", digits = 1)

Key Customer Segments and Churn Risk

Segment	Count	ChurnRate	AvgRiskScore
Senior Inactive French	63	81.0	74.5
Inactive High-Balance	686	30.3	31.1
Other	1731	15.7	15.1
French Female	519	15.4	16.3

# Visualize segment churn rates
ggplot(segment_analysis, aes(x = reorder(Segment, ChurnRate), y = ChurnRate, fill = ChurnRate)) +
  geom_bar(stat = "identity") +
  scale_fill_gradient(low = "#3498db", high = "#e74c3c") +
  geom_text(aes(label = paste0(round(ChurnRate, 1), "%")), hjust = -0.1) +
  coord_flip() +
  labs(title = "Churn Rate by Customer Segment",
       x = "Segment",
       y = "Churn Rate (%)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

Our XGBoost model enables a granular risk assessment of the customer base, categorizing customers into four risk tiers:

• Low Risk (0-25%): 73% of customers with minimal churn concern

• Medium Risk (25-50%): 13% of customers with moderate churn risk

• High Risk (50-75%): 7% of customers requiring proactive attention

• Very High Risk (75-100%): 7% of customers in critical need of retention intervention

This risk stratification creates a clear prioritization framework for retention efforts, allowing the bank to focus resources on the ~14% of customers with high or very high churn probability.

Examining high-risk customers reveals distinct profiles with the following characteristics:

• Average age of ~50 years (significantly higher than bank average)

• Substantial average balance of ~€90,000-94,000

• 75-80% are inactive members

• Female-skewed (53-60%)

• Higher geographic concentration in France

Our analysis identifies four high-value customer segments with distinct risk profiles:

1. Senior Inactive French Customers (81.0% churn rate):
   • Extremely high-risk segment with 81% churn
   • Combination of age, inactivity, and French market creates critical retention challenge
2. Inactive High-Balance Customers (30.3% churn rate):
   • Substantial segment with significant financial value at risk
   • Nearly one-third leave despite their high value to the bank
3. French Female Customers (15.4% churn rate):
   • Moderate risk segment
   • May benefit from targeted female-focused retention strategies
4. Other Customers (15.7% churn rate):
   • Largest segment with varied characteristics
   • Requires further subsegmentation for effective targeting

These actionable segments provide a clear framework for developing targeted retention strategies, with prioritization based on both churn probability and segment value.

6. Business Recommendations

6.1 Key Findings and Implications

Our comprehensive analysis and predictive modeling have revealed several critical insights with significant business implications:

1. Age is the Dominant Predictor: Age consistently emerged as the most important predictor across all models, with customers over 50 showing churn rates up to 55% (for the 50-60 age group). This strong age effect suggests targeted age-specific offerings are essential for retention.

2. Product Portfolio is Critical: Number of products showed unexpected importance, with a clear U-shaped relationship to churn. Specifically, customers with 2 products demonstrate optimal retention (8% churn), while those with 3-4 products show extremely high churn rates (80-100%), indicating potential product overloading.

3. Activity Status Drives Engagement: Inactive members are more than twice as likely to churn compared to active customers. This engagement metric represents a high-leverage opportunity for intervention through reactivation strategies.

4. Geographic Variation is Significant: German customers exhibit substantially higher churn rates (~33%) compared to France and Spain (both ~15%), suggesting targeted market-specific approaches are needed.

5. High-Value Customer Risk: Inactive customers with high balances represent a particularly vulnerable segment with 30.3% churn rate, putting significant deposit value at risk.

6. Segment Combinations Matter: When key risk factors combine (older age, inactive status, specific geography), churn rates approach 81% for segments like “Senior Inactive French Customers.”

7. Model Performance: Our XGBoost model achieved an AUC of 0.866, outperforming all other approaches and providing reliable risk stratification capability.

6.2 Recommended Retention Strategies

6.2.1 Segment-Specific Interventions

Segment	Risk Profile	Recommended Actions
Senior Inactive French Customers	Very High Risk (81% churn rate)	1. Dedicated relationship managers for personalized service 2. Special rate offers to incentivize activity 3. Specialized product offerings for retirement planning 4. French-specific customer experience improvements
Inactive High-Balance Customers	High Risk (30.3% churn rate)	1. Proactive wealth management consultations 2. Premium account services activation 3. Investment portfolio reviews 4. Competitive rate alerts for comparable products
Customers with 3-4 Products	Very High Risk (80-100% churn rate)	1. Product portfolio simplification consultations 2. Fee consolidation or reduction programs 3. Service quality review and enhancement 4. Dedicated support for complex banking relationships
German Customers	High Risk (33% churn rate)	1. Market-specific competitive analysis and response 2. Enhanced product offerings tailored to German market 3. Service quality improvements at German branches 4. Customer satisfaction monitoring and response system

6.2.2 General Retention Framework

Data-Driven Risk Monitoring:

• Implement real-time churn risk scoring using our optimized XGBoost model (AUC 0.866)

• Configure the prediction threshold at 0.6 for optimal balance of precision and recall

• Create automated alert systems when customers cross risk thresholds

• Develop dashboards showing customer risk levels and key contributing factors

Product Strategy Optimization:

• Review and potentially simplify offerings for customers with 3-4 products, as this represents a critical risk factor

• Focus cross-selling efforts on reaching exactly 2 products per customer (optimal retention point)

• Develop age-appropriate product bundles with compelling value propositions for different life stages

• Create special retention offers for customers identified in the highest risk segments

Engagement Enhancement:

• Prioritize reactivation campaigns for inactive customers, particularly those with high balances

• Develop digital engagement strategies that appeal to different age segments

• Create specific engagement metrics and goals for relationship managers

• Implement a structured customer contact program based on risk level and segment

Personalized Customer Experience:

• Develop tailored service models for identified high-risk segments

• Train customer-facing staff on segment-specific retention strategies

• Create specialized retention teams focused on the highest-risk customer groups

• Implement feedback mechanisms to continuously refine segment-based approaches

6.3 Implementation Roadmap

Immediate Actions (1-3 Months):

• Deploy the XGBoost prediction model as a scoring system integrated with CRM

• Configure real-time risk alerts for customers crossing the 0.6 probability threshold

• Launch targeted retention campaigns for the “Senior Inactive French” segment

• Begin portfolio review for customers with 3-4 products

Medium-Term Initiatives (3-6 Months):

• Implement comprehensive retention strategies for all identified high-risk segments

• Develop specialized product offerings for key age groups and risk segments

• Roll out training programs for relationship managers on using the prediction model

• Establish regular review process for segment performance and strategy effectiveness

Long-Term Strategy (6-12 Months):

• Integrate churn prediction into all customer touchpoints and decision systems

• Develop a comprehensive customer lifetime value optimization framework

• Create an automated testing program for retention strategies

• Implement continuous model retraining and refinement based on new churn data

6.4 Expected Business Impact

Quantitative Projections:

• Reduce overall churn rate from 20.4% to 14-16% within 12 months

• Increase retention of high-value customers by 30% from current levels

• Improve product balance optimization with 25% reduction in customers with 3-4 products

• Achieve 15% improvement in inactive customer reactivation

Financial Impact:

• Estimated revenue preservation of €3-5 million from reduced churn

• Balance retention improvement of €15-20 million from high-balance customer strategies

• 10-15% improvement in customer lifetime value

• 25-30% increase in marketing ROI through targeted retention versus acquisition

Strategic Benefits:

• Enhanced competitive position through data-driven customer management

• Improved product development informed by retention risk factors

• Creation of a continuous learning cycle for customer behavior understanding

• Development of predictive capabilities applicable to other business challenges

7. Conclusion

This analysis has demonstrated the substantial value of advanced predictive modeling for bank customer retention. By applying gradient boosting techniques to comprehensive customer data, we’ve developed a highly accurate model (AUC 0.866) capable of identifying customers at risk of churning before they leave.

Our modeling approach revealed that while traditional demographic factors like age remain critical predictors, behavioral factors like product portfolio composition and account activity are equally influential in determining churn risk. The identification of specific high-risk customer segments—particularly “Senior Inactive French Customers” with an 81% churn rate and “Inactive High-Balance Customers” with a 30.3% churn rate—provides clear priorities for intervention.

The threshold optimization analysis demonstrates that with a 0.6 probability threshold, the bank can achieve 86.3% accuracy in churn prediction while maintaining a favorable balance between precision (86.4%) and recall (98.4%). This configuration creates a practical implementation framework for deploying the model in production environments.

By implementing the segment-specific retention strategies outlined in this analysis, the bank has the opportunity to significantly reduce customer attrition while optimizing marketing expenditure. The potential 25-30% reduction in churn would translate to substantial preserved revenue and enhanced customer lifetime value.

In the evolving competitive landscape of retail banking, this predictive approach to customer retention represents a significant advancement over traditional reactive methods. The model not only identifies who is likely to leave but provides insights into why they might leave, enabling truly proactive relationship management.