Credit Card Users Churn Prediction

Problem Statement

Business Context

The Thera Bank recently saw a steep decline in the number of users of their credit card, credit cards are a good source of income for the bank because of different kinds of fees can be charged by the bank such as annual fees, balance transfer fees, cash advance fees, late payment fees, foreign transaction fees, and others. Some fees are charged to every user irrespective of usage, while others are charged under specified circumstances.

Customers’ leaving credit cards services would lead bank to loss, so Thera Bank wants to analyze the data of customers and identify who will leave their credit card services in advance and reason for that – so that bank could improve upon those areas.

You, as a Data scientist at Thera bank, need to come up with a classification model that will help the bank improve its services so that customers do not renounce their credit cards.

💡 This is a binary classification problem.

Data Description

• CLIENTNUM: Client number. Unique identifier for the customer holding the account
• Attrition_Flag: Internal event (customer activity) variable - if the account is closed then "Attrited Customer" else "Existing Customer"
• Customer_Age: Age in Years
• Gender: Gender of the account holder
• Dependent_count: Number of dependents
• Education_Level: Educational Qualification of the account holder - Graduate, High School, Unknown, Uneducated, College(refers to college student), Post-Graduate, Doctorate
• Marital_Status: Marital Status of the account holder
• Income_Category: Annual Income Category of the account holder
• Card_Category: Type of Card
• Months_on_book: Period of relationship with the bank (in months)
• Total_Relationship_Count: Total no. of products held by the customer
• Months_Inactive_12_mon: No. of months inactive in the last 12 months
• Contacts_Count_12_mon: No. of Contacts in the last 12 months
• Credit_Limit: Credit Limit on the Credit Card
• Total_Revolving_Bal: Total Revolving Balance on the Credit Card
• Avg_Open_To_Buy: Open to Buy Credit Line (Average of last 12 months)
• Total_Amt_Chng_Q4_Q1: Change in Transaction Amount (Q4 over Q1)
• Total_Trans_Amt: Total Transaction Amount (Last 12 months)
• Total_Trans_Ct: Total Transaction Count (Last 12 months)
• Total_Ct_Chng_Q4_Q1: Change in Transaction Count (Q4 over Q1)
• Avg_Utilization_Ratio: Average Card Utilization Ratio

What Is a Revolving Balance?

• If a customer doesn't pay the balance of the revolving credit account in full every month, the unpaid portion carries over to the next month. That's called a Revolving Balance.

What is the Average Open to buy?

• Open to Buy means the amount left on your credit card to use. Now, this column represents the average of this value for the last 12 months.

What is the Average Utilization Ratio?

• The Avg_Utilization_Ratio represents how much of the available credit the customer spent. This is useful for calculating credit scores.

Relation b/w Avg_Open_To_Buy, Credit_Limit and Avg_Utilization_Ratio:

• ( Avg_Open_To_Buy / Credit_Limit ) + Avg_Utilization_Ratio = 1

Importing necessary libraries

# For Google Colab
# !pip install mlxtend==0.22.0 scikit-learn==1.2.2 imbalanced-learn==0.11.0 seaborn==0.13.1 matplotlib==3.8.0 numpy==1.25.2 pandas==2.2.3 xgboost==2.0.3 -q --user

   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/1.4 MB ? eta -:--:--
   ━━╺━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.1/1.4 MB 2.4 MB/s eta 0:00:01
   ━━━━╺━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.1/1.4 MB 2.4 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━━━━━━━━━━ 0.5/1.4 MB 4.7 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸━━━━━ 1.2/1.4 MB 8.1 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.4/1.4 MB 7.9 MB/s eta 0:00:00

Note: After running the above cell once, restart the notebook kernel and continue below.

# Libraries to help with reading and manipulating data
import numpy as np
import pandas as pd

# Libraries to help with data visualization
import matplotlib.pyplot as plt
import seaborn as sns

# Scaling and One-hot encoding
from sklearn.preprocessing import StandardScaler, MinMaxScaler, OneHotEncoder

# SimpleImputer for missing values
from sklearn.impute import SimpleImputer
# k-Nearest Neighbors imputer for advanced missing value treatment
from sklearn.impute import KNNImputer

# Library to split data and tune hyperparameters
from sklearn.model_selection import train_test_split, RandomizedSearchCV

# Libraries for building classifier models
from sklearn.tree import DecisionTreeClassifier
# Logistic Regression (simplest classification model) used as a baseline for comparison when evaluating the performance
from sklearn.linear_model import LogisticRegression

# Libraries for building ensemble classifier models
from sklearn.ensemble import BaggingClassifier, RandomForestClassifier, AdaBoostClassifier, GradientBoostingClassifier
from xgboost import XGBClassifier  # This is a separate library, it is not sklearn

# Libraries for model evaluation
from sklearn import metrics
from sklearn.metrics import confusion_matrix, classification_report, roc_auc_score
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score

# SMOTE for oversampling & RandomUnderSampler for undersampling
from imblearn.over_sampling import SMOTE
from imblearn.under_sampling import RandomUnderSampler

# Supress scientific notations
pd.set_option("display.float_format", lambda x: "%.3f" % x)

# Suppress warnings
import warnings
warnings.filterwarnings('ignore')

Loading the Dataset

# Mounting Google Drive to current Colab Notebook
from google.colab import drive
drive.mount('/content/drive')

# Loading the original dataset
originaldataset = pd.read_csv("/content/drive/My Drive/Texas McCombs/Machine Learning/Project #3/churners-dataset.csv")

# Take a copy of the dataset
data = originaldataset.copy()

Mounted at /content/drive

Data Inspection and Cleaning

• Initial data inspection
• Sanity checks

# Have a look at the first 5 rows of the dataset
data.head()

	CLIENTNUM	Attrition_Flag	Customer_Age	Gender	Dependent_count	Education_Level	Marital_Status	Income_Category	Card_Category	Months_on_book	...	Months_Inactive_12_mon	Contacts_Count_12_mon	Credit_Limit	Total_Revolving_Bal	Avg_Open_To_Buy	Total_Amt_Chng_Q4_Q1	Total_Trans_Amt	Total_Trans_Ct	Total_Ct_Chng_Q4_Q1	Avg_Utilization_Ratio
0	768805383	Existing Customer	45	M	3	High School	Married	$60K - $80K	Blue	39	...	1	3	12691.000	777	11914.000	1.335	1144	42	1.625	0.061
1	818770008	Existing Customer	49	F	5	Graduate	Single	Less than $40K	Blue	44	...	1	2	8256.000	864	7392.000	1.541	1291	33	3.714	0.105
2	713982108	Existing Customer	51	M	3	Graduate	Married	$80K - $120K	Blue	36	...	1	0	3418.000	0	3418.000	2.594	1887	20	2.333	0.000
3	769911858	Existing Customer	40	F	4	High School	NaN	Less than $40K	Blue	34	...	4	1	3313.000	2517	796.000	1.405	1171	20	2.333	0.760
4	709106358	Existing Customer	40	M	3	Uneducated	Married	$60K - $80K	Blue	21	...	1	0	4716.000	0	4716.000	2.175	816	28	2.500	0.000

5 rows × 21 columns

# Have a look at the last 5 rows
data.tail()

	CLIENTNUM	Attrition_Flag	Customer_Age	Gender	Dependent_count	Education_Level	Marital_Status	Income_Category	Card_Category	Months_on_book	...	Months_Inactive_12_mon	Contacts_Count_12_mon	Credit_Limit	Total_Revolving_Bal	Avg_Open_To_Buy	Total_Amt_Chng_Q4_Q1	Total_Trans_Amt	Total_Trans_Ct	Total_Ct_Chng_Q4_Q1	Avg_Utilization_Ratio
10122	772366833	Existing Customer	50	M	2	Graduate	Single	$40K - $60K	Blue	40	...	2	3	4003.000	1851	2152.000	0.703	15476	117	0.857	0.462
10123	710638233	Attrited Customer	41	M	2	NaN	Divorced	$40K - $60K	Blue	25	...	2	3	4277.000	2186	2091.000	0.804	8764	69	0.683	0.511
10124	716506083	Attrited Customer	44	F	1	High School	Married	Less than $40K	Blue	36	...	3	4	5409.000	0	5409.000	0.819	10291	60	0.818	0.000
10125	717406983	Attrited Customer	30	M	2	Graduate	NaN	$40K - $60K	Blue	36	...	3	3	5281.000	0	5281.000	0.535	8395	62	0.722	0.000
10126	714337233	Attrited Customer	43	F	2	Graduate	Married	Less than $40K	Silver	25	...	2	4	10388.000	1961	8427.000	0.703	10294	61	0.649	0.189

5 rows × 21 columns

# Understand the shape of the dataset
data.shape

(10127, 21)

# Check data types and number of non-null values for each column
data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 10127 entries, 0 to 10126
Data columns (total 21 columns):
 #   Column                    Non-Null Count  Dtype  
