import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split, StratifiedKFold, GridSearchCV
from sklearn.preprocessing import StandardScaler, OneHotEncoder, LabelEncoder
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import (accuracy_score, f1_score, roc_auc_score,
                             classification_report, confusion_matrix, ConfusionMatrixDisplay,
                             RocCurveDisplay)
import warnings
warnings.filterwarnings('ignore')

# Plotting style
sns.set_theme(style='whitegrid', palette='muted', font_scale=1.1)
plt.rcParams['figure.dpi'] = 110

url = "https://raw.githubusercontent.com/IBM/telco-customer-churn-on-icp4d/master/data/Telco-Customer-Churn.csv"
df = pd.read_csv(url)

print(f"Shape: {df.shape}")
print(f"\nColumns:\n{list(df.columns)}")
df.head(3)

# --- 3.1 Basic inspection ---
print("Dtypes and nulls:")
print(df.dtypes)
print(f"\nMissing values per column:\n{df.isnull().sum()}")

# --- 3.2 Cleaning ---

# TotalCharges is stored as object due to whitespace entries; fix it
df['TotalCharges'] = pd.to_numeric(df['TotalCharges'], errors='coerce')

# Drop the 11 rows where TotalCharges couldn't be parsed (new customers, tenure=0)
df.dropna(subset=['TotalCharges'], inplace=True)

# Drop customerID — it's a unique identifier with no predictive value
df.drop(columns=['customerID'], inplace=True)

# SeniorCitizen is already 0/1 — leave it as numeric
# Encode target to binary int
df['Churn'] = (df['Churn'] == 'Yes').astype(int)

print(f"Clean shape: {df.shape}")
print(f"\nChurn distribution:\n{df['Churn'].value_counts()}")
print(f"Churn rate: {df['Churn'].mean():.1%}")

# --- 3.3 Define feature types ---

target = 'Churn'

numeric_features = ['tenure', 'MonthlyCharges', 'TotalCharges']

categorical_features = [
    'gender', 'SeniorCitizen', 'Partner', 'Dependents',
    'PhoneService', 'MultipleLines', 'InternetService',
    'OnlineSecurity', 'OnlineBackup', 'DeviceProtection',
    'TechSupport', 'StreamingTV', 'StreamingMovies',
    'Contract', 'PaperlessBilling', 'PaymentMethod'
]

X = df[numeric_features + categorical_features]
y = df[target]

print(f"Features: {X.shape[1]} | Target classes: {y.value_counts().to_dict()}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

# Countplot
churn_counts = y.value_counts()
axes[0].bar(['Stayed', 'Churned'], churn_counts.values, color=['#4C72B0', '#DD8452'], edgecolor='white', linewidth=1.2)
axes[0].set_title('Class Balance: Churn vs. Stayed', fontweight='bold')
axes[0].set_ylabel('Customer Count')
for i, v in enumerate(churn_counts.values):
    axes[0].text(i, v + 50, f'{v:,}\n({v/len(y):.0%})', ha='center', fontsize=10)

# Monthly charges by churn
df_plot = df.copy()
df_plot['Churn_Label'] = df_plot['Churn'].map({0: 'Stayed', 1: 'Churned'})
axes[1].hist(df_plot[df_plot['Churn']==0]['MonthlyCharges'], bins=30, alpha=0.7, label='Stayed', color='#4C72B0')
axes[1].hist(df_plot[df_plot['Churn']==1]['MonthlyCharges'], bins=30, alpha=0.7, label='Churned', color='#DD8452')
axes[1].set_title('Monthly Charges Distribution by Churn', fontweight='bold')
axes[1].set_xlabel('Monthly Charges ($)')
axes[1].set_ylabel('Count')
axes[1].legend()

plt.tight_layout()
plt.savefig('viz1_class_balance.png', bbox_inches='tight')
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Churn rate by contract
churn_by_contract = df.groupby('Contract')['Churn'].mean().sort_values(ascending=False)
axes[0].bar(churn_by_contract.index, churn_by_contract.values * 100,
            color=['#DD8452', '#4C72B0', '#55A868'], edgecolor='white', linewidth=1.2)
axes[0].set_title('Churn Rate by Contract Type', fontweight='bold')
axes[0].set_ylabel('Churn Rate (%)')
for i, v in enumerate(churn_by_contract.values):
    axes[0].text(i, v * 100 + 0.5, f'{v:.0%}', ha='center', fontsize=10, fontweight='bold')

# Tenure vs Churn (box)
df_plot_tenure = df.copy()
df_plot_tenure['Churn_Label'] = df_plot_tenure['Churn'].map({0: 'Stayed', 1: 'Churned'})
axes[1].boxplot(
    [df_plot_tenure[df_plot_tenure['Churn']==0]['tenure'],
     df_plot_tenure[df_plot_tenure['Churn']==1]['tenure']],
    labels=['Stayed', 'Churned'],
    patch_artist=True,
    boxprops=dict(facecolor='#4C72B0', alpha=0.7),
    medianprops=dict(color='black', linewidth=2)
)
axes[1].set_title('Tenure Distribution by Churn', fontweight='bold')
axes[1].set_ylabel('Tenure (months)')

plt.tight_layout()
plt.savefig('viz2_contract_tenure.png', bbox_inches='tight')
plt.show()

# --- Stratified 80/20 split ---
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.20, random_state=42, stratify=y
)

