Telcom Churn Prediction¶

# Download the dataset manually if necessary and place it in the 'data/' directory
# Then, load the dataset

import pandas as pd

# Load data
df = pd.read_csv("https://drive.google.com/uc?export=download&id=1A3MUldrs0z09DlYR6Y1utfySwKNO9Qsz", index_col=0)

# Display first few rows
df.head()

# Check for missing values and data types
print(df.info())

plt.figure(figsize=(17, 9))
plt.title('Comparison between tenure and Monthly Charges (Churn=color)')
sns.scatterplot(x=df["tenure"], y=df['MonthlyCharges'], hue=df['Churn'], s=50, alpha=0.4);

# create tenure_new by adding randomly +/- 0.4
df2 = df.copy()
df2['tenure_new'] = df2['tenure'] + np.random.uniform(-0.4, 0.4, size=len(df2))
plt.figure(figsize=(17, 9))
plt.title('Comparison between tenure and Monthly Charges (Churn=color)')
sns.scatterplot(x=df2["tenure_new"], y=df2['MonthlyCharges'], hue=df2['Churn'], s=50, alpha=0.4);

# Use one-hot encoding for categorical features
df_encoded = pd.get_dummies(df.drop('customerID', axis=1), drop_first=True)


df_encoded.drop(columns=['MonthlyCharges', 'TotalCharges', 'tenure']).mean().plot(
    kind='bar', figsize=(12, 6), title='Feature Distribution (percentage)'
)

print("Mean values of encoded features:")
print(df_encoded.mean())

import numpy as np
# Feature engineering

# Convert 'SeniorCitizen' to boolean
df_encoded["SeniorCitizen"] = df_encoded["SeniorCitizen"].astype(bool)

# Avg_MonthlyCharges: handle division by zero
df_encoded["Avg_MonthlyCharges"] = df_encoded["TotalCharges"] / df_encoded["tenure"].replace(0, np.nan)

# Count of active services (engagement depth)
df_encoded["Services_Count"] = df_encoded[
    [
        "PhoneService_Yes",
        "MultipleLines_Yes",
        "InternetService_Fiber optic",
        "OnlineSecurity_Yes",
        "OnlineBackup_Yes",
        "DeviceProtection_Yes",
        "TechSupport_Yes",
        "StreamingTV_Yes",
        "StreamingMovies_Yes",
    ]
].sum(axis=1)

# Protection bundle score
df_encoded["TechServices_Used"] = df_encoded[["OnlineSecurity_Yes", "OnlineBackup_Yes", "DeviceProtection_Yes", "TechSupport_Yes"]].sum(
    axis=1
)

# Heavy streamer
df_encoded["Is_Heavy_Streamer"] = (df_encoded["StreamingTV_Yes"] == 1) & (df_encoded["StreamingMovies_Yes"] == 1)

# Has multiple services (3+)
df_encoded["Has_MultipleServices"] = df_encoded["Services_Count"] >= 3

# Has tech support only
df_encoded["Has_TechSupport_Only"] = (
    (df_encoded["TechSupport_Yes"] == 1)
    & (df_encoded["OnlineSecurity_Yes"] == 0)
    & (df_encoded["OnlineBackup_Yes"] == 0)
    & (df_encoded["DeviceProtection_Yes"] == 0)
)

# What the customer would have paid if charged current monthly rate all along
df_encoded["UniformRateTotal"] = df_encoded["MonthlyCharges"] * df_encoded["tenure"]

# Difference between uniform-rate total and actual historical payments
df_encoded["HistoricalRateGap"] = df_encoded["UniformRateTotal"] - df_encoded["TotalCharges"]

# Historical average monthly rate
df_encoded["HistoricalAvgMonthlyRate"] = df_encoded["TotalCharges"] / df_encoded["tenure"]

# Difference between current monthly rate and historical average monthly rate
df_encoded["CurrentChargeOverHistorical"] = df_encoded["MonthlyCharges"] - df_encoded["HistoricalAvgMonthlyRate"]