---  ------                    --------------  -----  
 0   CLIENTNUM                 10127 non-null  int64  
 1   Attrition_Flag            10127 non-null  object 
 2   Customer_Age              10127 non-null  int64  
 3   Gender                    10127 non-null  object 
 4   Dependent_count           10127 non-null  int64  
 5   Education_Level           8608 non-null   object 
 6   Marital_Status            9378 non-null   object 
 7   Income_Category           10127 non-null  object 
 8   Card_Category             10127 non-null  object 
 9   Months_on_book            10127 non-null  int64  
 10  Total_Relationship_Count  10127 non-null  int64  
 11  Months_Inactive_12_mon    10127 non-null  int64  
 12  Contacts_Count_12_mon     10127 non-null  int64  
 13  Credit_Limit              10127 non-null  float64
 14  Total_Revolving_Bal       10127 non-null  int64  
 15  Avg_Open_To_Buy           10127 non-null  float64
 16  Total_Amt_Chng_Q4_Q1      10127 non-null  float64
 17  Total_Trans_Amt           10127 non-null  int64  
 18  Total_Trans_Ct            10127 non-null  int64  
 19  Total_Ct_Chng_Q4_Q1       10127 non-null  float64
 20  Avg_Utilization_Ratio     10127 non-null  float64
dtypes: float64(5), int64(10), object(6)
memory usage: 1.6+ MB

• Notice that there are total of 21 columns and 10,127 rows in the dataset.
• Data type is either integer, float or the object type. Object types cannot be used for model building, we will convert them to categories.
• Education_Level and Marital_Status columns have null values. We will further check this by using isna() method.

💡 Converting data types from "object" to "category" reduces the memory usage of the DataFrame.

# Check for duplicate values in the data
data.duplicated().sum()

0

• There is no duplicate rows in the data.
  👉 It is always a good idea to check for duplicate values after removing the index column (e.g., CLIENTNUM). If the index column is unique, the duplicate-checking will always return 0 prior to its removal, as the uniqueness of the index ensures no duplicates.

# Check the number of unique values in each column
data.nunique()

	0
CLIENTNUM	10127
Attrition_Flag	2
Customer_Age	45
Gender	2
Dependent_count	6
Education_Level	6
Marital_Status	3
Income_Category	6
Card_Category	4
Months_on_book	44
Total_Relationship_Count	6
Months_Inactive_12_mon	7
Contacts_Count_12_mon	7
Credit_Limit	6205
Total_Revolving_Bal	1974
Avg_Open_To_Buy	6813
Total_Amt_Chng_Q4_Q1	1158
Total_Trans_Amt	5033
Total_Trans_Ct	126
Total_Ct_Chng_Q4_Q1	830
Avg_Utilization_Ratio	964

dtype: int64

# Check the number of Null values in each column
data.isna().sum()

	0
CLIENTNUM	0
Attrition_Flag	0
Customer_Age	0
Gender	0
Dependent_count	0
Education_Level	1519
Marital_Status	749
Income_Category	0
Card_Category	0
Months_on_book	0
Total_Relationship_Count	0
Months_Inactive_12_mon	0
Contacts_Count_12_mon	0
Credit_Limit	0
Total_Revolving_Bal	0
Avg_Open_To_Buy	0
Total_Amt_Chng_Q4_Q1	0
Total_Trans_Amt	0
Total_Trans_Ct	0
Total_Ct_Chng_Q4_Q1	0
Avg_Utilization_Ratio	0

dtype: int64

• Education_Level and Marital_Status columns have significant amount of null values. 1519 and 749 respectively.
• No other column has any missing value.

• Attrition_Flag, Gender, Education_Level, Marital_Status, Income_Category, Card_Category columns have a manageable number of unique values (ranging from 2 to 6), making them suitable for encoding and analysis as categorical features.
• Credit_Limit, Total_Revolving_Bal, Avg_Open_To_Buy, Total_Amt_Chng_Q4_Q1, Total_Trans_Amt, Total_Trans_Ct, Total_Ct_Chng_Q4_Q1, Avg_Utilization_Ratio columns have a large number of unique values, confirming they are continuous numerical features that may benefit from scaling or normalization.

Dropping Irrelevant or Redundant Columns

• CLIENTNUM: Since this is a unique identifier for each customer, it doesn’t add any predictive value.
• Avg_Open_To_Buy: This column represents the difference between Credit_Limit and Total_Revolving_Bal, so it may be redundant.

# Dropping two columns from the dataframe
data.drop(columns=['CLIENTNUM','Avg_Open_To_Buy'], inplace=True)

Handling Missing Values

• Education_Level: Since this is a categorical feature with a significant number of missing values (1519 out of 10127), we will use imputation to replace missing values by filling them based on a related column Income_Category.
• Marital_Status: This column has fewer missing values (749 out of 10127), so imputing with the mode is likely effective.

# Display unique values in Income_Category and Education_Level
print("Unique values in Income_Category:", data['Income_Category'].unique())
print("Unique values in Education_Level:", data['Education_Level'].unique())
print("Unique values in Marital_Status:", data['Marital_Status'].unique())

Unique values in Income_Category: ['$60K - $80K' 'Less than $40K' '$80K - $120K' '$40K - $60K' '$120K +'
 'abc']
Unique values in Education_Level: ['High School' 'Graduate' 'Uneducated' nan 'College' 'Post-Graduate'
 'Doctorate']
Unique values in Marital_Status: ['Married' 'Single' nan 'Divorced']

• We will consider 'abc' in Income_Category as missing value.
• 'nan' in Education_Level is missing value.

# Replace the placeholder value 'abc' in Income_Category with NaN to indicate missing data
data['Income_Category'] = data['Income_Category'].replace('abc', np.nan)

# Fill missing values in Marital_Status with its mode
data['Marital_Status'].fillna(data['Marital_Status'].mode()[0], inplace=True)

# Impute missing Education_Level values based on the mode in each Income_Category group
data['Education_Level'] = data['Education_Level'].fillna(data.groupby('Income_Category')['Education_Level'].transform(lambda x: x.mode()[0] if not x.mode().empty else "High School"))

💡We could use sklearn's KNNImputer here to handle missing values by imputing them based on the k-nearest neighbors in the dataset, which replaces each missing value with a weighted average (or the most frequent value in the case of categorical data) of its k nearest neighbors, calculated using a specified distance metric. But, we took a manual approach above.

Data Type Corrections

• We will convert the following columns to categorical types: Attrition_Flag, Gender, Education_Level, Marital_Status, Income_Category, Card_Category
• This conversion will allow for more efficient memory usage and simplify encoding for machine learning models.

# Convert specified columns to categorical types
categorical_columns = ['Attrition_Flag', 'Gender', 'Education_Level', 'Marital_Status', 'Income_Category', 'Card_Category']
data[categorical_columns] = data[categorical_columns].astype('category')

# Let's see data types and non-null counts one more time
data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 10127 entries, 0 to 10126
Data columns (total 19 columns):
 #   Column                    Non-Null Count  Dtype   
---  ------                    --------------  -----   
 0   Attrition_Flag            10127 non-null  category
 1   Customer_Age              10127 non-null  int64   
 2   Gender                    10127 non-null  category
 3   Dependent_count           10127 non-null  int64   
 4   Education_Level           9982 non-null   category
 5   Marital_Status            10127 non-null  category
 6   Income_Category           9015 non-null   category
 7   Card_Category             10127 non-null  category
 8   Months_on_book            10127 non-null  int64   
 9   Total_Relationship_Count  10127 non-null  int64   
 10  Months_Inactive_12_mon    10127 non-null  int64   
 11  Contacts_Count_12_mon     10127 non-null  int64   
 12  Credit_Limit              10127 non-null  float64 
 13  Total_Revolving_Bal       10127 non-null  int64   
 14  Total_Amt_Chng_Q4_Q1      10127 non-null  float64 
 15  Total_Trans_Amt           10127 non-null  int64   
 16  Total_Trans_Ct            10127 non-null  int64   
 17  Total_Ct_Chng_Q4_Q1       10127 non-null  float64 
 18  Avg_Utilization_Ratio     10127 non-null  float64 
dtypes: category(6), float64(4), int64(9)
memory usage: 1.1 MB

• No more object types. Every feature is either categorical or numerical.
  💡 Notice that the memory usage has decreased from 1.6+ MB to 1.1 MB, this technique is generally useful for bigger datasets.
  💡 We have 19 columns now, including 6 categorical columns and no more missing value.

Data Summary

Statistical Summary of the Numerical Features

# Summary of continuous columns (numerical values only)
data[['Customer_Age','Dependent_count','Months_on_book','Total_Relationship_Count','Months_Inactive_12_mon','Contacts_Count_12_mon','Credit_Limit','Total_Revolving_Bal','Total_Amt_Chng_Q4_Q1','Total_Trans_Amt','Total_Trans_Ct','Total_Ct_Chng_Q4_Q1','Avg_Utilization_Ratio']].describe().T

	count	mean	std	min	25%	50%	75%	max
