Breast Cancer Wisconsin Diagnostic — ML Analysis

import os
import sys
import warnings
import pickle

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import shap

from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import (
    accuracy_score, precision_score, recall_score,
    f1_score, roc_auc_score, confusion_matrix,
    roc_curve, precision_recall_curve,
)

warnings.filterwarnings('ignore')

# Add src to path
SRC_DIR = os.path.join(os.path.dirname(os.path.abspath('.')), 'breast-cancer-ml-analysis-2', 'src')
BASE_DIR = '/home/gabriele/cancer-research/breast-cancer-ml-analysis-2'
sys.path.insert(0, os.path.join(BASE_DIR, 'src'))

from preprocessing import load_data, encode_target, split_features_target, get_train_test_split, fit_scaler, apply_scaler
from evaluation import compute_metrics, build_metrics_table, plot_confusion_matrix, plot_roc_curves, plot_precision_recall_curve, plot_threshold_analysis

RANDOM_STATE = 42
INPUT_PATH = os.path.join(BASE_DIR, 'input', 'wdbc.csv')
PLOTS_DIR = os.path.join(BASE_DIR, 'outputs', 'plots')
MODELS_DIR = os.path.join(BASE_DIR, 'outputs', 'models')
PLOT_STYLE = 'seaborn-v0_8-whitegrid'
plt.style.use(PLOT_STYLE)

print('Libraries loaded successfully')

# Load and inspect raw data
df_raw = load_data(INPUT_PATH, id_cols=['id'])
print(f'Shape: {df_raw.shape}')
print(f'Columns: {df_raw.columns.tolist()}')
df_raw.head()

# Data types and basic statistics
print('Data types:')
print(df_raw.dtypes)
print(f'\nMissing values: {df_raw.isnull().sum().sum()}')

# Summary statistics for numeric features
df_raw.describe().round(3)

# Class balance
class_counts = df_raw['diagnosis'].value_counts()
class_pct = df_raw['diagnosis'].value_counts(normalize=True) * 100
print('Class distribution:')
for label in ['B', 'M']:
    print(f'  {label}: {class_counts[label]} ({class_pct[label]:.1f}%)')
print(f'\nClass ratio (M/B): {class_counts["M"] / class_counts["B"]:.3f}')
print('Dataset is moderately imbalanced (63% benign, 37% malignant).')
print('This makes Recall, F1, and ROC-AUC more informative than raw Accuracy.')

# Encode target for plotting
df = encode_target(df_raw, target_col='diagnosis')
X, y = split_features_target(df, target_col='diagnosis')
feature_names = X.columns.tolist()

# Class distribution bar chart
plt.style.use(PLOT_STYLE)
fig, ax = plt.subplots(figsize=(5, 4))
counts = y.value_counts()
ax.bar(['Benign (0)', 'Malignant (1)'], [counts[0], counts[1]], color=['steelblue', 'tomato'])
ax.set_title('Class Distribution')
ax.set_ylabel('Count')
for i, v in enumerate([counts[0], counts[1]]):
    ax.text(i, v + 3, str(v), ha='center', fontweight='bold')
plt.tight_layout()
fig.savefig(os.path.join(PLOTS_DIR, 'class_distribution.png'), dpi=150, bbox_inches='tight')
plt.show()
print('Class distribution: 357 Benign (62.7%), 212 Malignant (37.3%)')

# Feature distributions by class (first 6 mean features)
mean_features = [c for c in feature_names if c.endswith('_mean')][:6]
fig, axes = plt.subplots(2, 3, figsize=(14, 8))
axes = axes.flatten()
for i, feat in enumerate(mean_features):
    for label, color, name in [(0, 'steelblue', 'Benign'), (1, 'tomato', 'Malignant')]:
        axes[i].hist(df[df['diagnosis'] == label][feat], bins=30, alpha=0.6,
                     color=color, label=name, density=True)
    axes[i].set_title(feat)
    axes[i].legend(fontsize=8)
