Telco Customer Churn Analysis

Overview

A fictional telco company in California aims to predict customer churn based on customer attributes in order to proactively retain at-risk customers through targeted retention programs. From a business perspective, the primary objective is to correctly identify customers who are highly likely to churn so that timely interventions can be made. Therefore, in building the predictive classification model, the focus is on minimizing false negatives — that is, avoiding situations where a customer is incorrectly predicted to stay when they are actually at risk of leaving. To achieve this, the model should prioritize high recall (sensitivity) for the churn class, even if this comes at the expense of slightly lower precision. This approach ensures that most potential churners are captured for retention outreach, reducing overall customer loss while allowing the business to manage campaign costs through subsequent segmentation and targeting.

In this document, we begin by conducting an exploratory data analysis to understand customer churn patterns across various demographic, behavioral, and service-related attributes. This analysis highlights which factors most strongly influence churn behavior. We then examine customer tenure and estimate the associated cost of retention to understand the financial implications of customer loss. Finally, we develop and contrast predictive classification models designed to identify high-risk customers with a focus on minimizing false negatives — ensuring that customers most likely to churn are accurately detected and prioritized for targeted retention strategies.

Key Assumptions & Findings:

• The analysis identified key churn drivers — most notably tenure, contract type, and payment method — using statistical (Cramér’s V, Point-Biserial Correlation) and machine learning (Random Forest, XGBoost) methods; by combining these insights, customers were segmented into high, medium, and low churn risk levels, enabling targeted retention strategies focused on high-risk customers to reduce overall churn and maximize lifetime value.
• Overall churn rate: There are 7,032 telco customers in California during Q3 with an overall churn rate of 26.6% (or 1,869 customers).
• Churned customers: 1,869 (26.6% of total), resulting in an estimated total financial loss of $509,282 (including $139,130 in monthly revenue loss, $56,070 in acquisition costs, and $314,082 in lifetime value loss) — approximately 55% of total customer value — underscoring the need to prioritize recall optimization to proactively identify and retain at-risk customers.
• Targeted strategy: Focus retention campaigns on 587 High Risk customers for maximum impact.

Exploratory Data Analysis

The telco customer churn data contains information about a fictional telco company that provided home phone and Internet services to 7032 customers in California in Q3. It indicates which customers have left, stayed, or signed up for their service.

The dataset is structured such that each row represents an individual customer, while each column corresponds to a specific customer attribute described in the metadata. It includes information on customers who have churned within the last month, identified in the Churn column. The dataset also details the services each customer has subscribed to, including phone, multiple lines, internet, online security, online backup, device protection, technical support, and streaming services (TV and movies). In addition, it captures account-level information such as tenure, contract type, payment method, paperless billing preference, monthly charges, and total charges. Lastly, it provides demographic details about customers, including gender, age range, and family status (partner and dependents). See Figure 1.

Figure 1: Data Definitions

Overall churn rate

We first calculate the overall churn rate is 26.6% for 7,032 customers in California during Q3 as shown in Figure 2.

Figure 2: Overall Churn Rate

Churn rate by categorical features

Next, we compute the conditional churn rate for each categorical feature, as illustrated in Figure 3. For every categorical variable, we compare the churn rates across its individual attribute levels and benchmark them against the overall churn rate estimated in the previous section. This comparison enables us to assess whether customer churn behavior varies meaningfully across categories within a feature.

To formally test the presence of such relationships, we conduct a Chi-square test of independence for each categorical feature. The null hypothesis assumes that the churn rate is independent of the categorical feature. We reject the null hypothesis when the Chi-square test statistic exceeds 10, indicating a statistically significant difference in churn rates across the feature’s attribute levels. Features meeting this criterion are therefore considered to have a statistically significant association with customer churn and are likely to contribute predictive value in subsequent modeling stages.

Figure 3 shows that attributes within gender, streaming TV, streaming movies, and phone service type exhibit minimal variation in churn rates and are statistically insignificant. In contrast, all other categorical features display meaningful and statistically significant differences — particularly payment method, contract type, and key service-related variables such as internet service, online security, online backup, device protection, and technical support — suggesting strong associations with customer churn behavior.

Figure 3: Churn Rate by Categorical Features

Cramér’s V, which measures the strength of association between categorical variables (ranging from 0 for no association to 1 for perfect association), supports these findings in Figure 3. As shown in Figure 4, contract type, internet service, and payment method exhibit the strongest associations with churn, underscoring their critical role in customer retention. In contrast, gender, streaming TV, streaming movies, and phone service type show the weakest associations, indicating limited influence on churn behavior.

Figure 4: Cramer’s V for Churn

Churn rate by numerical features

Similar to the analysis of categorical features, we examine whether numerical features exhibit distinct distributional differences between customers who churned and those who did not, as illustrated in Figure 5. This comparison helps identify whether continuous variables, such as tenure, monthly charges, and total charges, provide discriminatory power in distinguishing churn behavior.

To formally assess these differences, we apply the Mann–Whitney U test, a nonparametric alternative to the t-test suitable when data are not normally distributed. The test ranks all observations across both groups and evaluates whether their distributions are statistically different. The null hypothesis assumes identical distributions for churned and non-churned customers, while the alternative hypothesis posits a difference between them. Results show that tenure, monthly charges, and total charges are all statistically significant, with p-values < 0.05, indicating that these features have distinct distributions and carry predictive information about customer churn.

Figure 5: Churn Rate by Numerical Features

In addition to statistical significance testing, we compute Cohen’s |d| to quantify the effect size, or the magnitude of difference between churned and non-churned groups, independent of sample size. Cohen’s |d| measures the difference between two group means in units of pooled standard deviation, with thresholds of approximately 0.2 (small), 0.5 (medium), and 0.8 (large). This metric provides insight into the practical significance of observed differences beyond p-values.