# Paperless + autopay = loyal?
df_encoded["Is_Paperless_AutoPay"] = (df_encoded["PaperlessBilling_Yes"] == 1) & (df_encoded["PaymentMethod_Credit card (automatic)"] == 1)

# High-risk payment group
df_encoded["High_Risk_Payment"] = (df_encoded["PaperlessBilling_Yes"] == 1) & (df_encoded["PaymentMethod_Electronic check"] == 1)

# Social support features
df_encoded["Is_SeniorAlone"] = (df_encoded["SeniorCitizen"] == 1) & (df_encoded["Partner_Yes"] == 0) & (df_encoded["Dependents_Yes"] == 0)

df_encoded["Is_Young_Family"] = (df_encoded["SeniorCitizen"] == 0) & (
    (df_encoded["Partner_Yes"] == 1) | (df_encoded["Dependents_Yes"] == 1)
)

# Categoricl features for tenure and monthly charges
df_encoded["Is_NewCostumer"] = df_encoded["tenure"] == 1
df_encoded["Is_Longterm"] = df_encoded["tenure"] > 20
df_encoded["Is_MaximumTenure"] = df_encoded["tenure"] == 72
df_encoded["CheapFlatrate_yes"] = df_encoded["MonthlyCharges"] < 27
df_encoded["HighCost_g60"] = df_encoded["MonthlyCharges"] > 60

# Categorical features for monthly charges


# Hardcoded tariff bin boundaries based on  clustering
boundaries = [-np.inf, 25.53625, 34.4975, 48.31625, 66.9925, 88.14775, np.inf]
labels = ["VeryLow", "Low", "Medium", "MediumHigh", "High", "Premium"]

# Assign tariff band
df_encoded["TariffBand"] = pd.cut(df_encoded["MonthlyCharges"], bins=boundaries, labels=labels)

# One-hot encode as boolean features
df_encoded = pd.get_dummies(df_encoded, columns=["TariffBand"], prefix="Tariff", dtype=bool)


df_encoded.mean()

df_encoded.info()

df_encoded.describe()

# show the distribution of the numerical features
df_encoded.select_dtypes(include='number').hist(figsize=(15, 12), bins=60, grid=True)

from catboost import CatBoostClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import confusion_matrix, classification_report, precision_recall_curve, ConfusionMatrixDisplay
from sklearn.model_selection import StratifiedKFold
from imblearn.over_sampling import RandomOverSampler
import shap, pandas as pd, numpy as np, matplotlib.pyplot as plt
from IPython.display import display

# 📌 Oversample original training set
print("📌 Applying oversampling to balance full training set...")
ros = RandomOverSampler(random_state=42)
X_resampled, y_resampled = ros.fit_resample(X_train, y_train)
print("➡️ Class distribution after resampling:\n", pd.Series(y_resampled).value_counts())

# 📌 Standardize data for Logistic Regression
print("📌 Scaling features...")
scaler = StandardScaler()
X_resampled_scaled = scaler.fit_transform(X_resampled)
X_test_scaled = scaler.transform(X_test)

# 📊 Cross-validation setup
print("📌 Performing Stratified 5-Fold CV...")
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# ⏱️ Store results
cat_preds, lr_preds = [], []
cat_probs, lr_probs = [], []
cat_true, lr_true = [], []