print(f"Train size : {X_train.shape[0]:,} | Churn rate: {y_train.mean():.1%}")
print(f"Test  size : {X_test.shape[0]:,}  | Churn rate: {y_test.mean():.1%}")

# --- Cross-validation strategy ---
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
# StratifiedKFold preserves churn rate in every fold — essential for imbalanced data.

# Build a reusable ColumnTransformer
preprocessor = ColumnTransformer(transformers=[
    ('num', StandardScaler(), numeric_features),
    ('cat', OneHotEncoder(handle_unknown='ignore', sparse_output=False), categorical_features)
])

# Fit on train, transform both (no leakage)
X_train_proc = preprocessor.fit_transform(X_train)
X_test_proc  = preprocessor.transform(X_test)

print(f"Processed feature matrix shape (train): {X_train_proc.shape}")

# Helper function for clean metric reporting
def evaluate(name, model, X_tr, y_tr, X_te, y_te):
    """Fit model, return dict with key metrics."""
    model.fit(X_tr, y_tr)
    y_pred = model.predict(X_te)
    y_prob = model.predict_proba(X_te)[:, 1] if hasattr(model, 'predict_proba') else None
    
    acc   = accuracy_score(y_te, y_pred)
    f1    = f1_score(y_te, y_pred, average='weighted')
    auc   = roc_auc_score(y_te, y_prob) if y_prob is not None else np.nan
    
    print(f"\n{'='*40}")
    print(f"  {name}")
    print(f"{'='*40}")
    print(f"  Accuracy : {acc:.4f}")
    print(f"  F1 (wtd) : {f1:.4f}")
    print(f"  ROC AUC  : {auc:.4f}")
    
    return {'Model': name, 'Accuracy': acc, 'F1_weighted': f1, 'ROC_AUC': auc,
            'Best_Params': 'defaults', 'Notes': ''}

baseline_lr = LogisticRegression(max_iter=1000, class_weight='balanced', random_state=42)
baseline_results = evaluate(
    'Logistic Regression (Baseline)', baseline_lr,
    X_train_proc, y_train, X_test_proc, y_test
)

knn_pipe = Pipeline([('knn', KNeighborsClassifier())])

knn_grid = {
    'knn__n_neighbors': [3, 5, 7, 11, 15, 21],
    'knn__weights': ['uniform', 'distance'],
    'knn__p': [1, 2]
}

knn_search = GridSearchCV(
    KNeighborsClassifier(),
    {k.replace('knn__', ''): v for k, v in knn_grid.items()},
    cv=cv, scoring='roc_auc', n_jobs=-1, verbose=0
)

knn_search.fit(X_train_proc, y_train)

print(f"Best KNN params : {knn_search.best_params_}")
print(f"Best CV AUC     : {knn_search.best_score_:.4f}")

# Evaluate best KNN on test set
knn_best = knn_search.best_estimator_
y_pred_knn = knn_best.predict(X_test_proc)
y_prob_knn = knn_best.predict_proba(X_test_proc)[:, 1]

knn_results = {
    'Model': 'KNN',
    'Accuracy': accuracy_score(y_test, y_pred_knn),
    'F1_weighted': f1_score(y_test, y_pred_knn, average='weighted'),
    'ROC_AUC': roc_auc_score(y_test, y_prob_knn),
    'Best_Params': str(knn_search.best_params_),
    'Notes': 'Fast inference; sensitive to scale (handled)'
}