Results indicate medium effect sizes for Monthly Charges (d = 0.445) and Total Charges (d = –0.461), and a large effect size for Tenure (d = –0.857), reinforcing that these numerical features exhibit meaningful and practically significant distinctions between customers who churned and those who did not.

Finally, we compute the Point-Biserial Correlation to measure the strength of the relationship between the binary churn variable and each continuous feature (ranging from –1 to +1). As shown in Figure 6, tenure shows the strongest correlation with churn, further reinforcing its predictive importance.

Figure 6: Point-Biserial Correlation for Churn

Cost of Churn

In this section, we evaluate both the immediate and long-term financial impact of customer churn. As shown in Figure 7, 1,869 customers churned during Q3, representing 26.6% of all customers in California. This translates to an estimated monthly revenue loss of $139,130 and an acquisition cost of $56,070, assuming a customer acquisition cost of $30 per customer per month based on industry benchmarks. With an approximate gross margin of 60%, the lifetime value (LTV) of these churned customers is estimated at $314,082. Thus, the short-term loss (monthly revenue plus acquisition cost) totals $195,200, while the long-term loss (LTV plus short-term loss) reaches $509,282 — equivalent to roughly 55% of the total customer value. These results emphasize the importance of optimizing recall when building the predictive churn model. In practice, this means prioritizing the accurate identification of customers likely to churn, as failing to retain them is more costly than mistakenly targeting customers who would have stayed.

Figure 7: Lifetime Value of Churned and Retained Customers

Predictive Modeling

From the feature analysis, we identified several attributes strongly associated with customer churn. While Cramér’s V and Point-Biserial Correlation quantify the strength of individual relationships between categorical and numerical features and churn, they do not capture interactions among variables. To address this and align with the business goal of maximizing recall — the proportion of actual churners correctly identified — we employ two predictive classifiers: Random Forest and XGBoost, and use them to assess feature importance.

Get Chris H’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

The Random Forest classifier constructs multiple decision trees and aggregates their predictions through majority voting. It evaluates feature relevance using Permutation Importance, which measures the drop in model performance when a feature’s values are randomly shuffled — larger declines indicate higher importance. In parallel, the XGBoost classifier, a gradient boosting algorithm that sequentially builds trees to correct residual errors from prior iterations, typically achieves higher predictive accuracy and robustness. Feature importance in XGBoost is derived using the Gain metric, which quantifies the improvement in model accuracy contributed by each feature during data splits; higher gains indicate stronger predictive influence.

We first present the confusion matrices for the Random Forest and XGBoost classifiers in Figure 8 to evaluate which model best aligns with the business objective of maximizing recall — capturing as many true churners as possible. The Random Forest classifier correctly identifies 44.0% of churned customers, while the XGBoost model achieves a substantially higher 71.5% recall. Although XGBoost exhibits a moderate trade-off in precision (i.e., a higher number of false positives of 22.1%), this is acceptable given the higher business cost of missing true churners compared to targeting some non-churners.

Figure 8: Random Forest versus XGBoost Confusion Matrix

Overall, XGBoost delivers stronger predictive performance, as reflected in its superior recall and ROC-AUC, making it the preferred model for proactive customer retention strategies as shown in Figure 9.

Figure 8: Random Forest versus XGBoost ROC Curve

Finally, Figure 9 presents the feature importance rankings for the Random Forest and XGBoost classifiers. The results reveal notable differences in how each model prioritizes features, helping explain the Random Forest’s weaker predictive performance relative to XGBoost. For example, the Random Forest model ranks gender among its top predictors, despite this feature showing a weak association with churn in Figure 4 (Cramér’s V analysis). In contrast, XGBoost assigns the highest importance to variables such as contract type, which aligns more closely with the statistically significant churn drivers identified in earlier analyses.

Figure 9: Random Forest versus XGBoost Feature Importance

Feature Importance & Customer Churn Risk Level

From the feature and predictive analyses, we identified several variables strongly associated with customer churn. The next step is to synthesize these insights by combining the information from Figures 4 (Cramér’s V), 6 (Point-Biserial Correlation), and 9 (Random Forest & XGBoost Classifier feature importance). We normalize the feature importance scores from each method and compute their average to derive a final composite feature importance score. This consolidated metric allows us to rank features consistently and identify the top 10 drivers of churn, which can then be used to flag customers at the highest risk of leaving.

Based on the average normalized importance scores, Figure 5 shows that tenure (0.897) and contract type (0.842) are the strongest predictors of churn, followed by other customer service and payment-related features.

Figure 10: Top 10 Churn Drivers by Composite Feature Importance Score

Based on Figure 10, we categorize customers into three churn risk levels, as illustrated in Figure 11. This segmentation is derived by assigning scores to feature attributes according to their associated churn rates — attributes with higher churn rates receive higher scores. For example, within the Contract feature, month-to-month, one-year, and two-year contracts are assigned scores of 3, 2, and 1 respectively. Each customer’s total churn risk score is then calculated by summing the scores across all relevant features. Customers are subsequently grouped into High, Medium, and Low Risk tiers, with those in the High Risk segment having the highest cumulative scores. Note that we convert numerical features Tenure, Total Charges, and Monthly Charges in categorical features to calculate scores.

Using this approach, we identify 587 customers who have not yet churned but fall into the High Risk category — indicating a strong likelihood of churn and highlighting a key opportunity for targeted retention strategies.

Figure 11: Customer Churn Risk Level

Conclusion

The analysis reveals several key predictors that drive customer churn. To identify the top 10 churn drivers, we integrated insights from Figure 4 (Cramér’s V), Figure 6 (Point-Biserial Correlation), and Figure 9 (Random Forest and XGBoost feature importance). By combining and normalizing these metrics, we derived a comprehensive feature importance ranking. Finally, using these top predictors, we segmented customers into high, medium, and low churn risk groups — enabling the company to prioritize retention efforts and implement targeted strategies to reduce churn.