Customer_Age	10127.000	46.326	8.017	26.000	41.000	46.000	52.000	73.000
Dependent_count	10127.000	2.346	1.299	0.000	1.000	2.000	3.000	5.000
Months_on_book	10127.000	35.928	7.986	13.000	31.000	36.000	40.000	56.000
Total_Relationship_Count	10127.000	3.813	1.554	1.000	3.000	4.000	5.000	6.000
Months_Inactive_12_mon	10127.000	2.341	1.011	0.000	2.000	2.000	3.000	6.000
Contacts_Count_12_mon	10127.000	2.455	1.106	0.000	2.000	2.000	3.000	6.000
Credit_Limit	10127.000	8631.954	9088.777	1438.300	2555.000	4549.000	11067.500	34516.000
Total_Revolving_Bal	10127.000	1162.814	814.987	0.000	359.000	1276.000	1784.000	2517.000
Total_Amt_Chng_Q4_Q1	10127.000	0.760	0.219	0.000	0.631	0.736	0.859	3.397
Total_Trans_Amt	10127.000	4404.086	3397.129	510.000	2155.500	3899.000	4741.000	18484.000
Total_Trans_Ct	10127.000	64.859	23.473	10.000	45.000	67.000	81.000	139.000
Total_Ct_Chng_Q4_Q1	10127.000	0.712	0.238	0.000	0.582	0.702	0.818	3.714
Avg_Utilization_Ratio	10127.000	0.275	0.276	0.000	0.023	0.176	0.503	0.999

• Customer_Age: The average age is ~46 years, with a median of 46, indicating a fairly symmetric distribution. Ages range from 26 to 73.
• Dependent_count: The average number of dependents is ~2.35, w/ most customers having between 1 and 3. The range is from 0 to 5.
• Months_on_book: Customers have an average of 36 months on record with the bank, w/ a median of 36. Most customers fall within 31 to 40 months, with a maximum of 56 months.
• Total_Relationship_Count: This feature, with an average of around 3.8, shows that customers generally have 3 to 5 relationships with the bank, indicating customer engagement levels.
• Months_Inactive_12_mon: The mean and median are both around 2, with most customers having between 2 to 3 inactive months in the past year, suggesting moderate account inactivity.
• Contacts_Count_12_mon: With a mean of 2.46 and a median of 2, most customers have been in contact 2 to 3 times in the last year, showing low customer support engagement.
• Credit_Limit: This variable has a high mean (around 8632) but with substantial variability (standard deviation ~9089), indicating some customers have exceptionally high credit limits. The range extends from 1438 to 34,516.
• Total_Revolving_Bal: Customers have an average revolving balance of ~1163, with a moderate distribution between 359 and 1784, and a maximum of 2517.
• Total_Amt_Chng_Q4_Q1: The average change in transaction amount from Q4 to Q1 is ~0.76, with a median of 0.74, suggesting most customers have relatively stable spending, though some variability exists.
• Total_Trans_Amt: The mean transaction amount is ~4400, though with a high standard deviation, indicating a wide range of spending behavior among customers, from 510 to 18,484.
• Total_Trans_Ct: Customers have an average transaction count of around 65, with a median of 67, which suggests that most are moderately active. The range is from 10 to 139 transactions.
• Total_Ct_Chng_Q4_Q1: The average change in transaction count from Q4 to Q1 is about 0.71, with a median of 0.70, indicating relatively stable transaction frequency for most customers.
• Avg_Utilization_Ratio: With an average of around 0.27 and a median of 0.18, utilization ratios are generally low, but there’s variability with some customers reaching up to nearly full credit usage (0.999).

Summary of the Categorical Features

# Summary for the non-nemerical features
data.describe(include=['category']).T

	count	unique	top	freq
Attrition_Flag	10127	2	Existing Customer	8500
Gender	10127	2	F	5358
Education_Level	9982	6	Graduate	4502
Marital_Status	10127	3	Married	5436
Income_Category	9015	5	Less than $40K	3561
Card_Category	10127	4	Blue	9436

• Attrition_Flag: Out of 10,127 records, the majority of customers (8,500) are classified as "Existing Customers," indicating that most customers have retained their accounts.
• Gender: The gender distribution shows a slight majority of female (F) customers, with 5,358 out of the total.
• Education_Level: Most customers (4,647) are listed as "Graduate," suggesting that the largest group of customers has a graduate-level education.
• Marital_Status: A significant proportion of customers (5,436) are married, indicating that marital status may correlate with other factors in this dataset.
• Income_Category: The most common income level among customers is "Less than $40K," with 3,561 individuals, showing that a considerable number of customers fall into a lower income category.
• Card_Category: The vast majority of customers (9,436) hold the "Blue" card, suggesting this card type is the most common, possibly the entry-level or standard card offered by the bank.

Target Variable Identification

OUR TARGET VARIABLE IS: Attrition_Flag

• The goal of the analysis is to predict which customers are likely to leave Thera Bank's credit card services, as this attrition could lead to a loss of revenue from various fees associated with card use.
• By using Attrition_Flag as the target variable, the classification model can help Thera Bank identify customers at risk of leaving, allowing the bank to intervene and improve service areas that may prevent customer attrition.

# Continue with this dataframe
df = data.copy()

Exploratory Data Analysis (EDA)

We will seek answers to the following questions through exploratory data analysis (EDA).

Questions:

1. How is the total transaction amount distributed?
2. What is the distribution of the level of education of customers?
3. What is the distribution of the level of income of customers?
4. How does the change in transaction amount between Q4 and Q1 (total_ct_change_Q4_Q1) vary by the customer's account status (Attrition_Flag)?
5. How does the number of months a customer was inactive in the last 12 months (Months_Inactive_12_mon) vary by the customer's account status (Attrition_Flag)?
6. What are the attributes that have a strong correlation with each other?

Functions Needed for Exploratory Data Analysis

# Function to plot a boxplot and a histogram along the same scale

def histogram_boxplot(data, feature, figsize=(12, 7), kde=False, bins=None):
    """
    Boxplot and histogram combined

    data: dataframe
    feature: dataframe column
    figsize: size of figure (default (12,7))
    kde: whether to the show density curve (default False)
    bins: number of bins for histogram (default None)
    """
    f2, (ax_box2, ax_hist2) = plt.subplots(
        nrows=2,  # number of rows of the subplot grid= 2
        sharex=True,  # x-axis will be shared among all subplots
        gridspec_kw={"height_ratios": (0.25, 0.75)},
        figsize=figsize,
    )  # creating the 2 subplots
    sns.boxplot(
        data=data, x=feature, ax=ax_box2, showmeans=True, color="violet"
    )  # boxplot will be created and a triangle will indicate the mean value of the column
    sns.histplot(
        data=data, x=feature, kde=kde, ax=ax_hist2, bins=bins, palette="winter"
    ) if bins else sns.histplot(
        data=data, x=feature, kde=kde, ax=ax_hist2
    )  # for histogram
    ax_hist2.axvline(
        data[feature].mean(), color="green", linestyle="--"
    )  # add mean to the histogram
    ax_hist2.axvline(
        data[feature].median(), color="black", linestyle="-"
    )  # add median to the histogram

# Function to create labeled barplots

def labeled_barplot(data, feature, perc=False, n=None):
    """
    Barplot with percentage at the top

    data: dataframe
    feature: dataframe column
    perc: whether to display percentages instead of count (default is False)
    n: displays the top n category levels (default is None, i.e., display all levels)
    """

    total = len(data[feature])  # length of the column
    count = data[feature].nunique()
    if n is None:
        plt.figure(figsize=(count + 1, 5))
    else:
        plt.figure(figsize=(n + 1, 5))

    plt.xticks(rotation=90, fontsize=15)
    ax = sns.countplot(
        data=data,
        x=feature,
        palette="Paired",
        order=data[feature].value_counts().index[:n].sort_values(),
    )

    for p in ax.patches:
        if perc == True:
            label = "{:.1f}%".format(
                100 * p.get_height() / total
            )  # percentage of each class of the category
        else:
            label = p.get_height()  # count of each level of the category

        x = p.get_x() + p.get_width() / 2  # width of the plot
        y = p.get_height()  # height of the plot

        ax.annotate(
            label,
            (x, y),
            ha="center",
            va="center",
            size=12,
            xytext=(0, 5),
            textcoords="offset points",
        )  # annotate the percentage

    plt.show()  # show the plot

Univariate Analysis

# Visualize the distribution of the total transaction amount using a histogram and boxplot
histogram_boxplot(df,'Total_Trans_Amt')

• The total transaction amount follows a right-skewed distribution, with most observations ranging between 1,200 and 5,000.
• The average is approximately 4,400, while the median is around 3,900.
• There are many outliers present in this variable.
• There is a wide range of spending behavior among customers.

# Visualize the distribution of the credit limit using a histogram and boxplot
histogram_boxplot(df,'Credit_Limit')

• The distribution of the Credit_Limit is skewed to the right.
• 50% of the customers of the bank have a credit limit of less than 5000.
• Many customers have a maximum Credit Limit of 35,000, likely due to a cap in the data. Customers with higher limits might all be grouped under this cap, suggesting the presence of limits beyond 35,000 that are not explicitly recorded.

# Visualize the distribution of customers' education levels using a labeled bar plot
labeled_barplot(df,'Education_Level', perc=True)

• Customers with graduate-level education make up the highest percentage of observations, at 45.9%.
• High school graduates and customers with no education account for 19.9% and 14.7% of the observations, respectively.
• A total of 85.3% of customers have some level of education.

# Visualize the distribution of customers' income level using a labeled bar plot
labeled_barplot(df,'Income_Category', perc=True)