plt.suptitle('Feature Distributions by Class (Mean Features)', fontsize=13, y=1.02)
plt.tight_layout()
fig.savefig(os.path.join(PLOTS_DIR, 'feature_distributions.png'), dpi=150, bbox_inches='tight')
plt.show()

# Correlation heatmap
fig, ax = plt.subplots(figsize=(16, 13))
corr = df[feature_names].corr()
mask = np.triu(np.ones_like(corr, dtype=bool))
sns.heatmap(corr, mask=mask, cmap='coolwarm', center=0, annot=False,
            linewidths=0.5, ax=ax, vmin=-1, vmax=1)
ax.set_title('Feature Correlation Heatmap', fontsize=14)
plt.tight_layout()
fig.savefig(os.path.join(PLOTS_DIR, 'correlation_heatmap.png'), dpi=150, bbox_inches='tight')
plt.show()
print('Many highly correlated feature groups (e.g., radius/perimeter/area).')
print('This multi-collinearity will penalize linear models but is handled well by tree ensembles.')

# Pairplot of top 4 discriminative features
top4 = ['radius_mean', 'texture_mean', 'concave points_mean', 'area_mean']
pair_df = df[top4 + ['diagnosis']].copy()
pair_df['diagnosis'] = pair_df['diagnosis'].map({0: 'Benign', 1: 'Malignant'})
g = sns.pairplot(pair_df, hue='diagnosis',
                 palette={'Benign': 'steelblue', 'Malignant': 'tomato'},
                 plot_kws={'alpha': 0.5}, diag_kind='hist')
g.fig.suptitle('Pairplot of Top 4 Features', y=1.02, fontsize=13)
g.fig.savefig(os.path.join(PLOTS_DIR, 'pairplot_top_features.png'), dpi=150, bbox_inches='tight')
plt.show()
print('Strong linear separation visible in radius_mean vs concave_points_mean.')
print('This suggests linear models (Logistic Regression) may perform well.')

X_train, X_test, y_train, y_test = get_train_test_split(
    X, y, test_size=0.2, random_state=RANDOM_STATE
)

# Fit scaler ONLY on train data — prevents data leakage
scaler = fit_scaler(X_train)
X_train_sc = apply_scaler(scaler, X_train)
X_test_sc = apply_scaler(scaler, X_test)

print(f'Training set: {X_train_sc.shape} | Test set: {X_test_sc.shape}')
print(f'Train class balance: {y_train.mean():.3f} malignant')
print(f'Test class balance:  {y_test.mean():.3f} malignant')
print('Stratification successful: class ratio preserved in both splits.')

# Define baseline models
models = {
    'Logistic Regression': LogisticRegression(max_iter=1000, random_state=RANDOM_STATE),
    'Decision Tree': DecisionTreeClassifier(random_state=RANDOM_STATE),
    'Random Forest': RandomForestClassifier(n_estimators=100, random_state=RANDOM_STATE),
    'Gradient Boosting': GradientBoostingClassifier(n_estimators=100, random_state=RANDOM_STATE),
}

baseline_results = {}
baseline_preds = {}
baseline_probs = {}

for name, model in models.items():
    model.fit(X_train_sc, y_train)
    y_pred = model.predict(X_test_sc)
    y_prob = model.predict_proba(X_test_sc)[:, 1]
    metrics = compute_metrics(y_test, y_pred, y_prob)
    baseline_results[name] = metrics
    baseline_preds[name] = y_pred
    baseline_probs[name] = y_prob
    print(f'{name}: Acc={metrics["Accuracy"]:.4f}, F1={metrics["F1"]:.4f}, AUC={metrics["ROC-AUC"]:.4f}')

# Baseline metrics table
baseline_df = build_metrics_table(baseline_results)
print('Baseline Model Performance:')
baseline_df