Appendix

Libraries

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from google.colab import drive
drive.mount('/content/gdrive')

import kagglehub
from kagglehub import KaggleDatasetAdapter

import math
from scipy.stats import norm
from scipy.stats import chi2_contingency
from scipy.stats import mannwhitneyu

pd.set_option('display.max_columns', None)

Functions

def _read_data():
  """
  
  read in telco data from kaggle
  
  """
  # Set the path to the file you'd like to load
  import kagglehub

  # Download latest version
  path = kagglehub.dataset_download("blastchar/telco-customer-churn")
  file_path = path+"/WA_Fn-UseC_-Telco-Customer-Churn.csv"

  # Load the latest version
  df = pd.read_csv(file_path)

  # Recast TotalCharges
  df['TotalCharges'] = pd.to_numeric(df['TotalCharges'], errors='coerce')
  df['TotalCharges'] = df['TotalCharges'].astype('float')

  # Drop NA
  df.dropna(inplace=True)

  # Additional Transformations
  # --- SeniorCitizen Mapping
  senior_map = {0:"No", 1:"Yes"}
  df['SeniorCitizen'] = df['SeniorCitizen'].map(senior_map)

  # --- PhoneServiceType
  df["PhoneServiceType"] = (
      np.where((df["PhoneService"] == "Yes") & (df["MultipleLines"] == "Yes"),"Multiple Lines",
          np.where(df["PhoneService"] == "Yes","Phone Only",
                  "No Phone"
                )
      )
  )

  return df

def _chi_square_test(df,col,confidence = 0.95):
  contingency_table = pd.crosstab(df[col], df["Churn"])
  chi2, p_value, dof, expected = chi2_contingency(contingency_table)

  return chi2, p_value, contingency_table

def _churn_rate_by_group_calc(df,col=None,confidence=0.95):
  """

  calculate churn rate by group

  """

  if col == None:
    churn_rate = (
        df.groupby(["Churn"])
          .size()
            .reset_index(name="count")
        )

    counts = len(df)

    churn_rate["total_count"] = counts
    churn_rate["churn_rate"] = (churn_rate["count"] / churn_rate["total_count"]) * 100

    return churn_rate

  else:
    # Calculate churn count by group
    churn_rate_group = (
        df.groupby([col, "Churn"])
          .size()
            .reset_index(name="count")
        )

    # Calculate the total count for each gender
    group_counts = df.groupby(col).size().reset_index(name="total_count")

    # Merge the churn data with gender counts
    churn_rate_group = pd.merge(churn_rate_group, group_counts, on=[col])

    # Calculate the churn rate
    churn_rate_group["churn_rate"] = (churn_rate_group["count"] / churn_rate_group["total_count"] * 100)#.round(2)

    # Filter only for churn rate
    z = norm.ppf(1-(1-confidence)/2)
    churn_rate_group["error"] = z * (((churn_rate_group.churn_rate/100.0 * (1 - churn_rate_group.churn_rate/100.0))/ churn_rate_group.total_count)**(0.5))
    churn_rate_group["CI_lower"] = churn_rate_group.churn_rate - churn_rate_group.error*100.0
    churn_rate_group["CI_uppper"] = churn_rate_group.churn_rate + churn_rate_group.error*100.0

    return churn_rate_group#, chi2, p_value, contingency_table

def _risk_reduction_calc(df_churn_rate,col):
  """
  
  risk reduction calculation for additional services

  """

  # Filter for churn rate with no service and yes churn
  churn_without_service = df_churn_rate[(df_churn_rate[col] == "No") & (df_churn_rate["Churn"] == "Yes")]["churn_rate"].values[0]

  # Filter for churn rate with yes service and yes churn
  churn_with_service = df_churn_rate[(df_churn_rate[col] == "Yes") & (df_churn_rate["Churn"] == "Yes")]["churn_rate"].values[0]

  # Calculate risk reduction
  risk_reduction = (1 - churn_with_service / churn_without_service) * 100

  return risk_reduction

def _odds_calc(contingency_table, col):

  #index_values = sorted(contingency_table.index.values.tolist(), reverse=True)
  # Calculate odds ratio

  if col != "gender":
    # group and churned
    a = contingency_table.loc["Yes","Yes"]
    # group and not churned
    b = contingency_table.loc["Yes","No"]
    # not group and churned
    c = contingency_table.loc["No","Yes"]
    # not group and not churned
    d = contingency_table.loc["No","No"]

    odds_basis = "Yes"
    odds_ratio = (a * d) / (b * c)
    odds_ratio = [col,odds_basis,odds_ratio]

  else:
    # group and churned
    a = contingency_table.loc["Male","Yes"]
    # group and not churned
    b = contingency_table.loc["Male","No"]
    # not group and churned
    c = contingency_table.loc["Female","Yes"]
    # not group and not churned
    d = contingency_table.loc["Female","No"]

    odds_basis = "Male"
    odds_ratio = (a * d) / (b * c)
    odds_ratio = [col,odds_basis,odds_ratio]

  return odds_ratio

def cohens_d(group1, group2):
    """Calculate Cohen's d for effect size"""
    from math import sqrt

    # Mean difference
    mean_diff = group1.mean() - group2.mean()

    # Pooled standard deviation
    n1, n2 = len(group1), len(group2)
    var1, var2 = group1.var(), group2.var()
    pooled_std = sqrt(((n1-1)*var1 + (n2-1)*var2) / (n1 + n2 - 2))

    return mean_diff / pooled_std

def interpret_cohens_d(d):
    #if abs(d) < 0.2: return "negligible"
    #elif abs(d) < 0.5: return "small"
    #elif abs(d) < 0.8: return "medium"
    #else: return "large"
    if abs(d) < 0.2: return "small"
    elif abs(d) < 0.5: return "medium"
    elif abs(d) < 0.8: return "large"
    else: return "large"