fold = 0
for train_idx, val_idx in skf.split(X_resampled, y_resampled):
    fold += 1
    print(f"\n🔁 Fold {fold}")

    X_tr, X_val = X_resampled.iloc[train_idx], X_resampled.iloc[val_idx]
    y_tr, y_val = y_resampled.iloc[train_idx], y_resampled.iloc[val_idx]

    X_tr_scaled = scaler.fit_transform(X_tr)
    X_val_scaled = scaler.transform(X_val)

    # Train CatBoost
    cat = CatBoostClassifier(verbose=0, random_state=42, iterations=600, learning_rate=0.05)
    cat.fit(X_tr, y_tr)
    cat_preds.append(cat.predict(X_val))
    cat_probs.append(cat.predict_proba(X_val)[:, 1])
    cat_true.append(y_val)

    # Train Logistic Regression
    lr = LogisticRegression(max_iter=5000, random_state=42)
    lr.fit(X_tr_scaled, y_tr)
    lr_preds.append(lr.predict(X_val_scaled))
    lr_probs.append(lr.predict_proba(X_val_scaled)[:, 1])
    lr_true.append(y_val)

# Concatenate results
y_true_cat = pd.concat(cat_true)
y_pred_cat = np.concatenate(cat_preds)
y_prob_cat = np.concatenate(cat_probs)

y_true_lr = pd.concat(lr_true)
y_pred_lr = np.concatenate(lr_preds)
y_prob_lr = np.concatenate(lr_probs)

# 📋 Classification Reports
print("\n📋 Classification Report – CatBoost:")
print(classification_report(y_true_cat, y_pred_cat, target_names=["No", "Yes"]))

print("\n📋 Classification Report – Logistic Regression:")
print(classification_report(y_true_lr, y_pred_lr, target_names=["No", "Yes"]))

# 📊 Confusion Matrices
print("📊 Confusion Matrices...")
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
ConfusionMatrixDisplay.from_predictions(y_true_cat, y_pred_cat, ax=ax[0], display_labels=["No", "Yes"])
ax[0].set_title("CatBoost Confusion Matrix")
ConfusionMatrixDisplay.from_predictions(y_true_lr, y_pred_lr, ax=ax[1], display_labels=["No", "Yes"])
ax[1].set_title("LogReg Confusion Matrix")
plt.tight_layout()
plt.show()

# 🧠 SHAP values on sample from full training set
print("📊 Computing SHAP values on 1000-sample from full training set...")
X_sample_cat = X_resampled.sample(1000, random_state=42)
X_sample_lr = pd.DataFrame(scaler.transform(X_sample_cat), columns=X_sample_cat.columns, index=X_sample_cat.index)

explainer_cat = shap.Explainer(cat)
shap_vals_cat = explainer_cat(X_sample_cat)

# ✅ Use modern masker API to avoid warning
masker = shap.maskers.Independent(X_resampled_scaled)
explainer_lr = shap.LinearExplainer(lr, masker)
shap_vals_lr = explainer_lr(X_sample_lr)

# 🔍 SHAP Summary Plots
shap.summary_plot(shap_vals_cat, X_sample_cat, max_display=15, show=False, alpha=0.2)
plt.title("CatBoost – SHAP Feature Importance (Top 15)")
plt.tight_layout()
plt.show()

shap.summary_plot(shap_vals_lr, X_sample_cat, max_display=15, show=False, alpha=0.2)
plt.title("LogReg – SHAP Feature Importance (Top 15)")
plt.tight_layout()
plt.show()

# 📊 Built-in feature importances
print("📊 Gathering feature importances...")
cat_imp = pd.Series(cat.feature_importances_, index=X_train.columns)
lr_imp = pd.Series(np.abs(lr.coef_[0]), index=X_train.columns)

shap_mean_cat = pd.Series(np.abs(shap_vals_cat.values).mean(axis=0), index=X_sample_cat.columns)
shap_mean_lr = pd.Series(np.abs(shap_vals_lr.values).mean(axis=0), index=X_sample_cat.columns)

all_importances = pd.DataFrame(
    {"CatBoost_SHAP": shap_mean_cat, "LogReg_SHAP": shap_mean_lr, "CatBoost_Builtin": cat_imp, "LogReg_Builtin": lr_imp}
).sort_values("CatBoost_SHAP", ascending=False)

print("📋 Combined Feature Importances:")
display(all_importances)