# GridSearchCV — Random Forest
rf_param_grid = {
    'n_estimators': [100, 200],
    'max_depth': [None, 5, 10],
    'min_samples_split': [2, 5],
}
rf_grid = GridSearchCV(
    RandomForestClassifier(random_state=RANDOM_STATE),
    rf_param_grid, cv=5, scoring='roc_auc', n_jobs=-1, verbose=0
)
rf_grid.fit(X_train_sc, y_train)
print(f'Best RF params: {rf_grid.best_params_}')
print(f'Best RF CV ROC-AUC: {rf_grid.best_score_:.4f}')

# GridSearchCV — Gradient Boosting
gb_param_grid = {
    'n_estimators': [100, 200],
    'learning_rate': [0.05, 0.1],
    'max_depth': [3, 5],
}
gb_grid = GridSearchCV(
    GradientBoostingClassifier(random_state=RANDOM_STATE),
    gb_param_grid, cv=5, scoring='roc_auc', n_jobs=-1, verbose=0
)
gb_grid.fit(X_train_sc, y_train)
print(f'Best GB params: {gb_grid.best_params_}')
print(f'Best GB CV ROC-AUC: {gb_grid.best_score_:.4f}')

# Evaluate tuned models
tuned_models = {
    'Random Forest (Tuned)': rf_grid.best_estimator_,
    'Gradient Boosting (Tuned)': gb_grid.best_estimator_,
}
tuned_results = {}
tuned_preds = {}
tuned_probs = {}

for name, model in tuned_models.items():
    y_pred = model.predict(X_test_sc)
    y_prob = model.predict_proba(X_test_sc)[:, 1]
    metrics = compute_metrics(y_test, y_pred, y_prob)
    tuned_results[name] = metrics
    tuned_preds[name] = y_pred
    tuned_probs[name] = y_prob
    print(f'{name}: Acc={metrics["Accuracy"]:.4f}, F1={metrics["F1"]:.4f}, AUC={metrics["ROC-AUC"]:.4f}')

# Combined results
all_results = {**baseline_results, **tuned_results}
all_preds = {**baseline_preds, **tuned_preds}
all_probs = {**baseline_probs, **tuned_probs}
all_df = build_metrics_table(all_results)
print('All Models — Complete Metrics Table:')
all_df

# Select best model by ROC-AUC
best_model_name = all_df['ROC-AUC'].idxmax()
print(f'Best model by ROC-AUC: {best_model_name}')
print(f'ROC-AUC: {all_df.loc[best_model_name, "ROC-AUC"]:.4f}')

model_obj_map = {
    'Logistic Regression': models['Logistic Regression'],
    'Decision Tree': models['Decision Tree'],
    'Random Forest': models['Random Forest'],
    'Gradient Boosting': models['Gradient Boosting'],
    'Random Forest (Tuned)': rf_grid.best_estimator_,
    'Gradient Boosting (Tuned)': gb_grid.best_estimator_,
}
best_model = model_obj_map[best_model_name]
best_probs = all_probs[best_model_name]
best_preds = all_preds[best_model_name]

# Save best model
with open(os.path.join(MODELS_DIR, 'best_model.pkl'), 'wb') as f:
    pickle.dump(best_model, f)
print(f'Best model saved to outputs/models/best_model.pkl')

# Display complete metrics table
print('Complete Metrics — All Models (Baseline + Tuned):')
all_df.style.highlight_max(color='lightgreen', axis=0)

# Confusion matrices for all models
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
axes = axes.flatten()
model_names = list(all_results.keys())
class_names = ['Benign (0)', 'Malignant (1)']

for idx, name in enumerate(model_names):
    cm = confusion_matrix(y_test, all_preds[name])
    sns.heatmap(cm, annot=True, fmt='d', cmap='Blues',
                xticklabels=class_names, yticklabels=class_names, ax=axes[idx])
    axes[idx].set_title(name, fontsize=10)
    axes[idx].set_xlabel('Predicted')
    axes[idx].set_ylabel('True')