Read data

df = _read_data()

# setup columns
cols = df.columns.tolist()
tar_cols = ["Churn"]
cus_cols = ["customerID"]
num_cols = [col for col in df.columns if (df[col].dtypes != 'O')]
rm_cat_cols = ['PhoneService','MultipleLines']
cat_cols = [col for col in df.columns if (df[col].dtypes == 'O') & (col not in tar_cols + cus_cols + rm_cat_cols)]

Exploratory Data Analysis

# unconditional churn rate

churn_rate = _churn_rate_by_group_calc(df)




# conditional churn rate - categorical variables

cat_cols_2 = ["PhoneServiceType", "PaymentMethod", "InternetService", "Contract", "gender"]

# --- storage ---
churn_rate_res = {}
chi2_res = {}
p_value_res = {}
contingency_res = {}
odds_res = {}
risk_res = {}

for c in cat_cols:
  print(c)
  if c not in cat_cols_2:
    churn_rate = _churn_rate_by_group_calc(df,col = c)
    churn_rate_res[c] = churn_rate[churn_rate[c].isin(["Yes","No"])]
    contingency_table = churn_rate_res[c][[c,"Churn","count"]].pivot(index = c, columns = "Churn", values = "count")
    contingency_res[c] = contingency_table
    chi2, p_value, dof, expected = chi2_contingency(contingency_table)
    chi2_res[c] = chi2
    p_value_res[c] = p_value
    odds_ratio = _odds_calc(contingency_table, c)
    odds_res[c] = odds_ratio
    risk = _risk_reduction_calc(churn_rate_res[c],c)
    risk_res[c] = risk
  else:
    churn_rate = _churn_rate_by_group_calc(df,col = c)
    churn_rate_res[c] = churn_rate
    contingency_table = churn_rate_res[c][[c,"Churn","count"]].pivot(index = c, columns = "Churn", values = "count")
    contingency_res[c] = contingency_table
    chi2, p_value, dof, expected = chi2_contingency(contingency_table)
    chi2_res[c] = chi2
    p_value_res[c] = p_value

# Convert scalar values to lists before creating DataFrame
chi2_res_list = {k: [v] for k, v in chi2_res.items()}
df_chi2 = pd.DataFrame(chi2_res_list).T.reset_index().rename(columns = {"index": "dimension", 0: "chi2"})

churn_rate_dim = pd.DataFrame()

for c in cat_cols:
  churn_rate_temp = churn_rate_res[c].rename(columns = {c: "dimension"})
  churn_rate_temp = churn_rate_temp[churn_rate_temp["Churn"].isin(["Yes"])].copy()
  churn_rate_temp["category"] = c
  churn_rate_dim = pd.concat([churn_rate_dim,churn_rate_temp])
  #globals()[f"churn_rate_{c}"] = churn_rate_temp[churn_rate_temp["Churn"].isin(["Yes"])]

# Get unique categories
categories = churn_rate_dim['category'].unique()

# Calculate overall churn rate for 'Yes'
overall_churn_rate_yes = df['Churn'].value_counts(normalize=True)['Yes'] * 100

# Create subplots
fig, axes = plt.subplots(3, 5, figsize=(20, 18))
axes = axes.flatten() # Flatten the 2D array of axes for easy iteration

# Create a bar plot for each category
for i, category in enumerate(categories):
    df_subset = churn_rate_dim[churn_rate_dim['category'] == category]
    # Sort by churn rate
    df_subset = df_subset.sort_values(by='churn_rate', ascending=False)

    # Get Chi-Square value for the category
    chi2_value = df_chi2[df_chi2['dimension'] == category]['chi2'].values[0]

    axes[i].bar(df_subset['dimension'], df_subset['churn_rate'], yerr=df_subset['error']*100, capsize=5)
    axes[i].set_title(f'{category} (Chi2: {chi2_value:.2f})')
    #axes[i].set_xlabel('Dimension')
    axes[i].set_ylabel('Churn Rate (%)')
    axes[i].tick_params(axis='x', rotation=45) # Rotate x-axis labels for better readability
    axes[i].grid(axis='y', linestyle='--', alpha=0.6)

    # Add horizontal line for overall churn rate
    axes[i].axhline(y=overall_churn_rate_yes, color='red', linestyle='--')

    # Annotate the overall churn rate line
    axes[i].annotate(f'Overall Churn Rate: {overall_churn_rate_yes:.2f}%',
                     xy=(axes[i].get_xlim()[1], overall_churn_rate_yes), # Position the annotation at the right edge of the subplot
                     xytext=(-10, 5), # Offset text slightly from the point
                     textcoords='offset points',
                     horizontalalignment='right', verticalalignment='bottom',
                     color='red')


# Remove any unused subplots
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

plt.tight_layout()
plt.show()





# conditional churn rate - numerical features
# --- Split data into churned and retained groups

df_churned = df[df['Churn'] == 'Yes'][num_cols]
df_retained = df[df['Churn'] == 'No'][num_cols]

# --- Plot ksdensity for churned and retained with 1 row and shading

fig, axes = plt.subplots(1, len(num_cols), figsize=(15, 5))

for i, col in enumerate(num_cols):
    df_churned[col].plot(kind='kde', ax=axes[i], label='Churned')
    df_retained[col].plot(kind='kde', ax=axes[i], label='Retained')

    # Add vertical lines for mean and median
    axes[i].axvline(df_churned[col].mean(), color='blue', linestyle='dashed', linewidth=1, label=f'Churned Mean: {df_churned[col].mean():.2f}')
    axes[i].axvline(df_churned[col].median(), color='blue', linestyle='dotted', linewidth=1, label=f'Churned Median: {df_churned[col].median():.2f}')
    axes[i].axvline(df_retained[col].mean(), color='orange', linestyle='dashed', linewidth=1, label=f'Retained Mean: {df_retained[col].mean():.2f}')
    axes[i].axvline(df_retained[col].median(), color='orange', linestyle='dotted', linewidth=1, label=f'Retained Median: {df_retained[col].median():.2f}')

    axes[i].set_title(f'KDE Plot of {col}')
    axes[i].set_xlabel(col)
    axes[i].set_ylabel('Density')
    axes[i].legend()