print(f"\nKNN Test Accuracy : {knn_results['Accuracy']:.4f}")
print(f"KNN Test F1       : {knn_results['F1_weighted']:.4f}")
print(f"KNN Test ROC AUC  : {knn_results['ROC_AUC']:.4f}")

# Diagnostic: Confusion Matrix for KNN
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
ConfusionMatrixDisplay.from_predictions(
    y_test, y_pred_knn, display_labels=['Stayed', 'Churned'],
    colorbar=False, ax=ax, cmap='Blues'
)
ax.set_title('KNN — Confusion Matrix (Test Set)', fontweight='bold')
plt.tight_layout()
plt.savefig('cm_knn.png', bbox_inches='tight')
plt.show()

dt_grid = {
    'max_depth': [3, 5, 7, 10, None],
    'min_samples_split': [10, 20, 50],
    'criterion': ['gini', 'entropy']
}

dt_search = GridSearchCV(
    DecisionTreeClassifier(class_weight='balanced', random_state=42),
    dt_grid, cv=cv, scoring='roc_auc', n_jobs=-1
)

dt_search.fit(X_train_proc, y_train)

print(f"Best DT params : {dt_search.best_params_}")
print(f"Best CV AUC    : {dt_search.best_score_:.4f}")

dt_best = dt_search.best_estimator_
y_pred_dt = dt_best.predict(X_test_proc)
y_prob_dt = dt_best.predict_proba(X_test_proc)[:, 1]

dt_results = {
    'Model': 'Decision Tree',
    'Accuracy': accuracy_score(y_test, y_pred_dt),
    'F1_weighted': f1_score(y_test, y_pred_dt, average='weighted'),
    'ROC_AUC': roc_auc_score(y_test, y_prob_dt),
    'Best_Params': str(dt_search.best_params_),
    'Notes': 'Interpretable; prone to overfit without depth cap'
}

print(f"\nDecision Tree Test Accuracy : {dt_results['Accuracy']:.4f}")
print(f"Decision Tree Test F1       : {dt_results['F1_weighted']:.4f}")
print(f"Decision Tree Test ROC AUC  : {dt_results['ROC_AUC']:.4f}")

# Diagnostic: Confusion Matrix for Decision Tree
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
ConfusionMatrixDisplay.from_predictions(
    y_test, y_pred_dt, display_labels=['Stayed', 'Churned'],
    colorbar=False, ax=ax, cmap='Oranges'
)
ax.set_title('Decision Tree — Confusion Matrix (Test Set)', fontweight='bold')
plt.tight_layout()
plt.savefig('cm_dt.png', bbox_inches='tight')
plt.show()

rf_grid = {
    'n_estimators': [100, 200],
    'max_depth': [5, 10, None],
    'max_features': ['sqrt', 'log2'],
    'min_samples_leaf': [1, 5]
}

rf_search = GridSearchCV(
    RandomForestClassifier(class_weight='balanced', random_state=42),
    rf_grid, cv=cv, scoring='roc_auc', n_jobs=-1
)

rf_search.fit(X_train_proc, y_train)

print(f"Best RF params : {rf_search.best_params_}")
print(f"Best CV AUC    : {rf_search.best_score_:.4f}")

rf_best = rf_search.best_estimator_
y_pred_rf = rf_best.predict(X_test_proc)
y_prob_rf = rf_best.predict_proba(X_test_proc)[:, 1]

rf_results = {
    'Model': 'Random Forest',
    'Accuracy': accuracy_score(y_test, y_pred_rf),
    'F1_weighted': f1_score(y_test, y_pred_rf, average='weighted'),
    'ROC_AUC': roc_auc_score(y_test, y_prob_rf),
    'Best_Params': str(rf_search.best_params_),
    'Notes': 'Best AUC; feature importances actionable for business'
}

print(f"\nRandom Forest Test Accuracy : {rf_results['Accuracy']:.4f}")
print(f"Random Forest Test F1       : {rf_results['F1_weighted']:.4f}")
print(f"Random Forest Test ROC AUC  : {rf_results['ROC_AUC']:.4f}")

# Diagnostic 1: Confusion Matrix
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
ConfusionMatrixDisplay.from_predictions(
    y_test, y_pred_rf, display_labels=['Stayed', 'Churned'],
    colorbar=False, ax=ax, cmap='Greens'
)
ax.set_title('Random Forest — Confusion Matrix (Test Set)', fontweight='bold')
plt.tight_layout()
plt.savefig('cm_rf.png', bbox_inches='tight')
plt.show()