• Customers with an income of less than 40K represent the highest percentage of observations, at 35.2%.
• A significant portion of customers fall into the lower income category.
• Only 7.2% of customers earn above 120K.
• Also, 11% of customers have missing income category information, which is not displayed in the chart above as it's been coded as NaN.

# Visualize the distribution of our classes using a labeled bar plot
labeled_barplot(df,'Attrition_Flag', perc=True)

Bivariate Analysis

Attrition vs Customer Transactions & Limits

cols = df[['Total_Ct_Chng_Q4_Q1','Total_Trans_Amt','Total_Trans_Ct','Total_Amt_Chng_Q4_Q1','Credit_Limit', 'Total_Revolving_Bal']].columns.tolist()
plt.figure(figsize=(14,7))

for i, variable in enumerate(cols):
                     plt.subplot(2, 3, i+1)
                     sns.boxplot(x=df["Attrition_Flag"], y=df[variable], palette="PuBu")
                     plt.tight_layout()
                     plt.title(variable)
plt.show()

• Customers with low change in transaction count between Q4 and Q1 are more likely to attrite.
• Lower total transaction amount and fewer transactions over the past year are also associated with higher attrition.
• Customers holding a smaller total revolving balance are especially likely to attrite.
• Although the impact is less significant, lower assigned credit limits are slightly linked to increased attrition.

Attrition vs Customer Engagement

cols = df[['Months_Inactive_12_mon','Contacts_Count_12_mon','Months_on_book','Total_Relationship_Count','Avg_Utilization_Ratio']].columns.tolist()
plt.figure(figsize=(14,7))

for i, variable in enumerate(cols):
                     plt.subplot(2, 3, i+1)
                     sns.boxplot(x=df["Attrition_Flag"], y=df[variable], palette="PuBu")
                     plt.tight_layout()
                     plt.title(variable)
plt.show()

• Customers with an average number of inactive months are more likely to stay.
• Customers with a slightly lower average contact count in the last 12 months are more likely to stay.
• Customers with a low average number of total relationships are more likely to attrite, which seems to contradict the first two observations.
• Customers with a low average utilization ratio are more likely to attrite.
• The length of the relationship with the bank (in months) does not appear to impact customer attrition.

Attrition vs Customer Demographics

cols = df[['Customer_Age','Dependent_count']].columns.tolist()
plt.figure(figsize=(14,7))

for i, variable in enumerate(cols):
                     plt.subplot(2, 3, i+1)
                     sns.boxplot(x=df["Attrition_Flag"], y=df[variable], palette="PuBu")
                     plt.tight_layout()
                     plt.title(variable)
plt.show()

• Younger customers are slightly more likely to stay, although the difference is almost insignificant.
• Customers with one dependent are more likely to stay.

Correlation Check

plt.figure(figsize=(15,11))
sns.heatmap(df.corr(numeric_only=True), annot=True, vmin=-1, vmax=1, fmt='.2f', cmap="Spectral")
plt.show()

• Total transaction amount and total transaction count are highly positively correlated.
• Months on book and customer age show a strong positive correlation.
• Total revolving balance and average utilization ratio are positively correlated.
• Total amount change between Q4 and Q1 and total count change between Q4 and Q1 have a moderate correlation.
• Average utilization ratio and credit limit have a negative correlation.
• Total transaction amount and total relationship count are also negatively correlated.

👉 Given the correlation ratios above, it’s reasonable to consider dropping highly correlated pairs (0.80 and above) to reduce redundancy in the machine learning model:

• Total Transaction Amount and Total Transaction Count (0.81): These metrics are both strong indicators of customer activity. Retaining both may add little value, so keeping only one or creating a composite feature could be more effective.
• Months on Book and Customer Age (0.79): These variables overlap as indicators of the customer's tenure or history with the bank. Representing them with a single variable may capture the necessary information while reducing redundancy.
• Customer age showed no significant effect on the Attrition Flag in the bivariate analysis above.

# Drop the columns 'Total_Trans_Ct' and 'Customer_Age' from the DataFrame
df = df.drop(['Total_Trans_Ct', 'Customer_Age'], axis=1)

Data Preprocessing

Outlier Detection and Treatment

# Outlier detection using boxplots
numeric_columns = df.select_dtypes(include=np.number).columns.tolist()

plt.figure(figsize=(13, 11))

for i, variable in enumerate(numeric_columns):
    plt.subplot(6, 4, i+1)
    plt.boxplot(df[variable], whis=1.5)
    plt.tight_layout()
    plt.title(variable)

plt.show()

• There are many outliers in the Change in Transaction Amount, Change in Transaction Count between Q4 and Q1, and Credit Limit data.
• However, we will not treat these outliers, as they are considered valid values and will contribute to the model.

Data Prep for Model Building

# Encode the Attrition_Flag column by replacing Existing and Attrited customers to 0 and 1 respectively
df["Attrition_Flag"].replace("Existing Customer", 0, inplace=True)
df["Attrition_Flag"].replace("Attrited Customer", 1, inplace=True)

# Separating the features and the target variable
X = df.drop(['Attrition_Flag'], axis=1)   # independent variables
y = df['Attrition_Flag']   # dependent variable

# Splitting the data into training, validation and test sets (A 70/15/15 split)

# Step 1: Split the data into training (70%) and temporary (30%)
X_train, X_temp, y_train, y_temp = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)

# Step 2: Split the temporary set into validation (50% of temp -> 15% of total) and test (50% of temp -> 15% of total)
X_val, X_test, y_val, y_test = train_test_split(X_temp, y_temp, test_size=0.5, random_state=0, stratify=y_temp)

print(X_train.shape, X_val.shape, X_test.shape)

(7088, 16) (1519, 16) (1520, 16)

• We have 7,088 observations in the train set, 1,519 observations in the validation set, and 1,520 observations in the test set.
• This is a 70% train, 15% validation, and 15% test distribution.

# Check class balance for whole data, train set, validation set, and test set

print("Target value ratio in y")
print(y.value_counts(1))
print("*" * 80)
print("Target value ratio in y_train")
print(y_train.value_counts(1))
print("*" * 80)
print("Target value ratio in y_val")
print(y_val.value_counts(1))
print("*" * 80)
print("Target value ratio in y_test")
print(y_test.value_counts(1))
print("*" * 80)

Target value ratio in y
Attrition_Flag
0   0.839
1   0.161
Name: proportion, dtype: float64
********************************************************************************
Target value ratio in y_train
Attrition_Flag
0   0.839
1   0.161
Name: proportion, dtype: float64
********************************************************************************
Target value ratio in y_val
Attrition_Flag
0   0.839
1   0.161
Name: proportion, dtype: float64
********************************************************************************
Target value ratio in y_test
Attrition_Flag
0   0.839
1   0.161
Name: proportion, dtype: float64
********************************************************************************

• The same proportion of each class in both the training and test sets (84% and 16%)

Encoding Categorical Features for All Three Sets

# Converting categorical variables into dummy variables
X_train = pd.get_dummies(X_train, drop_first=True, dtype=int)
X_val = pd.get_dummies(X_val, drop_first=True, dtype=int)
X_test = pd.get_dummies(X_test, drop_first=True, dtype=int)
print(X_train.shape, X_val.shape, X_test.shape)

(7088, 26) (1519, 26) (1520, 26)

• Notice that the number of columns increased from 16 to 26 with the dummies addition.

X_train.head()

	Dependent_count	Months_on_book	Total_Relationship_Count	Months_Inactive_12_mon	Contacts_Count_12_mon	Credit_Limit	Total_Revolving_Bal	Total_Amt_Chng_Q4_Q1	Total_Trans_Amt	Total_Ct_Chng_Q4_Q1	...	Education_Level_Uneducated	Marital_Status_Married	Marital_Status_Single	Income_Category_$40K - $60K	Income_Category_$60K - $80K	Income_Category_$80K - $120K	Income_Category_Less than $40K	Card_Category_Gold	Card_Category_Platinum	Card_Category_Silver
1798	1	55	5	1	4	2661.000	1466	0.631	1354	0.316	...	0	1	0	0	0	0	1	0	0	0
4696	2	41	5	2	3	14902.000	0	0.312	2038	0.625	...	0	1	0	0	0	1	0	0	0	0
6448	2	35	3	3	2	2793.000	1517	0.712	3485	0.895	...	0	1	0	0	0	0	1	0	0	0
3230	0	36	3	3	2	6853.000	1679	0.975	3543	1.000	...	0	0	0	0	0	0	0	0	0	0
6293	3	33	3	3	3	2542.000	0	0.587	4087	0.780	...	0	1	0	0	0	0	0	0	0	0

5 rows × 26 columns

• This is how our dataset looks like after dummies

Model Building

Model evaluation criterion

The nature of predictions made by the classification model will translate as follows:

• True positives (TP) are failures correctly predicted by the model.
• False negatives (FN) are real failures in a generator where there is no detection by model.
• False positives (FP) are failure detections in a generator where there is no failure.

Which metric to optimize?

• The bank’s primary goal is to prevent customers from leaving by identifying those who are likely to cancel their credit cards.
• In this case, the best evaluation metric would be Recall.
• High Recall (minimizing False Negatives) would mean the model is able to identify as many at-risk customers as possible. If we fail to identify these customers (FN), the bank cannot take steps to retain them, which is costly.
• If we also want to consider Precision (avoiding too many False Positives, where customers are incorrectly flagged as at-risk), then F1 Score becomes useful. But, in this case, we will go with Recall.