plt.tight_layout()
plt.show()

u_monthly, p_monthly = mannwhitneyu(
    df_churned['MonthlyCharges'],
    df_retained['MonthlyCharges'],
    alternative='two-sided'
)

# Mann-Whitney U test for TotalCharges
u_total, p_total = mannwhitneyu(
    df_churned['TotalCharges'],
    df_retained['TotalCharges'],
    alternative='two-sided'
)

# Mann-Whitney U test for tenure
u_tenure, p_tenure = mannwhitneyu(
    df_churned['tenure'],
    df_retained['tenure'],
    alternative='two-sided'
)

print("Mann-Whitney U Test Results:")
print(f"MonthlyCharges: U = {u_monthly:.0f}, p = {p_monthly:.6f}")
print(f"TotalCharges: U = {u_total:.0f}, p = {p_total:.6f}")
print(f"Tenure: U = {u_tenure:.0f}, p = {p_tenure:.6f}")

# Calculate Cohen's d for both variables
d_monthly = cohens_d(df_churned['MonthlyCharges'], df_retained['MonthlyCharges'])
d_total = cohens_d(df_churned['TotalCharges'], df_retained['TotalCharges'])
d_tenure = cohens_d(df_churned['tenure'], df_retained['tenure'])

print(f"\nCohen's d Effect Size:")

print(f"MonthlyCharges: d = {d_monthly:.3f}")
print(f"MonthlyCharges effect: {interpret_cohens_d(d_monthly)}")

print(f"TotalCharges: d = {d_total:.3f}")
print(f"TotalCharges effect: {interpret_cohens_d(d_total)}")

print(f"Tenure: d = {d_tenure:.3f}")
print(f"Tenure effect: {interpret_cohens_d(d_tenure)}")

Customer Cost Analysis

# --- customers
# ----- churned
df_churned = df[df['Churn'] == 'Yes']
n_churned = len(df_churned)
# ----- retained
df_retained = df[df['Churn'] == 'No']
n_retained = len(df_retained)

# --- monthly revenue
# ----- churned/loss
monthly_revenue_churned = df_churned['MonthlyCharges'].sum()
# ----- retained/gain
monthly_revenue_retained = df_retained['MonthlyCharges'].sum()

# --- customer acquisition cost per month per customer (industry average)
cac = 30 #300
churned_acquisition_cost = cac * n_churned
retained_acquisition_cost = cac * n_retained

# --- lifetime value
gross_margin = 0.6 # 60% margin
# ----- churned/loss
churned_monthly_charge_avg = df_churned["MonthlyCharges"].mean()
monthly_churn_rate = n_churned / (n_churned + n_retained)
ltv_loss_avg = (churned_monthly_charge_avg * gross_margin) / monthly_churn_rate
ltv_loss_total = ltv_loss_avg * n_churned
# ------ retained/gain
retained_monthly_charge_avg = df_retained["MonthlyCharges"].mean()
monthly_retain_rate = n_retained / (n_churned + n_retained)
ltv_gain_avg = (retained_monthly_charge_avg * gross_margin) / monthly_retain_rate
ltv_gain_total = ltv_gain_avg * n_retained

# --- total loss/gain (immediate + long-term)
# ----- churned/loss
total_immediate_loss = monthly_revenue_churned + churned_acquisition_cost
total_long_term_loss = ltv_loss_total + total_immediate_loss
# ----- retained/gain
total_immediate_gain = monthly_revenue_retained - retained_acquisition_cost
total_long_term_gain = ltv_gain_total + total_immediate_gain

print("CHURN LOSS ANALYSIS")
print("=" * 50)
print(f"Customers Churned: {n_churned:,}")
print(f"1. Loss Monthly Revenue: ${monthly_revenue_churned:,.2f}")
print(f"2. Acquisition Cost Wasted: ${churned_acquisition_cost:,.2f}")
print(f"3. Lifetime Value Loss: ${ltv_loss_total:,.2f}")
print("-" * 50)
print(f"TOTAL IMMEDIATE LOSS: ${total_immediate_loss:,.2f}")
print(f"TOTAL LONG-TERM LOSS: ${total_long_term_loss:,.2f}")

print("")
print("")

print("RETAINED GAIN ANALYSIS")
print("=" * 50)
print(f"Customers Retained: {n_retained:,}")
print(f"1. Gain Monthly Revenue: ${monthly_revenue_retained:,.2f}")
print(f"2. Acquisition Cost: ${retained_acquisition_cost:,.2f}")
print(f"3. Lifetime Value Gain: ${ltv_gain_total:,.2f}")
print("-" * 50)
print(f"TOTAL IMMEDIATE GAIN: ${total_immediate_gain:,.2f}")
print(f"TOTAL LONG-TERM GAIN: ${total_long_term_gain:,.2f}")

print("")
print("")

print("TOTAL NET LOSS/GAIN ANALYSIS")
print("=" * 50)
print(f"Customers Retained: {n_retained+n_churned:,}")
print(f"1. Gain Monthly Revenue: ${monthly_revenue_retained-monthly_revenue_churned:,.2f}")
print(f"2. Acquisition Cost: ${retained_acquisition_cost+churned_acquisition_cost:,.2f}")
print(f"3. Lifetime Value Gain: ${ltv_gain_total-ltv_loss_total:,.2f}")
print("-" * 50)
print(f"TOTAL IMMEDIATE GAIN: ${total_immediate_gain-total_immediate_loss:,.2f}")
print(f"TOTAL LONG-TERM GAIN: ${total_long_term_gain-total_long_term_loss:,.2f}")