# Diagnostic 2: Feature Importances (top 15)
feature_names = (
    numeric_features +
    list(preprocessor.named_transformers_['cat'].get_feature_names_out(categorical_features))
)

importances = pd.Series(rf_best.feature_importances_, index=feature_names)
top15 = importances.nlargest(15).sort_values()

fig, ax = plt.subplots(figsize=(8, 5))
top15.plot(kind='barh', ax=ax, color='#4C72B0', edgecolor='white')
ax.set_title('Random Forest — Top 15 Feature Importances', fontweight='bold')
ax.set_xlabel('Importance Score')
plt.tight_layout()
plt.savefig('feature_importance_rf.png', bbox_inches='tight')
plt.show()

results_list = [baseline_results, knn_results, dt_results, rf_results]
results_df = pd.DataFrame(results_list)[[
    'Model', 'Best_Params', 'Accuracy', 'F1_weighted', 'ROC_AUC', 'Notes'
]]

# Format numeric columns
for col in ['Accuracy', 'F1_weighted', 'ROC_AUC']:
    results_df[col] = results_df[col].map('{:.4f}'.format)

results_df.columns = ['Model', 'Best Hyperparameters', 'Accuracy', 'F1 (Weighted)', 'ROC AUC', 'Notes']
print(results_df.to_string(index=False))

# ROC Curves — all models on one plot
fig, ax = plt.subplots(figsize=(8, 6))

for name, prob in [
    ('Logistic Regression (Baseline)', baseline_lr.predict_proba(X_test_proc)[:, 1]),
    ('KNN', y_prob_knn),
    ('Decision Tree', y_prob_dt),
    ('Random Forest', y_prob_rf),
]:
    RocCurveDisplay.from_predictions(y_test, prob, name=name, ax=ax)

ax.plot([0, 1], [0, 1], 'k--', lw=1, label='Random Classifier')
ax.set_title('ROC Curves — All Models (Test Set)', fontweight='bold')
ax.legend(loc='lower right', fontsize=9)
plt.tight_layout()
plt.savefig('roc_comparison.png', bbox_inches='tight')
plt.show()

Attribute	Value
Rows	7,043
Columns	21
Target	`Churn` (Yes / No)
Numeric features	tenure, MonthlyCharges, TotalCharges
Categorical features	gender, SeniorCitizen, Partner, Dependents, PhoneService, MultipleLines, InternetService, OnlineSecurity, OnlineBackup, DeviceProtection, TechSupport, StreamingTV, StreamingMovies, Contract, PaperlessBilling, PaymentMethod

Metric	Why I chose it
ROC AUC	The retention team needs a ranked list, not just yes/no. AUC measures how well the model orders customers from most to least risky.
F1 (weighted)	More honest than accuracy when classes are imbalanced. Accounts for both false positives and false negatives.
Accuracy	Included because it's on the rubric, but it's the least meaningful metric here.

Limitation	Why It Matters	What I'd Do Next
Static snapshot of data	Customer behavior changes over time — a model trained on old data might not work as well later	Retrain monthly and track AUC over time
No customer service data	Call center interactions are probably a strong churn signal that isn't in this dataset	Add complaint/support ticket data if available
Binary outcome only	Doesn't tell you when someone will churn, just that they might	Look into survival analysis for time-to-churn modeling
Default 0.5 threshold	The decision cutoff isn't optimized for the actual business cost of missing a churner	Tune the threshold based on the cost of a false negative vs. false positive
Limited feature engineering	I didn't try creating new features from combinations of existing ones	Explore interaction terms like high charges + no tech support

Predicting Customer Churn in Telecom¶

1. Introduction¶

The Problem¶

Research Question¶

Why This Actually Matters¶

2. Data Source & Description¶

Loading the Data¶

What's in the Dataset¶

3. Data Preparation & EDA¶

Visualization 1: Class Balance & Monthly Charges¶

Visualization 2: Churn Rate by Contract Type & Tenure¶

How EDA Shaped My Modeling Decisions¶

4. Train/Test Split & Evaluation Plan¶

5. Preprocessing Pipeline¶

6. Baseline Model – Logistic Regression¶

7. Model Tuning & Comparison¶

Model 1: K-Nearest Neighbors (KNN)¶

Model 2: Decision Tree¶

Model 3: Random Forest¶

8. Model Comparison¶

What the Results Tell Me¶

9. Story & Limitations¶

What I Found — Plain English¶

What a Company Could Actually Do With This¶

Limitations¶

Appendix: Reproducibility Checklist¶