Functions needed to output different metrics

# Defining a function to compute different metrics to check performance of a classification model built using sklearn
def model_performance_classification_sklearn(model, predictors, target):
    """
    Function to compute different metrics to check classification model performance

    model: classifier
    predictors: independent variables
    target: dependent variable
    """

    # Predicting using the independent variables
    pred = model.predict(predictors)

    acc = accuracy_score(target, pred)  # to compute Accuracy
    recall = recall_score(target, pred)  # to compute Recall
    precision = precision_score(target, pred)  # to compute Precision
    f1 = f1_score(target, pred)  # to compute F1-score

    # Creating a dataframe of metrics
    df_perf = pd.DataFrame(
        {
            "Accuracy": acc,
            "Recall": recall,
            "Precision": prec,
            "F1": f1
        },
        index=[0],
    )

    return df_perf

def confusion_matrix_sklearn(model, predictors, target):
    """
    To plot the confusion_matrix with percentages

    model: classifier
    predictors: independent variables
    target: dependent variable
    """
    y_pred = model.predict(predictors)
    cm = confusion_matrix(target, y_pred)
    labels = np.asarray(
        [
            ["{0:0.0f}".format(item) + "\n{0:.2%}".format(item / cm.flatten().sum())]
            for item in cm.flatten()
        ]
    ).reshape(2, 2)

    plt.figure(figsize=(6, 4))
    sns.heatmap(cm, annot=labels, fmt="")
    plt.ylabel("True label")
    plt.xlabel("Predicted label")

Building and Evaluating Classification Models

As a requirement for this project, we will build six classification models using three variations of the dataset: the original data (baseline data), oversampled data, and undersampled data. This results in a total of 18 models.

Classification Models to be Built:

1. Decision Tree
2. Bagging Classifier
3. Random Forest
4. AdaBoost
5. Gradient Boosting
6. XGBoost

After building these 18 models, we will:

1. Evaluate all models based on performance metrics.
2. Select the 3 best-performing models among them.
3. Hyper-tune the selected 3 models to optimize performance.

Model Building with baseline data

• We are adding Logistic Regression to our model list. Logistic Regression is the simplest classification model, which we will use as a baseline to evaluate other models' performance. It sets a benchmark for comparison, and if more complex models do not significantly outperform it, the added complexity may be unnecessary for the given problem."

# List to store all the models
models = []
models_trained = []

# Appending models to the list
models.append(("Logistic Regression", LogisticRegression(solver="newton-cg", class_weight="balanced", random_state=1)))
models.append(("Decision Tree", DecisionTreeClassifier(random_state=1)))
models.append(("Bagging", BaggingClassifier(random_state=1)))
models.append(("Random Forest", RandomForestClassifier(random_state=1)))
models.append(("AdaBoost", AdaBoostClassifier(random_state=1)))
models.append(("Gradient Boosting", GradientBoostingClassifier(random_state=1)))
models.append(("XGBoost", XGBClassifier(use_label_encoder=False, eval_metric='logloss', random_state=1)))

# Training and evaluating models on the training and validation data
print("\nTraining Performance (Recall) on Baseline Data:\n")
for name, model in models:
    model.fit(X_train, y_train)
    scores = recall_score(y_train, model.predict(X_train))
    models_trained.append((name, model))
    print(f"{name}: {scores:.4f}")

print("\nValidation Performance (Recall) on Baseline Data:\n")
for name, model in models_trained:
    scores_val = recall_score(y_val, model.predict(X_val))
    print(f"{name}: {scores_val:.4f}")

Training Performance (Recall) on Baseline Data:

Logistic Regression: 0.7700
Decision Tree: 1.0000
Bagging: 0.9763
Random Forest: 1.0000
AdaBoost: 0.7428
Gradient Boosting: 0.8218
XGBoost: 1.0000

Validation Performance (Recall) on Baseline Data:

Logistic Regression: 0.7459
Decision Tree: 0.7008
Bagging: 0.7172
Random Forest: 0.6885
AdaBoost: 0.7172
Gradient Boosting: 0.7541
XGBoost: 0.8074

• The high recall scores on training data (1.00 for Decision Tree, Random Forest, and XGBoost) indicate that these models are capturing all positive cases on the training set, suggesting a likelihood of overfitting.
• Models like Bagging and Gradient Boosting show slightly lower training recall (0.97 and 0.82), indicating they might generalize better compared to the models with perfect training recall.
• In the validation set, recall scores vary significantly:
  • XGBoost performs the best with a recall of 0.80, indicating strong performance on unseen data.
  • Gradient Boosting (0.75), AdaBoost and Bagging (0.71) also show competitive recall, with better generalization compared to Decision Tree or Random Forest.
  • Decision Tree and Random Forest exhibit a notable drop in recall from training to validation, confirming overfitting trends for these models.
• Only XGBoost and Gradient Boosting perform better than Logistic Regression on the validation dataset.
• Overall, XGBoost stands out for its combination of strong recall on both training and validation sets, making it a promising candidate for further optimization.

Model Building with Oversampled data

The idea of oversampling in data evaluation is to address issues related to class imbalance, where one class (usually the minority class —attrited customers in this case—) is underrepresented in the dataset compared to the other class (majority class).

# Synthetic Minority Over Sampling Technique (SMOTE)
sm = SMOTE(sampling_strategy=1, k_neighbors=5, random_state=1)
X_train_over, y_train_over = sm.fit_resample(X_train, y_train)

print("Before OverSampling, count of label '1': {}".format(sum(y_train == 1)))
print("Before OverSampling, count of label '0': {} \n".format(sum(y_train == 0)))

print("After OverSampling, count of label '1': {}".format(sum(y_train_over == 1)))
print("After OverSampling, count of label '0': {} \n".format(sum(y_train_over == 0)))

print("After OverSampling, the shape of train_X: {}".format(X_train_over.shape))
print("After OverSampling, the shape of train_y: {} \n".format(y_train_over.shape))

Before OverSampling, count of label '1': 1139
Before OverSampling, count of label '0': 5949 

After OverSampling, count of label '1': 5949
After OverSampling, count of label '0': 5949 

After OverSampling, the shape of train_X: (11898, 26)
After OverSampling, the shape of train_y: (11898,) 

• Synthetic data was generated to balance the class distribution after oversampling.

# List to store all the models
models = []
models_trained = []

# Appending models to the list
models.append(("Logistic Regression", LogisticRegression(solver="newton-cg", class_weight="balanced", random_state=1)))
models.append(("Decision Tree", DecisionTreeClassifier(random_state=1)))
models.append(("Bagging", BaggingClassifier(random_state=1)))
models.append(("Random Forest", RandomForestClassifier(random_state=1)))
models.append(("AdaBoost", AdaBoostClassifier(random_state=1)))
models.append(("Gradient Boosting", GradientBoostingClassifier(random_state=1)))
models.append(("XGBoost", XGBClassifier(use_label_encoder=False, eval_metric='logloss', random_state=1)))

# Training and evaluating models on the training (oversampled) and validation data
print("\nTraining Performance (Recall) on Oversampled Data:\n")
for name, model in models:
    model.fit(X_train_over, y_train_over)
    scores = recall_score(y_train_over, model.predict(X_train_over))
    models_trained.append((name, model))
    print(f"{name}: {scores:.4f}")

print("\nValidation Performance (Recall) on Oversampled Data:\n")
for name, model in models_trained:
    scores_val = recall_score(y_val, model.predict(X_val))
    print(f"{name}: {scores_val:.4f}")

Training Performance (Recall) on Oversampled Data:

Logistic Regression: 0.8546
Decision Tree: 1.0000
Bagging: 0.9965
Random Forest: 1.0000
AdaBoost: 0.9365
Gradient Boosting: 0.9660
XGBoost: 0.9998

Validation Performance (Recall) on Oversampled Data:

Logistic Regression: 0.5697
Decision Tree: 0.7664
Bagging: 0.7869
Random Forest: 0.7664
AdaBoost: 0.7377
Gradient Boosting: 0.8238
XGBoost: 0.8361

• The perfect recall scores in training for Decision Tree, Random Forest, and XGBoost (1.00 or close) indicate that these models capture nearly all positive cases in the oversampled training data, raising concerns about overfitting for these models.
• Bagging, Gradient Boosting, and AdaBoost exhibit slightly lower recall in training (0.99 to 0.93), suggesting a more balanced generalization compared to the models with perfect recall.
• On the validation set, recall scores for all models drop noticeably compared to training. This indicates these models may struggle more with unseen data, further supporting the overfitting concern.
• Gradient Boosting and XGBoost maintain the highest validation recall (0.82 & 0.83) with relatively smaller recall drops, making them strong candidates for further tuning and practical application.

Model Building with Undersampled data

Similarly, the idea of undersampling in data evaluation is to address issues related to class imbalance, where the majority class is overrepresented in the dataset compared to the minority class. By undersampling the majority class, we reduce the number of majority class samples to make the dataset more balanced.

# Random undersampler for under sampling the data (Random Undersampling)
rus = RandomUnderSampler(random_state=1, sampling_strategy=1)
X_train_un, y_train_un = rus.fit_resample(X_train, y_train)