#print(n_churned,monthly_revenue_churned,churned_acquisition_cost,total_immediate_loss,ltv_loss_total,total_long_term_loss)

Predictive Modeling

# Setup train and test data

# --- setup X and y
X = df.drop(tar_cols + rm_cat_cols + cus_cols + num_to_cat_cols, axis=1)
y = df[tar_cols[0]].map({'Yes': 1, 'No': 0}).astype(int)

# --- split X and y into train and test data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3, random_state = 0, stratify=y)

print(f"Shape of the feature training data {X_train.shape} and test data {X_test.shape}.")
print(f"Training set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print(f"Churn rate in training: {y_train.mean():.3f}")
print(f"Churn rate in test: {y_test.mean():.3f}")

# --- one hot encoding applied on categorical features
X_train_encoded = X_train.copy()
X_test_encoded = X_test.copy()
label_encoders = {}
for col in cat_cols:
    le = LabelEncoder()
    X_train_encoded[col] = le.fit_transform(X_train_encoded[col].astype(str))
    X_test_encoded[col] = le.fit_transform(X_test_encoded[col].astype(str))
    label_encoders[col] = le


print(f"Shape of the feature training data {X_train_encoded.shape} and test data {X_test_encoded.shape} after one hot encoding applied on categorical features.")

# random forest classifier

# --- initialize rfc model
n_estimators = 100
max_depth = 10
min_samples_leaf = 50
class_weight = "balanced"
random_state = 42
model_rfc = RandomForestClassifier(
    n_estimators = n_estimators
    #, max_depth = max_depth
    #, min_samples_leaf = min_samples_leaf
    , class_weight = class_weight
    , random_state = random_state

    )

# --- train model
model_rfc.fit(X_train_encoded, y_train)

# --- Cross-validation scores
rfc_cv_scores = cross_val_score(model_rfc, X_train_encoded, y_train, cv=5, scoring='roc_auc')

#print(f"Random Forest CV AUC: {rfc_cv_scores.mean():.3f} (+/- {rfc_cv_scores.std() * 2:.3f})")

# --- calculate prediction & metrics
y_pred_rfc = model_rfc.predict(X_test_encoded)
y_pred_proba_rfc = model_rfc.predict_proba(X_test_encoded)[:, 1]

accuracy_rfc = accuracy_score(y_test, y_pred_rfc)
precision_rfc = precision_score(y_test, y_pred_rfc)
recall_rfc = recall_score(y_test, y_pred_rfc)
f1_rfc = f1_score(y_test, y_pred_rfc)
auc_roc_rfc = roc_auc_score(y_test, y_pred_proba_rfc)

print(f"Random forest classifier performance with {n_estimators} decision-trees:")
print("="*50)
print(f"Model accuracy score : {accuracy_rfc:0.4f}")
print(f"Model precision score : {precision_rfc:0.4f}")
print(f"Model recall score : {recall_rfc:0.4f}")
print(f"Model F1 score : {f1_rfc:0.4f}")
print(f"Model AUC-ROC score : {auc_roc_rfc:0.4f}")

# --- classification report
print(f"\nClassification Report:")
print("="*50)
print(classification_report(y_test, y_pred_rfc, target_names=['Not Churn', 'Churn']))

# --- calculate ROC curves
fpr_rf, tpr_rf, _ = roc_curve(y_test, y_pred_proba_rfc)

# --- calculate AUC scores
auc_rf = roc_auc_score(y_test, y_pred_proba_rfc)

# --- calculate Precision-Recall curves
precision_rf, recall_rf, _ = precision_recall_curve(y_test, y_pred_proba_rfc)

# --- calculate AUC-PR
pr_auc_rf = auc(recall_rf, precision_rf)

# xgboost

# --- initialize rfc model
n_estimators = 100
max_depth = 10
min_samples_leaf = 50
class_weight = "balanced"
random_state = 42
learning_rate = 0.1
model_xgb = xgb.XGBClassifier(
    n_estimators = n_estimators
    #, max_depth = max_depth
    #, min_samples_leaf = min_samples_leaf
    #, class_weight = class_weight
    , learning_rate = learning_rate
    , random_state = random_state
    , scale_pos_weight=len(y_train[y_train==0]) / len(y_train[y_train==1])
    )

# --- train model
model_xgb.fit(X_train_encoded, y_train)

# --- Cross-validation scores
xgb_cv_scores = cross_val_score(model_xgb, X_train_encoded, y_train, cv=5, scoring='roc_auc')

#print(f"Random Forest CV AUC: {rfc_cv_scores.mean():.3f} (+/- {rfc_cv_scores.std() * 2:.3f})")

# --- calculate prediction & metrics
y_pred_xgb = model_xgb.predict(X_test_encoded)
y_pred_proba_xgb = model_xgb.predict_proba(X_test_encoded)[:, 1]

accuracy_xgb = accuracy_score(y_test, y_pred_xgb)
precision_xgb = precision_score(y_test, y_pred_xgb)
recall_xgb = recall_score(y_test, y_pred_xgb)
f1_xgb = f1_score(y_test, y_pred_xgb)
auc_roc_xgb = roc_auc_score(y_test, y_pred_proba_xgb)

print(f"XGboost performance with {n_estimators} decision-trees:")
print("="*50)
print(f"Model accuracy score : {accuracy_xgb:0.4f}")
print(f"Model precision score : {precision_xgb:0.4f}")
print(f"Model recall score : {recall_xgb:0.4f}")
print(f"Model F1 score : {f1_xgb:0.4f}")
print(f"Model AUC-ROC score : {auc_roc_xgb:0.4f}")

# --- classification report
print(f"\nClassification Report:")
print("="*50)
print(classification_report(y_test, y_pred_xgb, target_names=['Not Churn', 'Churn']))

# --- calculate ROC curves
fpr_xgb, tpr_xgb, _ = roc_curve(y_test, y_pred_proba_xgb)

# --- calculate AUC scores
auc_xgb = roc_auc_score(y_test, y_pred_proba_xgb)

# --- calculate Precision-Recall curves
precision_xgb, recall_xgb, _ = precision_recall_curve(y_test, y_pred_proba_xgb)

# --- calculate AUC-PR
pr_auc_xgb = auc(recall_xgb, precision_xgb)

# plot auc-roc curves

plt.figure(figsize=(6, 5))
plt.plot(fpr_rf, tpr_rf, label=f'Random Forest (AUC = {auc_rf:.3f})', linewidth=2)
plt.plot(fpr_xgb, tpr_xgb, label=f'XGBoost (AUC = {auc_xgb:.3f})', linewidth=2)
plt.plot([0, 1], [0, 1], 'k--', label='Random Classifier')

plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve - Churn Prediction Models')
plt.legend()
plt.grid(True, alpha=0.3)

# Add performance benchmarks
#for threshold in [0.3, 0.5, 0.7]:
#    idx_rf = np.argmax(y_pred_proba_rfc >= threshold)
#    if idx_rf < len(fpr_rf):
#        plt.scatter(fpr_rf[idx_rf], tpr_rf[idx_rf], s=100,
#                   label=f'RF threshold {threshold}', alpha=0.7)

plt.tight_layout()
plt.show()

def plot_confusion_matrix_heatmap(y_true, y_pred, model_name, ax):
    """Plot confusion matrix as heatmap"""
    cm = confusion_matrix(y_true, y_pred, labels=[1, 0]) # Positive (1: Churn) first, then Negative (0: Not Churn)
    cm_percentage = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis] * 100

    sns.heatmap(cm_percentage, annot=True, fmt='.1f', cmap='Blues',
                xticklabels=['Predicted Churn', 'Predicted Not Churn'],
                yticklabels=['Actual Churn', 'Actual Not Churn'],
                cbar_kws={'label': 'Percentage'}, ax=ax)

    ax.set_title(f'{model_name}\nConfusion Matrix (%)')
    ax.set_ylabel('Actual')
    ax.set_xlabel('Predicted')