# 🎯 Precision/Recall Curve
def plot_precision_recall_vs_threshold(y_true, y_scores, label):
    precisions, recalls, thresholds = precision_recall_curve(y_true, y_scores)
    plt.plot(thresholds, precisions[:-1], label=f'{label} – Precision', linestyle='-')
    plt.plot(thresholds, recalls[:-1], label=f'{label} – Recall', linestyle='--')


print("📊 Precision and Recall vs Threshold...")
plt.figure(figsize=(10, 6))
plot_precision_recall_vs_threshold(y_true_cat, y_prob_cat, label="CatBoost")
plot_precision_recall_vs_threshold(y_true_lr, y_prob_lr, label="LogReg")
plt.xlabel("Threshold")
plt.ylabel("Score")
plt.title("Precision and Recall vs. Threshold")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

from catboost import CatBoostClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import confusion_matrix, classification_report, precision_recall_curve, ConfusionMatrixDisplay
from sklearn.model_selection import StratifiedKFold
from imblearn.over_sampling import RandomOverSampler
import shap, pandas as pd, numpy as np, matplotlib.pyplot as plt
from IPython.display import display

# 📌 Setup CV
print("📌 Setting up Stratified 5-Fold CV...")
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# ⏱️ Store results
cat_preds, lr_preds = [], []
cat_probs, lr_probs = [], []
cat_true, lr_true = [], []

# 🧪 CV loop with oversampling inside
fold = 0
scaler = StandardScaler()  # Create scaler once
for train_idx, val_idx in skf.split(X_train, y_train):
    fold += 1
    print(f"\n🔁 Fold {fold}")

    # Split and oversample
    X_tr, y_tr = X_train.iloc[train_idx], y_train.iloc[train_idx]
    X_val, y_val = X_train.iloc[val_idx], y_train.iloc[val_idx]
    ros = RandomOverSampler(random_state=42)
    X_tr_resampled, y_tr_resampled = ros.fit_resample(X_tr, y_tr)

    # Scale for LogReg
    X_tr_scaled = scaler.fit_transform(X_tr_resampled)
    X_val_scaled = scaler.transform(X_val)

    # Train CatBoost
    cat = CatBoostClassifier(verbose=0, random_state=42, iterations=600, learning_rate=0.05)
    cat.fit(X_tr_resampled, y_tr_resampled)
    cat_preds.append(cat.predict(X_val))
    cat_probs.append(cat.predict_proba(X_val)[:, 1])
    cat_true.append(y_val)

    # Train Logistic Regression
    lr = LogisticRegression(max_iter=5000, random_state=42)
    lr.fit(X_tr_scaled, y_tr_resampled)
    lr_preds.append(lr.predict(X_val_scaled))
    lr_probs.append(lr.predict_proba(X_val_scaled)[:, 1])
    lr_true.append(y_val)

# 📊 Concatenate results
y_true_cat = pd.concat(cat_true)
y_pred_cat = np.concatenate(cat_preds)
y_prob_cat = np.concatenate(cat_probs)

y_true_lr = pd.concat(lr_true)
y_pred_lr = np.concatenate(lr_preds)
y_prob_lr = np.concatenate(lr_probs)

# 📋 Classification Reports
print("\n📋 Classification Report – CatBoost:")
print(classification_report(y_true_cat, y_pred_cat, target_names=["No", "Yes"]))

print("\n📋 Classification Report – Logistic Regression:")
print(classification_report(y_true_lr, y_pred_lr, target_names=["No", "Yes"]))

# 📊 Confusion Matrices
print("📊 Confusion Matrices...")
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
ConfusionMatrixDisplay.from_predictions(y_true_cat, y_pred_cat, ax=ax[0], display_labels=["No", "Yes"])
ax[0].set_title("CatBoost Confusion Matrix")
ConfusionMatrixDisplay.from_predictions(y_true_lr, y_pred_lr, ax=ax[1], display_labels=["No", "Yes"])
ax[1].set_title("LogReg Confusion Matrix")
plt.tight_layout()
plt.show()