print("Before Under Sampling, count of label '1': {}".format(sum(y_train == 1)))
print("Before Under Sampling, count of label '0': {} \n".format(sum(y_train == 0)))

print("After Under Sampling, count of label '1': {}".format(sum(y_train_un == 1)))
print("After Under Sampling, count of label '0': {} \n".format(sum(y_train_un == 0)))

print("After Under Sampling, the shape of train_X: {}".format(X_train_un.shape))
print("After Under Sampling, the shape of train_y: {} \n".format(y_train_un.shape))

Before Under Sampling, count of label '1': 1139
Before Under Sampling, count of label '0': 5949 

After Under Sampling, count of label '1': 1139
After Under Sampling, count of label '0': 1139 

After Under Sampling, the shape of train_X: (2278, 26)
After Under Sampling, the shape of train_y: (2278,) 

• The data was balanced by undersampling the majority class, reducing its count to match the minority class for an equal distribution.

# List to store all the models
models = []
models_trained = []

# Appending models to the list
models.append(("Logistic Regression", LogisticRegression(solver="newton-cg", class_weight="balanced", random_state=1)))
models.append(("Decision Tree", DecisionTreeClassifier(random_state=1)))
models.append(("Bagging", BaggingClassifier(random_state=1)))
models.append(("Random Forest", RandomForestClassifier(random_state=1)))
models.append(("AdaBoost", AdaBoostClassifier(random_state=1)))
models.append(("Gradient Boosting", GradientBoostingClassifier(random_state=1)))
models.append(("XGBoost", XGBClassifier(use_label_encoder=False, eval_metric='logloss', random_state=1)))

# Training and evaluating models on the training (undersampled) and validation data
print("\nTraining Performance (Recall) on Undersampled Data:\n")
for name, model in models:
    model.fit(X_train_un, y_train_un)
    scores_train = recall_score(y_train_un, model.predict(X_train_un))
    models_trained.append((name, model))
    print(f"{name}: {scores_train:.4f}")

print("\nValidation Performance (Recall) on Undersampled Data:\n")
for name, model in models_trained:
    scores = recall_score(y_val, model.predict(X_val))
    print(f"{name}: {scores:.4f}")

Training Performance (Recall) on Undersampled Data:

Logistic Regression: 0.7709
Decision Tree: 1.0000
Bagging: 0.9895
Random Forest: 1.0000
AdaBoost: 0.9350
Gradient Boosting: 0.9701
XGBoost: 1.0000

Validation Performance (Recall) on Undersampled Data:

Logistic Regression: 0.7582
Decision Tree: 0.8852
Bagging: 0.9057
Random Forest: 0.9139
AdaBoost: 0.9016
Gradient Boosting: 0.9303
XGBoost: 0.9344

• The training recall for Decision Tree, Random Forest, and XGBoost is perfect (1.00), indicating these models are fully capturing positive cases in the undersampled training data, but it also raises a concern of overfitting.
• Bagging, Gradient Boosting, and AdaBoost exhibit slightly lower but still high training recall (0.98 to 0.93), suggesting they may generalize better than the models with perfect recall.
• On the validation set, recall values are strong across all models (ranging from 0.88 to 0.93), with Gradient Boosting and XGBoost performing the best (0.93).
• The smaller gap between training and validation recall for Gradient Boosting and XGBoost makes these models strong candidates for further evaluation and tuning, as they balance performance and generalization effectively.

Selecting the Best Three Models

1. Random Forest is selected for its consistency and strong generalization across various datasets.
2. Gradient Boosting is chosen for its balanced performance and ability to handle class imbalance effectively.
3. XGBoost is selected for its high recall performance across both training and validation data, as well as its robustness to overfitting.

These models are chosen based on their ability to handle class imbalance, generalize well to unseen data, and maintain a good balance between recall and overfitting, making them ideal candidates for further fine-tuning and optimization.

# List to store the selected best-performing models
best_models = []

# Appending the best-performing models to the list
best_models.append(("Random Forest", RandomForestClassifier(random_state=1)))
best_models.append(("Gradient Boosting", GradientBoostingClassifier(random_state=1)))
best_models.append(("XGBoost", XGBClassifier(use_label_encoder=False, eval_metric='logloss', random_state=1)))

def evaluate_models(models, X_train, y_train, X_val, y_val, description=""):
    """
    Function to train and evaluate models on given data.
    Prints the classification report and recall score for each model.
    """
    print(f"\nValidation Performance ({description}):\n")
    for name, model in models:
        # Train the model
        model.fit(X_train, y_train)

        # Predict on the validation data
        y_pred = model.predict(X_val)

        # Compute recall score and print classification report
        recall = recall_score(y_val, y_pred, average='binary')
        print(f"{name} - Recall: {recall:.4f}")
        print(f"Classification Report for {name}:\n{classification_report(y_val, y_pred)}")

# Evaluate models on baseline data
evaluate_models(best_models, X_train, y_train, X_val, y_val, description="Baseline Data")

# Evaluate models on oversampled data
evaluate_models(best_models, X_train_over, y_train_over, X_val, y_val, description="Oversampled Data")

# Evaluate models on undersampled data
evaluate_models(best_models, X_train_un, y_train_un, X_val, y_val, description="Undersampled Data")

Validation Performance (Baseline Data):

Random Forest - Recall: 0.6885
Classification Report for Random Forest:
              precision    recall  f1-score   support

           0       0.94      0.99      0.97      1275
           1       0.92      0.69      0.79       244

    accuracy                           0.94      1519
   macro avg       0.93      0.84      0.88      1519
weighted avg       0.94      0.94      0.94      1519

Gradient Boosting - Recall: 0.7541
Classification Report for Gradient Boosting:
              precision    recall  f1-score   support

           0       0.95      0.99      0.97      1275
           1       0.93      0.75      0.83       244

    accuracy                           0.95      1519
   macro avg       0.94      0.87      0.90      1519
weighted avg       0.95      0.95      0.95      1519

XGBoost - Recall: 0.8074
Classification Report for XGBoost:
              precision    recall  f1-score   support

           0       0.96      0.99      0.97      1275
           1       0.92      0.81      0.86       244

    accuracy                           0.96      1519
   macro avg       0.94      0.90      0.92      1519
weighted avg       0.96      0.96      0.96      1519


Validation Performance (Oversampled Data):

Random Forest - Recall: 0.7664
Classification Report for Random Forest:
              precision    recall  f1-score   support

           0       0.96      0.97      0.96      1275
           1       0.83      0.77      0.80       244

    accuracy                           0.94      1519
   macro avg       0.89      0.87      0.88      1519
weighted avg       0.94      0.94      0.94      1519

Gradient Boosting - Recall: 0.8238
Classification Report for Gradient Boosting:
              precision    recall  f1-score   support

           0       0.97      0.96      0.97      1275
           1       0.81      0.82      0.82       244

    accuracy                           0.94      1519
   macro avg       0.89      0.89      0.89      1519
weighted avg       0.94      0.94      0.94      1519

XGBoost - Recall: 0.8361
Classification Report for XGBoost:
              precision    recall  f1-score   support

           0       0.97      0.98      0.97      1275
           1       0.88      0.84      0.86       244

    accuracy                           0.95      1519
   macro avg       0.92      0.91      0.91      1519
weighted avg       0.95      0.95      0.95      1519


Validation Performance (Undersampled Data):

Random Forest - Recall: 0.9139
Classification Report for Random Forest:
              precision    recall  f1-score   support

           0       0.98      0.93      0.96      1275
           1       0.73      0.91      0.81       244

    accuracy                           0.93      1519
   macro avg       0.86      0.92      0.88      1519
weighted avg       0.94      0.93      0.93      1519

Gradient Boosting - Recall: 0.9303
Classification Report for Gradient Boosting:
              precision    recall  f1-score   support

           0       0.99      0.92      0.95      1275
           1       0.70      0.93      0.80       244

    accuracy                           0.93      1519
   macro avg       0.84      0.93      0.88      1519
weighted avg       0.94      0.93      0.93      1519

XGBoost - Recall: 0.9344
Classification Report for XGBoost:
              precision    recall  f1-score   support

           0       0.99      0.93      0.96      1275
           1       0.73      0.93      0.82       244

    accuracy                           0.93      1519
   macro avg       0.86      0.93      0.89      1519
weighted avg       0.94      0.93      0.94      1519

• XGBoost consistently demonstrates the highest accuracy, recall, precision, and F1-score across all data sets, making it the top performer in this evaluation.
• Gradient Boosting and Random Forest also perform strongly, with Random Forest excelling in recall for the positive class and Gradient Boosting offering a well-balanced trade-off between precision and recall.
• While XGBoost provides the most balanced and reliable performance, all three models—XGBoost, Gradient Boosting and Random Forest are excellent choices for tasks that prioritize high recall for the positive class.

Hyperparameter Tuning

Sample Parameter Grids used for Hyperparameter Tuning

• Sample parameter grids are provided below for hyperparameter tuning aim to balance model performance improvement with execution time. It's possible to adjust the grid size based on system configuration and the acceptable execution time.
  

⚠️ Please note that if the parameter grid is extended to improve the model performance further, the execution time will increase significantly.