# Create subplots for both models
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

plot_confusion_matrix_heatmap(y_test, y_pred_rfc, 'Random Forest', ax1)
plot_confusion_matrix_heatmap(y_test, y_pred_xgb, 'XGBoost', ax2)

plt.tight_layout()
plt.show()

# Print detailed confusion matrix analysis
def analyze_confusion_matrix(y_true, y_pred, model_name):
    cm = confusion_matrix(y_true, y_pred)
    tn, fp, fn, tp = cm.ravel()

    print(f"\n{model_name} Confusion Matrix Analysis:")
    print(f"True Negatives:  {tn} (Correctly predicted not to churn)")
    print(f"False Positives: {fp} (Predicted to churn but didn't)")
    print(f"False Negatives: {fn} (Predicted not to churn but did)")
    print(f"True Positives:  {tp} (Correctly predicted to churn)")
    print(f"False Positive Rate: {fp/(fp+tn):.3f}")
    print(f"False Negative Rate: {fn/(fn+tp):.3f}")

analyze_confusion_matrix(y_test, y_pred_rfc, "Random Forest")
analyze_confusion_matrix(y_test, y_pred_xgb, "XGBoost")

# Random Forest Feature Importance
rf_importance = pd.DataFrame({
    'feature': X.columns,
    'importance': model_rfc.feature_importances_
}).sort_values('importance', ascending=False)

# XGBoost Feature Importance
xgb_importance = pd.DataFrame({
    'feature': X.columns,
    'importance': model_xgb.feature_importances_
}).sort_values('importance', ascending=False)

# Plot feature importance comparison
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))

# Random Forest
ax1.barh(rf_importance['feature'].head(10), rf_importance['importance'].head(10))
ax1.set_title('Random Forest - Top 10 Feature Importance')
ax1.set_xlabel('Importance Score')

# XGBoost
ax2.barh(xgb_importance['feature'].head(10), xgb_importance['importance'].head(10))
ax2.set_title('XGBoost - Top 10 Feature Importance')
ax2.set_xlabel('Importance Score')

plt.tight_layout()
plt.show()

# Print top features
print("Top 10 Features - Random Forest:")
print(rf_importance.head(10).round(4))
print("\nTop 10 Features - XGBoost:")
print(xgb_importance.head(10).round(4))

Customer Churn Risk Level

# combine results
# Use pd.concat and groupby to combine the dataframes, avoiding the merge issue

# Random Forest Feature Importance
rf_importance = pd.DataFrame({
    'Feature': X.columns,
    'rfc_importance': model_rfc.feature_importances_
})

# XGBoost Feature Importance
xgb_importance = pd.DataFrame({
    'Feature': X.columns,
    'xgb_importance': model_xgb.feature_importances_
})

df_results = pd.concat([df_cramerV.set_index('Feature'),
                        df_pointbiserial_corr.set_index('Feature'),
                        rf_importance.set_index('Feature'),
                        xgb_importance.set_index('Feature')],
                          axis=1).reset_index()

# calculate normalized scores and overall rank
df_results['Cramers_V_Norm'] = normalize_scores(df_results['Cramers_V'])
df_results['PointBiserial_Corr_Norm'] = normalize_scores(df_results['PointBiserial_Corr'])
df_results['rfc_importance_Norm'] = normalize_scores(df_results['rfc_importance'])
df_results['xgb_importance_Norm'] = normalize_scores(df_results['xgb_importance'])

# calculate overall score -- average all scores

df_results['Overall_Score'] = (df_results['Cramers_V_Norm'].combine_first(df_results['PointBiserial_Corr_Norm']) +
                                df_results['rfc_importance_Norm'] +
                                df_results['xgb_importance_Norm']) / 3