plt.suptitle('Confusion Matrices — All Models', fontsize=14, y=1.02)
plt.tight_layout()
fig.savefig(os.path.join(PLOTS_DIR, 'confusion_matrices_all.png'), dpi=150, bbox_inches='tight')
plt.show()

# Also save individual confusion matrices
for name in model_names:
    safe = name.lower().replace(' ', '_').replace('(', '').replace(')', '')
    plot_confusion_matrix(y_test, all_preds[name], model_name=name,
                          save_path=os.path.join(PLOTS_DIR, f'confusion_matrix_{safe}.png'))
    plt.close('all')
print('Confusion matrices saved.')

# ROC curves — all models
plt.style.use(PLOT_STYLE)
fig, ax = plt.subplots(figsize=(8, 6))
colors = ['steelblue', 'darkorange', 'green', 'red', 'purple', 'brown']
for (name, y_prob), color in zip(all_probs.items(), colors):
    fpr, tpr, _ = roc_curve(y_test, y_prob)
    auc = roc_auc_score(y_test, y_prob)
    ls = '--' if 'Tuned' in name else '-'
    ax.plot(fpr, tpr, lw=2, linestyle=ls, color=color, label=f'{name} (AUC={auc:.4f})')
ax.plot([0, 1], [0, 1], 'k--', lw=1, label='Random classifier')
ax.set_xlabel('False Positive Rate')
ax.set_ylabel('True Positive Rate')
ax.set_title('ROC Curves — All Models')
ax.legend(loc='lower right', fontsize=8)
plt.tight_layout()
fig.savefig(os.path.join(PLOTS_DIR, 'roc_curves.png'), dpi=150, bbox_inches='tight')
plt.show()

# Precision-Recall curve
plt.style.use(PLOT_STYLE)
precision_vals, recall_vals, thresholds = precision_recall_curve(y_test, best_probs)
fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(recall_vals, precision_vals, lw=2, color='steelblue')
ax.set_xlabel('Recall')
ax.set_ylabel('Precision')
ax.set_title(f'Precision-Recall Curve — {best_model_name}')
ax.set_xlim([0, 1])
ax.set_ylim([0, 1.05])
plt.tight_layout()
fig.savefig(os.path.join(PLOTS_DIR, 'precision_recall_curve.png'), dpi=150, bbox_inches='tight')
plt.show()

# Precision and Recall vs Threshold
plt.style.use(PLOT_STYLE)
thresholds_full = np.append(thresholds, 1.0)
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(thresholds_full, precision_vals, lw=2, label='Precision', color='darkorange')
ax.plot(thresholds_full, recall_vals, lw=2, label='Recall', color='steelblue')
ax.axvline(x=0.3, color='red', linestyle='--', lw=1.5, label='Screening threshold (0.3)')
ax.axvline(x=0.5, color='green', linestyle='--', lw=1.5, label='Default threshold (0.5)')
ax.set_xlabel('Classification Threshold')
ax.set_ylabel('Score')
ax.set_title(f'Precision & Recall vs Threshold — {best_model_name}')
ax.legend()
ax.set_xlim([0, 1])
ax.set_ylim([0, 1.05])
plt.tight_layout()
fig.savefig(os.path.join(PLOTS_DIR, 'threshold_analysis.png'), dpi=150, bbox_inches='tight')
plt.show()

# Clinical threshold recommendation
clinical_threshold = 0.3
best_preds_clinical = (best_probs >= clinical_threshold).astype(int)
metrics_default = compute_metrics(y_test, best_preds, best_probs)
metrics_clinical = compute_metrics(y_test, best_preds_clinical, best_probs)

print(f'=== {best_model_name} ===')
print(f'\nDefault threshold (0.5):')
for k, v in metrics_default.items():
    print(f'  {k}: {v:.4f}')
print(f'\nScreening threshold (0.3):')
for k, v in metrics_clinical.items():
    print(f'  {k}: {v:.4f}')