• For Random Forest:

param_grid = {
    "n_estimators": [50,110,25],
    "min_samples_leaf": np.arange(1, 4),
    "max_features": [np.arange(0.3, 0.6, 0.1),'sqrt'],
    "max_samples": np.arange(0.4, 0.7, 0.1)
}

• For Gradient Boosting:

param_grid = {
    "init": [AdaBoostClassifier(random_state=1),DecisionTreeClassifier(random_state=1)],
    "n_estimators": np.arange(50,110,25),
    "learning_rate": [0.01,0.1,0.05],
    "subsample":[0.7,0.9],
    "max_features":[0.5,0.7,1],
}

• For XGBoost:

param_grid={'n_estimators':np.arange(50,110,25),
            'scale_pos_weight':[1,2,5],
            'learning_rate':[0.01,0.1,0.05],
            'gamma':[1,3],
            'subsample':[0.7,0.9]
}

Data and Method Preferences

💡 For hyperparameter tuning, it is generally recommended to use the oversampled data to address class imbalance. However, since our selected models showed the highest Recall on the Undersampled Data, we will use this dataset for tuning.
💡 For hyperparameter tuning, RandomizedSearchCV is more flexible, faster, and can often provide better results with less computational cost than GridSearchCV.

Tuning the Random Forest model on UnderSampled data using the RandomizedSearchCV method

%%time

# Choose the type of the classifier
model_rf2 = RandomForestClassifier(random_state=1)

# Parameters to pass in RandomSearchCV
parameters = {
    "n_estimators": [50,110,25],
    "min_samples_leaf": np.arange(1, 4),
    "max_features": [np.arange(0.3, 0.6, 0.1),'sqrt'],
    "max_samples": np.arange(0.4, 0.7, 0.1)
}

# Type of scoring used to compare parameter combinations
scorer = metrics.make_scorer(metrics.recall_score)

# Calling RandomizedSearchCV
grid_obj = RandomizedSearchCV(model_rf2, parameters, n_iter=30, scoring=scorer, cv=5, random_state=1, n_jobs=-1, verbose=2)
# Using n_iter=30 (by default, n_iter=10), so randomized search will try 30 different combinations of the hyperparameters

# Fitting parameters in RandomizedSearchCV
grid_obj.fit(X_train_un, y_train_un)

# Print the best combination of parameters
print("Best parameters are {} with CV score={}:" .format(grid_obj.best_params_, grid_obj.best_score_))

Fitting 5 folds for each of 30 candidates, totalling 150 fits
Best parameters are {'n_estimators': 110, 'min_samples_leaf': 1, 'max_samples': 0.6, 'max_features': 'sqrt'} with CV score=0.9025156503593786:
CPU times: user 936 ms, sys: 57.2 ms, total: 993 ms
Wall time: 17.7 s

# The best cross-validation score achieved with these parameters
grid_obj.best_score_

0.9025156503593786

• 90.2% Recall suggests the model has been well-tuned and is performing efficiently with the selected hyperparameters.

Build a Random Forest model with Obtained Best Parameters

# Set the best combination of parameters using grid_obj.best_params_
model_rf2_tuned = RandomForestClassifier(
    random_state=1,
    n_estimators=grid_obj.best_params_['n_estimators'],
    min_samples_leaf=grid_obj.best_params_['min_samples_leaf'],
    max_samples=grid_obj.best_params_['max_samples'],
    max_features=grid_obj.best_params_['max_features']
)

# Fit the best algorithm to the data
model_rf2_tuned.fit(X_train_un, y_train_un)

RandomForestClassifier(max_samples=0.6, n_estimators=110, random_state=1)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Check the Performance of the Random Forest Model on the Test Set

# Checking recall score on test set
rf2_recall = recall_score(y_test, model_rf2_tuned.predict(X_test))
print("Recall on test set:", rf2_recall)

# Checking precision score on test set
rf2_precision = precision_score(y_test, model_rf2_tuned.predict(X_test))
print("Precision on test set:", rf2_precision)

# Checking accuracy score on test set
rf2_accuracy = accuracy_score(y_test, model_rf2_tuned.predict(X_test))
print("Accuracy on test set:", rf2_accuracy)

Recall on test set: 0.9098360655737705
Precision on test set: 0.6851851851851852
Accuracy on test set: 0.9184210526315789

• The Random Forest model has a recall score of 0.91, indicating that it identifies a large proportion of positive cases, which is crucial for minimizing false negatives. Combined with a precision of 0.68, this suggests the model performs well in minimizing false positives, but it might have a relatively higher number of false positives compared to its recall.
  
  👉 Typically, test performance can be slightly lower than validation performance due to potential data leakage in the validation set.

See the Confusion Matrix for the Random Forest Model on the Test Set

# Predicting on the test set
y_pred = model_rf2_tuned.predict(X_test)

# Generating the confusion matrix
cm = confusion_matrix(y_test, y_pred)

# Plotting the confusion matrix using seaborn
plt.figure(figsize=(6, 5))
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', xticklabels=['Positive', 'Negative'], yticklabels=['Positive', 'Negative'])
plt.xlabel('Predicted')
plt.ylabel('Actual')
plt.title('Confusion Matrix')
plt.show()

• The Random Forest model is highly effective at identifying positive cases.
• The model also performs well in predicting negative cases, with only 22 False Positives.
• Overall, the model performs well, with a high percentage of correct predictions across both classes.

Tuning the Gradient Boosting model on UnderSampled data using the RandomizedSearchCV method

%%time

# Creating pipeline
model_gbm2 = GradientBoostingClassifier(random_state=1)

# Parameters to pass in RandomSearchCV
parameters = {
    "init": [AdaBoostClassifier(random_state=1), DecisionTreeClassifier(random_state=1)],
    "n_estimators": np.arange(50,110,25),
    "learning_rate": [0.01,0.1,0.05],
    "subsample":[0.7,0.9],
    "max_features":[0.5,0.7,1],
}

# Type of scoring used to compare parameter combinations
scorer = metrics.make_scorer(metrics.recall_score)

# Calling RandomizedSearchCV
grid_obj2 = RandomizedSearchCV(estimator=model_gbm2, param_distributions=parameters, n_iter=50, scoring=scorer, cv=5, random_state=1, n_jobs=-1)
# Using n_iter=50, so randomized search will try 50 different combinations of the hyperparameters

# Fitting parameters in RandomizedSearchCV
grid_obj2.fit(X_train_un, y_train_un)

# Print the best combination of parameters
print("Best parameters are {} with CV score={}:" .format(grid_obj2.best_params_, grid_obj2.best_score_))

Best parameters are {'subsample': 0.9, 'n_estimators': 100, 'max_features': 0.5, 'learning_rate': 0.1, 'init': AdaBoostClassifier(random_state=1)} with CV score=0.9358798979828427:
CPU times: user 2.54 s, sys: 261 ms, total: 2.8 s
Wall time: 1min 25s

# The best cross-validation score achieved with these parameters
grid_obj2.best_score_

0.9358798979828427

• 93.5% Recall suggests that the Gradient Boosting model has been well-tuned and is performing efficiently with the selected hyperparameters.

Build a Gradient Boosting model with Obtained Best Parameters

# Set the best combination of parameters using grid_obj2.best_params_
model_gbm2_tuned = GradientBoostingClassifier(
    random_state=1,
    subsample=grid_obj2.best_params_['subsample'],
    n_estimators=grid_obj2.best_params_['n_estimators'],
    max_features=grid_obj2.best_params_['max_features'],
    learning_rate=grid_obj2.best_params_['learning_rate'],
    init=AdaBoostClassifier(random_state=1),
)

# Fit the best algorithm to the data
model_gbm2_tuned.fit(X_train_un, y_train_un)

GradientBoostingClassifier(init=AdaBoostClassifier(random_state=1),
                           max_features=0.5, random_state=1, subsample=0.9)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Check the Performance of the Gradient Boosting Model on the Test Set

# Checking recall score on test set
gbm2_recall = recall_score(y_test, model_gbm2_tuned.predict(X_test))
print("Recall on test set:", gbm2_recall)

# Checking precision score on test set
gbm2_precision = precision_score(y_test, model_gbm2_tuned.predict(X_test))
print("Precision on test set:", gbm2_precision)

# Checking accuracy score on test set
gbm2_accuracy = accuracy_score(y_test, model_gbm2_tuned.predict(X_test))
print("Accuracy on test set:", gbm2_accuracy)

Recall on test set: 0.930327868852459
Precision on test set: 0.7206349206349206
Accuracy on test set: 0.930921052631579

• The Gradient Boosting model shows a high recall score of 0.93, indicating it effectively identifies most positive cases, which is crucial for minimizing false negatives. With a precision of 0.72, the model performs well in identifying positives but may produce a few false positives. The accuracy score of 0.93 demonstrates strong overall performance on the test set.

See the Confusion Matrix for the Gradient Boosting Model on the Test Set

# Predicting on the test set
y_pred = model_gbm2_tuned.predict(X_test)

# Generating the confusion matrix
cm = confusion_matrix(y_test, y_pred)