#df_results[["Feature","Overall_Score","Cramers_V","PointBiserial_Corr","RF_Importance","Mutual_Info"]].sort_values(by=["Overall_Score"], ascending=[False])
#df_results[["Feature","Overall_Score","Cramers_V","PointBiserial_Corr","rfc_importance","xgb_importance"]].sort_values(by=["Overall_Score"], ascending=[False]).head(10)
df_results[["Feature","Overall_Score"]].sort_values(by=["Overall_Score"], ascending=[False]).head(10)

# Sort the dataframe by Overall_Score in descending order
df_results_sorted = df_results.sort_values(by="Overall_Score", ascending=False).head(10)

# Create a horizontal bar chart
plt.figure(figsize=(12, 4))
bars = plt.barh(df_results_sorted["Feature"], df_results_sorted["Overall_Score"], color='skyblue')
plt.xlabel("Overall Score")
plt.ylabel("Feature")
plt.title("Top 10 Churn Drivers by Composite Feature Importance Score")
plt.gca().invert_yaxis() # Invert y-axis to have the highest score at the top
plt.grid(axis='x', linestyle='--', alpha=0.6)

# Add the score values on the bars
for bar in bars:
    width = bar.get_width()
    plt.text(width, bar.get_y() + bar.get_height()/2, f'{width:.3f}', ha='left', va='center')

plt.show()

df_results_sorted = df_results[["Feature","Overall_Score"]].sort_values(by=["Overall_Score"], ascending=[False]).head(10)

churn_rate_dim = pd.DataFrame()

for c in cat_cols + num_to_cat_cols:
  churn_rate_temp = churn_rate_res[c].rename(columns = {c: "dimension"})
  churn_rate_temp = churn_rate_temp[churn_rate_temp["Churn"].isin(["Yes"])].copy()
  churn_rate_temp["Feature"] = c
  churn_rate_dim = pd.concat([churn_rate_dim,churn_rate_temp])
  #globals()[f"churn_rate_{c}"] = churn_rate_temp[churn_rate_temp["Churn"].isin(["Yes"])]

churn_rate_dim["overall_churn_rate"] = overall_churn_rate_yes
churn_rate_dim['Score'] = churn_rate_dim.groupby('Feature')['churn_rate'].rank(method='dense', ascending=True).astype(int)
churn_rate_dim = churn_rate_dim.merge(df_results_sorted[["Feature"]], how = "inner", on = ["Feature"])

churn_rate_dim

# Select relevant columns from churn_rate_dim for scoring
churn_scores = churn_rate_dim[['dimension', 'Feature', 'Score']].copy()

# Rename 'dimension' to match the column names in df where applicable
# Create a mapping for renaming columns based on the original feature name
rename_map = {
    'TenureGroup': 'TenureGroup',
    'Contract': 'Contract',
    'PaymentMethod': 'PaymentMethod',
    'TechSupport': 'TechSupport',
    'OnlineSecurity': 'OnlineSecurity',
    'TotalChargesGroup': 'TotalChargesGroup',
    'InternetService': 'InternetService',
    'MonthlyChargesGroup': 'MonthlyChargesGroup',
    'OnlineBackup': 'OnlineBackup',
    'DeviceProtection': 'DeviceProtection'
}

# Merge the original dataframe with the churn scores based on feature and dimension
df_scored = df.copy()

# Iterate through the top features and merge the scores
for feature in churn_scores['Feature'].unique():
    if feature in rename_map:
        # Select the scores for the current feature
        feature_scores = churn_scores[churn_scores['Feature'] == feature].copy()
        # Rename the 'dimension' column to the actual feature name for merging
        feature_scores.rename(columns={'dimension': rename_map[feature]}, inplace=True)
        # Merge with the main dataframe
        df_scored = pd.merge(df_scored, feature_scores[[rename_map[feature], 'Score']], on=rename_map[feature], how='left', suffixes=('', f'_{feature}'))
        # Rename the merged score column
        df_scored.rename(columns={'Score': f'Score_{feature}'}, inplace=True)


# Calculate a combined risk score for each customer
# Sum the scores from the top features (handle potential NaNs if a customer doesn't have a value in a dimension)
score_columns = [f'Score_{feature}' for feature in df_results_sorted['Feature'].head(10) if feature in rename_map]
df_scored['churn_risk_score'] = df_scored[score_columns].sum(axis=1, skipna=True)

# Define risk levels based on the combined score (these thresholds might need tuning)
def categorize_risk(score):
    if score >= df_scored['churn_risk_score'].quantile(0.8): # Top 20% as High Risk
        return '1. High Risk (Top 20%)'
    elif score >= df_scored['churn_risk_score'].quantile(0.4): # Middle 40% as Medium Risk
        return '2. Medium Risk (Middle 40%)'
    else: # Bottom 40% as Low Risk
        return '3. Low Risk (Bottom 40%)'

df_scored['risk_level'] = df_scored['churn_risk_score'].apply(categorize_risk)

# Display the first few rows with the new risk score and level
#display(df_scored[['customerID', 'churn_risk_score', 'risk_level', 'Churn']].head())

# Optional: Analyze churn rate by risk level
# Calculate normalized churn rate by risk level
#churn_by_risk_level_normalized = df_scored.groupby('risk_level')['Churn'].value_counts(normalize=True).unstack() * 100

# Calculate raw counts by risk level
churn_by_risk_level_counts = df_scored.groupby('risk_level')['Churn'].value_counts().unstack()

# Combine normalized rates and raw counts
#churn_by_risk_level_combined = pd.concat([churn_by_risk_level_counts.add_suffix(' (Count)'), churn_by_risk_level_normalized.add_suffix(' (%)')], axis=1)

display(churn_by_risk_level_counts)

# should be concerned with potential churn when Churn = "No" and risk level 1