print('\nAt threshold 0.3: Recall improves at the cost of slightly lower Precision.')
print('Recommended for population screening to minimise missed malignancies.')

# Confusion matrix at clinical threshold
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
class_names = ['Benign (0)', 'Malignant (1)']
for ax, preds, title in [
    (axes[0], best_preds, f'{best_model_name}\n(Default threshold 0.5)'),
    (axes[1], best_preds_clinical, f'{best_model_name}\n(Screening threshold 0.3)'),
]:
    cm = confusion_matrix(y_test, preds)
    sns.heatmap(cm, annot=True, fmt='d', cmap='Blues',
                xticklabels=class_names, yticklabels=class_names, ax=ax)
    ax.set_title(title, fontsize=10)
    ax.set_xlabel('Predicted')
    ax.set_ylabel('True')
plt.tight_layout()
fig.savefig(os.path.join(PLOTS_DIR, 'confusion_matrix_clinical_threshold.png'), dpi=150, bbox_inches='tight')
plt.show()

# SHAP analysis for best model
print(f'Computing SHAP values for: {best_model_name}')

if 'Gradient Boosting' in best_model_name or 'Random Forest' in best_model_name:
    explainer = shap.TreeExplainer(best_model)
    shap_values = explainer.shap_values(X_test_sc)
    if isinstance(shap_values, list):
        sv = shap_values[1]
    else:
        sv = shap_values
else:
    explainer = shap.LinearExplainer(best_model, X_train_sc, feature_perturbation='interventional')
    shap_values = explainer.shap_values(X_test_sc)
    sv = shap_values

print(f'SHAP values shape: {sv.shape}')

# SHAP summary beeswarm plot
plt.figure(figsize=(10, 8))
shap.summary_plot(sv, X_test_sc, feature_names=feature_names, show=False)
plt.title(f'SHAP Summary Plot (Beeswarm) — {best_model_name}', fontsize=13)
plt.tight_layout()
plt.savefig(os.path.join(PLOTS_DIR, 'shap_summary.png'), dpi=150, bbox_inches='tight')
plt.show()
print('SHAP beeswarm: each dot is one test sample. Red = high feature value, blue = low.')
print('Positive SHAP value -> pushes prediction toward Malignant.')

# SHAP bar plot — mean absolute SHAP values (global feature importance)
plt.figure(figsize=(9, 7))
shap.summary_plot(sv, X_test_sc, feature_names=feature_names, plot_type='bar', show=False)
plt.title(f'SHAP Global Feature Importance (Mean |SHAP|) — {best_model_name}', fontsize=12)
plt.tight_layout()
plt.savefig(os.path.join(PLOTS_DIR, 'shap_bar.png'), dpi=150, bbox_inches='tight')
plt.show()

# Top 10 features by mean |SHAP|
mean_abs_shap = np.abs(sv).mean(axis=0)
shap_importance = pd.Series(mean_abs_shap, index=feature_names).sort_values(ascending=False)
print('Top 10 features by mean |SHAP| value:')
print(shap_importance.head(10).round(4).to_string())
print('\nClinical interpretation:')
print('  - concave_points_worst / concavity_worst: measure nuclear irregularity — higher in malignant cells')
print('  - radius_worst / perimeter_worst / area_worst: nuclear size — malignant nuclei are larger')
print('  - texture_worst: nuclear texture (standard deviation of grayscale values) — coarser in cancer')
print('These findings are consistent with established cytological criteria for malignancy.')

# Final summary printout
print('=' * 60)
print('FINAL RESULTS SUMMARY')
print('=' * 60)
print(f'Best model: {best_model_name}')
print(f'Test set metrics:')
for k, v in all_results[best_model_name].items():
    print(f'  {k}: {v:.4f}')
print(f'\nAll outputs saved to outputs/plots/ and outputs/models/')
print(f'Best model pickle: outputs/models/best_model.pkl')
print('=' * 60)