# Plotting the confusion matrix using seaborn
plt.figure(figsize=(6, 5))
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', xticklabels=['Positive', 'Negative'], yticklabels=['Positive', 'Negative'])
plt.xlabel('Predicted')
plt.ylabel('Actual')
plt.title('Confusion Matrix')
plt.show()

• The Gradient Boosting model is good at identifying positive cases, with only 88 missed positives.
• The model also performs well in predicting negative cases, with only 17 False Positives.
• Overall, the model demonstrates strong performance, with a high percentage of correct predictions across both classes.

Tuning the XGBoost model on UnderSampled data using the RandomizedSearchCV method

%%time

# Defining the model
model_xgb2 = XGBClassifier(random_state=1, eval_metric='logloss')

# Parameters to pass in RandomSearchCV
parameters = {'n_estimators':np.arange(50,110,25),
              'scale_pos_weight':[1,2,5],
              'learning_rate':[0.01,0.1,0.05],
              'gamma':[1,3],
              'subsample':[0.7,0.9]
}

# Type of scoring used to compare parameter combinations
scorer = metrics.make_scorer(metrics.recall_score)

# Calling RandomizedSearchCV
grid_obj3 = RandomizedSearchCV(estimator=model_xgb2, param_distributions=parameters, n_iter=50, scoring=scorer, cv=5, random_state=1, n_jobs=-1)
# Using n_iter=50, so randomized search will try 50 different combinations of the hyperparameters

# Fitting parameters in RandomizedSearchCV
grid_obj3.fit(X_train_un, y_train_un)

# Print the best combination of parameters
print("Best parameters are {} with CV score={}:" .format(grid_obj3.best_params_, grid_obj3.best_score_))

Best parameters are {'subsample': 0.7, 'scale_pos_weight': 5, 'n_estimators': 50, 'learning_rate': 0.01, 'gamma': 3} with CV score=0.9991228070175439:
CPU times: user 1.84 s, sys: 153 ms, total: 1.99 s
Wall time: 44.1 s

# The best cross-validation score achieved with these parameters
grid_obj3.best_score_

0.9991228070175439

• A 99.9% Recall indicates that the XGBoost model is highly effective at identifying positive cases, demonstrating excellent performance with the selected hyperparameters.

Build an XGBoost model with Obtained Best Parameters

# Set the best combination of parameters using grid_obj3.best_params_
model_xgb2_tuned = XGBClassifier(
    random_state=1,
    eval_metric="logloss",
    subsample=grid_obj3.best_params_['subsample'],
    scale_pos_weight=grid_obj3.best_params_['scale_pos_weight'],
    n_estimators=grid_obj3.best_params_['n_estimators'],
    learning_rate=grid_obj3.best_params_['learning_rate'],
    gamma=grid_obj3.best_params_['gamma'],
)

# Fit the best algorithm to the data
model_xgb2_tuned.fit(X_train_un, y_train_un)

XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=None, device=None, early_stopping_rounds=None,
              enable_categorical=False, eval_metric='logloss',
              feature_types=None, gamma=3, grow_policy=None,
              importance_type=None, interaction_constraints=None,
              learning_rate=0.01, max_bin=None, max_cat_threshold=None,
              max_cat_to_onehot=None, max_delta_step=None, max_depth=None,
              max_leaves=None, min_child_weight=None, missing=nan,
              monotone_constraints=None, multi_strategy=None, n_estimators=50,
              n_jobs=None, num_parallel_tree=None, random_state=1, ...)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Check the Performance of the XGBoost Model on the Test Set

# Checking recall score on test set
xgb2_recall = recall_score(y_test, model_xgb2_tuned.predict(X_test))
print("Recall on test set:", xgb2_recall)

# Checking precision score on test set
xgb2_precision = precision_score(y_test, model_xgb2_tuned.predict(X_test))
print("Precision on test set:", xgb2_precision)

# Checking accuracy score on test set
xgb2_accuracy = accuracy_score(y_test, model_xgb2_tuned.predict(X_test))
print("Accuracy on test set:", xgb2_accuracy)

Recall on test set: 1.0
Precision on test set: 0.23238095238095238
Accuracy on test set: 0.4697368421052632

• The XGBoost model shows a perfect recall score of 1.0, meaning it effectively identifies all positive cases, which is crucial for minimizing false negatives. However, with a low precision of 0.23, the model is likely generating many false positives, meaning it incorrectly classifies a large number of negatives as positives.
• The accuracy score of 0.46 suggests that the model's overall performance on the test set is relatively low, indicating room for improvement, especially in reducing false positives while maintaining recall.

See the Confusion Matrix for the XGBoost Model on the Test Set

# Predicting on the test set
y_pred = model_xgb2_tuned.predict(X_test)

# Generating the confusion matrix
cm = confusion_matrix(y_test, y_pred)

# Plotting the confusion matrix using seaborn
plt.figure(figsize=(6, 5))
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', xticklabels=['Positive', 'Negative'], yticklabels=['Positive', 'Negative'])
plt.xlabel('Predicted')
plt.ylabel('Actual')
plt.title('Confusion Matrix')
plt.show()

• The XGBoost model correctly identifies 470 actual positive cases, which is great for recall, but it misses 806 actual positive cases by predicting them as negative.
• The model does not predict any actual negatives as positives, with 0 false positives.
• Overall, the model struggles with identifying positive cases, as shown by the large number of false negatives (806), which lowers recall. However, there are no false positives, indicating it is cautious in predicting positives.

Model Comparison and Final Model Selection

# Comparing the performance of all three models
models_test_comp_df = pd.DataFrame(
    {
        "Random Forest": [rf2_recall, rf2_precision, rf2_accuracy],
        "Gradient Boosting": [gbm2_recall, gbm2_precision, gbm2_accuracy],
        "XGBoost": [xgb2_recall, xgb2_precision, xgb2_accuracy]
    },
    index=["Recall", "Precision", "Accuracy"]
)

# Display the performance comparison
print("Performance comparison:")
print(models_test_comp_df)

Performance comparison:
           Random Forest  Gradient Boosting  XGBoost
Recall             0.910              0.930    1.000
Precision          0.685              0.721    0.232
Accuracy           0.918              0.931    0.470

• Recall: The XGBoost model achieves perfect recall (1.0), indicating it identifies all positive cases, which is crucial for minimizing false negatives. However, this comes at the expense of precision.
• Precision: The Gradient Boosting model performs the best in terms of precision (0.72), striking a balance between minimizing false positives while still maintaining good recall.
• Accuracy: The Gradient Boosting model also has the highest accuracy (0.93), suggesting it performs well across both classes without being skewed by an imbalance.

The Best Model

🥇 Based on these results, Gradient Boosting appears to be the best model overall. It strikes a balance between recall, precision, and accuracy, making it the most reliable for this task.

Feature Importance

feature_names = X_train_un.columns
importances = model_gbm2_tuned.feature_importances_
indices = np.argsort(importances)

plt.figure(figsize=(11, 11))
plt.title("Feature Importances")
plt.barh(range(len(indices)), importances[indices], color="blue", align="center")
plt.yticks(range(len(indices)), [feature_names[i] for i in indices])
plt.xlabel("Related Importance")
plt.show()

Business Insights and Conclusions

Model Development:

1. Random Forest Classifier: This model showed a balanced performance with good recall (0.91), making it suitable for identifying customers likely to leave (i.e., minimizing false negatives). However, its precision (0.68) suggests there is still room for improvement in reducing false positives (identifying customers as likely to leave when they may not).

2. Gradient Boosting Classifier: This model had high recall (0.93), which is beneficial for capturing most of the customers at risk of leaving, with a reasonable precision of 0.71. It provided the best balance between performance metrics (accuracy of 0.93), making it a strong contender for deployment.

3. XGBoost Classifier: While it achieved perfect recall (1.0), the precision (0.23) was very low, meaning it flagged too many false positives. This model would not be ideal due to its low precision, despite its ability to identify all at-risk customers.

Model Selection:

After comparing the models, Gradient Boosting emerges as the best model due to its balanced performance. It strikes a good balance between recall, precision, and accuracy, making it reliable for identifying at-risk customers while minimizing unnecessary interventions.

Key Features Affecting Customer Attrition:

The analysis of customer features (such as customer transactions, customer engagement, and demographics) likely reveals the key drivers behind customer attrition. For instance, customers with higher late payment frequencies, low engagement with their accounts, or high outstanding balances are more likely to leave.

Business Impact:

• Proactive Customer Retention: By using this model, Thera Bank can identify customers who are at risk of leaving their credit card services and take corrective actions to retain them. This could include personalized offers, fee waivers, loyalty programs, or targeted communication that highlights the benefits of continued service.

• Revenue Protection: Retaining credit card customers would help stabilize the revenue streams from fees like annual charges, late payment fees, and foreign transaction fees, all of which are crucial for the bank's financial health.

• Customer Satisfaction and Service Improvement: By addressing the reasons behind customer churn (e.g., high fees, poor service, or lack of engagement), the bank can improve its offerings and customer satisfaction. This could lead to increased brand loyalty, higher customer lifetime value, and enhanced reputation in the market.

Conclusion:

The Gradient Boosting model is the most effective tool for predicting and preventing credit card attrition at Thera Bank. By identifying at-risk customers early, the bank can proactively implement retention strategies to minimize churn, protect revenues, and enhance customer satisfaction.

---