# 🔁 Re-train final models on full resampled training set for SHAP
print("📌 Retraining models on full balanced training set for SHAP...")
ros = RandomOverSampler(random_state=42)
X_resampled, y_resampled = ros.fit_resample(X_train, y_train)
X_resampled_scaled = scaler.fit_transform(X_resampled)

cat_final = CatBoostClassifier(verbose=0, random_state=42, iterations=600, learning_rate=0.05)
cat_final.fit(X_resampled, y_resampled)

lr_final = LogisticRegression(max_iter=5000, random_state=42)
lr_final.fit(X_resampled_scaled, y_resampled)

# 📊 SHAP values on sample
print("📊 Computing SHAP values on 1000-sample from full resampled set...")
X_sample_cat = X_resampled.sample(1000, random_state=42)
X_sample_scaled = pd.DataFrame(scaler.transform(X_sample_cat), columns=X_sample_cat.columns, index=X_sample_cat.index)

explainer_cat = shap.Explainer(cat_final)
shap_vals_cat = explainer_cat(X_sample_cat)

explainer_lr = shap.Explainer(lr_final, X_sample_scaled)
shap_vals_lr = explainer_lr(X_sample_scaled)

# 🔍 SHAP Summary Plots
shap.summary_plot(shap_vals_cat, X_sample_cat, max_display=15, show=False, alpha=0.2)
plt.title("CatBoost – SHAP Feature Importance (Top 15)")
plt.tight_layout()
plt.show()

shap.summary_plot(shap_vals_lr, X_sample_cat, max_display=15, show=False, alpha=0.2)
plt.title("LogReg – SHAP Feature Importance (Top 15)")
plt.tight_layout()
plt.show()

# 📊 Built-in feature importances
print("📊 Gathering feature importances...")
cat_importance = pd.Series(cat_final.feature_importances_, index=X_train.columns)
lr_coeff = pd.Series(np.abs(lr_final.coef_[0]), index=X_train.columns)

shap_mean_cat = pd.Series(np.abs(shap_vals_cat.values).mean(axis=0), index=X_sample_cat.columns)
shap_mean_lr = pd.Series(np.abs(shap_vals_lr.values).mean(axis=0), index=X_sample_cat.columns)

# 🧾 Combine into one DataFrame
all_importances = pd.DataFrame(
    {"CatBoost_SHAP": shap_mean_cat, "LogReg_SHAP": shap_mean_lr, "CatBoost_Builtin": cat_importance, "LogReg_Coefficient": lr_coeff}
).sort_values("CatBoost_SHAP", ascending=False)

print("📋 Combined Feature Importances:")
display(all_importances)


# 🎯 Precision/Recall Curve
def plot_precision_recall_vs_threshold(y_true, y_scores, label):
    precisions, recalls, thresholds = precision_recall_curve(y_true, y_scores)
    plt.plot(thresholds, precisions[:-1], label=f'{label} – Precision', linestyle='-')
    plt.plot(thresholds, recalls[:-1], label=f'{label} – Recall', linestyle='--')


print("📊 Precision and Recall vs Threshold...")
plt.figure(figsize=(10, 6))
plot_precision_recall_vs_threshold(y_true_cat, y_prob_cat, label="CatBoost")
plot_precision_recall_vs_threshold(y_true_lr, y_prob_lr, label="LogReg")
plt.xlabel("Threshold")
plt.ylabel("Score")
plt.title("Precision and Recall vs. Threshold")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# Example: Evaluate CatBoost on real test set
print("📊 Final Evaluation on X_test:")
cat_test_pred = cat_final.predict(X_test)
lr_test_pred = lr_final.predict(scaler.transform(X_test))

print("CatBoost:")
print(classification_report(y_test, cat_test_pred))

print("LogReg:")
print(classification_report(y_test, lr_test_pred))