03 · Model Selection and Evaluation¶

Model evaluation is the discipline that separates practitioners who ship reliable systems from those who ship impressive-looking notebooks. Every methodological shortcut here — leaking test data, optimizing a single metric, never plotting residuals — produces models that fail in production in ways that are embarrassingly avoidable. This file treats evaluation not as a checklist but as a statistical testing problem with explicit failure modes.

The Evaluation Hierarchy¶

Before any code: the mental model you need to operate cleanly.

ALL DATA
│
├── Training Set (~60-70%)     ← model sees gradients here
├── Validation Set (~15-20%)   ← hyperparameter selection, model comparison
└── Test Set (~15-20%)         ← touched ONCE, at the very end

The cardinal rule: The test set is a time capsule. You seal it before the project begins and open it exactly once. Any decision made after looking at test performance is a decision made with information from the future. That’s data leakage at the evaluation level — subtler than feature leakage but equally poisonous.

⚠️ What Most Practitioners Get Wrong: The p-hacking equivalent in ML is running dozens of model configurations, picking the best test score, and reporting it as the model’s performance. The test set has implicitly become a second validation set. Every model you try on test data costs you roughly 0.01–0.05 AUC in honest expected generalization performance (the exact number depends on variance and sample size). Stop touching it.

1. Train / Validation / Test Splits¶

The Basics¶

from sklearn.model_selection import train_test_split
import numpy as np

# Standard 70/15/15 split
X_train, X_temp, y_train, y_temp = train_test_split(
    X, y, test_size=0.30, random_state=42, stratify=y
)
X_val, X_test, y_val, y_test = train_test_split(
    X_temp, y_temp, test_size=0.50, random_state=42, stratify=y_temp
)

print(f"Train: {X_train.shape[0]}, Val: {X_val.shape[0]}, Test: {X_test.shape[0]}")

When Standard Splits Break¶

Scenario	Problem	Correct Approach
Time series data	Random split leaks future into past	Time-ordered split (no shuffling)
Grouped data (patients, users)	Same entity in train and test	GroupKFold — split by group
Class imbalance < 1%	Rare class disappears from val	StratifiedKFold
Small dataset (< 1000)	Any fixed split is high variance	k-fold cross-validation
Distribution shift expected	Val set may not reflect production	Covariate-shift-aware splitting

2. Cross-Validation¶

When you have limited data, a single train/val split has high variance — you might get lucky or unlucky. k-fold CV fixes this by using all data for both training and validation, in rotation.

k-Fold Cross-Validation¶

from sklearn.model_selection import KFold, StratifiedKFold, cross_val_score
from sklearn.ensemble import RandomForestClassifier
import numpy as np

# Basic k-fold
kf = KFold(n_splits=5, shuffle=True, random_state=42)
model = RandomForestClassifier(n_estimators=100, random_state=42)

scores = cross_val_score(model, X, y, cv=kf, scoring='roc_auc')
print(f"CV AUC: {scores.mean():.4f} ± {scores.std():.4f}")
# e.g. CV AUC: 0.8821 ± 0.0134

# Stratified k-fold (preserves class proportions in each fold)
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
scores_strat = cross_val_score(model, X, y, cv=skf, scoring='roc_auc')
print(f"Stratified CV AUC: {scores_strat.mean():.4f} ± {scores_strat.std():.4f}")

Nested Cross-Validation (for Hyperparameter Tuning + Evaluation)¶

The subtle trap: if you use the same CV loop to both tune hyperparameters AND estimate performance, you’re overfitting to the outer folds. Nested CV separates them.

from sklearn.model_selection import GridSearchCV, cross_val_score

# Inner loop: hyperparameter tuning
inner_cv = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
param_grid = {'n_estimators': [50, 100, 200], 'max_depth': [3, 5, None]}

grid_search = GridSearchCV(
    RandomForestClassifier(random_state=42),
    param_grid,
    cv=inner_cv,
    scoring='roc_auc'
)

# Outer loop: performance estimation
outer_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
nested_scores = cross_val_score(grid_search, X, y, cv=outer_cv, scoring='roc_auc')

print(f"Nested CV AUC: {nested_scores.mean():.4f} ± {nested_scores.std():.4f}")
# This is an honest estimate of generalization performance

CV for Time Series¶

from sklearn.model_selection import TimeSeriesSplit

tscv = TimeSeriesSplit(n_splits=5)
for fold, (train_idx, val_idx) in enumerate(tscv.split(X_time)):
    X_tr, X_v = X_time[train_idx], X_time[val_idx]
    y_tr, y_v = y_time[train_idx], y_time[val_idx]
    # Train on past, validate on future — always

3. The Bias-Variance Tradeoff (As a Mathematical Tension, Not a Buzzword)¶

Everyone cites bias-variance. Almost no one applies it. Here’s the actual decomposition.

The Decomposition¶

For any model f̂ trained on dataset D, the expected test error decomposes as:

E[(y - f̂(x))²] = Bias² + Variance + Irreducible Noise

Where:

Bias² = (E[f̂(x)] - f(x))² — how wrong the model is on average across all possible training sets
Variance = E[(f̂(x) - E[f̂(x)])²] — how much the model changes across different training sets
Irreducible Noise = inherent noise in the data generation process; you cannot reduce this

Visualizing It¶

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline

np.random.seed(42)
n_samples = 30
X_true = np.linspace(0, 1, 200)
y_true = np.sin(2 * np.pi * X_true)

def simulate_bias_variance(degree, n_repeats=50):
    """Train many models on different samples, measure bias and variance."""
    predictions = np.zeros((n_repeats, len(X_true)))
    
    for i in range(n_repeats):
        X_sample = np.random.uniform(0, 1, n_samples)
        y_sample = np.sin(2 * np.pi * X_sample) + np.random.normal(0, 0.3, n_samples)
        
        model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
        model.fit(X_sample.reshape(-1, 1), y_sample)
        predictions[i] = model.predict(X_true.reshape(-1, 1))
    
    mean_pred = predictions.mean(axis=0)
    bias_sq = np.mean((mean_pred - y_true) ** 2)
    variance = np.mean(predictions.var(axis=0))
    
    return bias_sq, variance, bias_sq + variance

degrees = [1, 2, 4, 9, 15]
results = []
for d in degrees:
    b, v, total = simulate_bias_variance(d)
    results.append({'degree': d, 'bias²': b, 'variance': v, 'total': total})
    print(f"Degree {d:2d}: Bias²={b:.4f}, Variance={v:.4f}, Total={total:.4f}")

# Output shows the U-shape: degree 1 high bias, degree 15 high variance

In Practice: Diagnosing Your Model¶

from sklearn.model_selection import learning_curve

train_sizes, train_scores, val_scores = learning_curve(
    model, X_train, y_train,
    cv=5,
    train_sizes=np.linspace(0.1, 1.0, 10),
    scoring='neg_mean_squared_error',
    n_jobs=-1
)

train_mean = -train_scores.mean(axis=1)
val_mean = -val_scores.mean(axis=1)

plt.figure(figsize=(10, 5))
plt.plot(train_sizes, train_mean, label='Training error')
plt.plot(train_sizes, val_mean, label='Validation error')
plt.xlabel('Training set size')
plt.ylabel('MSE')
plt.legend()
plt.title('Learning Curve — Diagnose Bias vs Variance')
plt.tight_layout()
plt.savefig('learning_curve.png', dpi=150)

Interpretation table:

Pattern	Diagnosis	Fix
Both errors high, converge together (high)	High bias (underfitting)	More complex model, more features, less regularization
Train error low, val error much higher	High variance (overfitting)	More data, regularization, simpler model, ensemble
Both errors low, converge together	Well-fit	Deploy it
Val error decreasing as data grows	Variance-dominated, get more data	More training data will help
Val error plateaued high	Bias-dominated	Model complexity issue, not a data issue

4. Evaluation Metrics¶

Choosing a metric is choosing what you optimize for. This is a business decision dressed as a technical one.

Classification Metrics¶

from sklearn.metrics import (
    accuracy_score, precision_score, recall_score, f1_score,
    roc_auc_score, average_precision_score, confusion_matrix,
    classification_report, brier_score_loss
)
import matplotlib.pyplot as plt
from sklearn.calibration import calibration_curve

def full_classification_report(y_true, y_pred, y_prob):
    """
    Complete evaluation for a binary classifier.
    y_prob: probability of positive class
    """
    print("=== Classification Metrics ===")
    print(f"Accuracy:          {accuracy_score(y_true, y_pred):.4f}")
    print(f"Precision:         {precision_score(y_true, y_pred):.4f}")
    print(f"Recall (Sens.):    {recall_score(y_true, y_pred):.4f}")
    print(f"F1 Score:          {f1_score(y_true, y_pred):.4f}")
    print(f"ROC-AUC:           {roc_auc_score(y_true, y_prob):.4f}")
    print(f"PR-AUC (Avg Prec): {average_precision_score(y_true, y_prob):.4f}")
    print(f"Brier Score:       {brier_score_loss(y_true, y_prob):.4f}")
    print()
    print(classification_report(y_true, y_pred))
    
    # Confusion matrix
    cm = confusion_matrix(y_true, y_pred)
    tn, fp, fn, tp = cm.ravel()
    print(f"TN={tn}, FP={fp}, FN={fn}, TP={tp}")
    print(f"Specificity (TNR): {tn/(tn+fp):.4f}")
    print(f"NPV:               {tn/(tn+fn):.4f}")

When to Use Which Metric¶

Metric	Use When	Misleading When
Accuracy	Classes balanced, equal misclass cost	Imbalanced classes (99% negative → 99% accuracy trivially)
Precision	FP is expensive (spam filter — false positive = deleted legit email)	You also care about FN
Recall	FN is expensive (cancer screening — missing a case is catastrophic)	You also care about FP
F1	Balance P and R; imbalanced classes	One error type is much worse than the other
ROC-AUC	Ranking quality; want threshold-agnostic metric	High class imbalance (PR-AUC is better here)
PR-AUC	Imbalanced classes; positive class is rare	You need absolute probability estimates
Brier Score	Probabilistic calibration quality	Pure ranking is enough
Log Loss	Training objective for calibrated classifiers	Outliers: a confident wrong prediction explodes log loss

ROC Curve vs PR Curve¶

from sklearn.metrics import roc_curve, precision_recall_curve
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# ROC Curve
fpr, tpr, thresholds_roc = roc_curve(y_test, y_prob)
auc_roc = roc_auc_score(y_test, y_prob)
axes[0].plot(fpr, tpr, label=f'ROC AUC = {auc_roc:.3f}')
axes[0].plot([0,1],[0,1], 'k--', label='Random classifier')
axes[0].set_xlabel('False Positive Rate')
axes[0].set_ylabel('True Positive Rate (Recall)')
axes[0].set_title('ROC Curve')
axes[0].legend()

# Precision-Recall Curve
precision, recall, thresholds_pr = precision_recall_curve(y_test, y_prob)
auc_pr = average_precision_score(y_test, y_prob)
axes[1].plot(recall, precision, label=f'PR AUC = {auc_pr:.3f}')
base = y_test.mean()  # prevalence
axes[1].axhline(y=base, color='k', linestyle='--', label=f'Random ({base:.3f})')
axes[1].set_xlabel('Recall')
axes[1].set_ylabel('Precision')
axes[1].set_title('Precision-Recall Curve')
axes[1].legend()

plt.tight_layout()
plt.savefig('roc_pr_curves.png', dpi=150)

Model Calibration — The Most Neglected Metric¶

A model can have a perfect ROC-AUC and still be completely miscalibrated. If your model outputs 0.9 probability for events that happen 50% of the time, it’s useless for any downstream decision that relies on probabilities.

from sklearn.calibration import calibration_curve, CalibratedClassifierCV
from sklearn.ensemble import GradientBoostingClassifier
import matplotlib.pyplot as plt

# Calibration curve
prob_true, prob_pred = calibration_curve(y_test, y_prob, n_bins=10)

plt.figure(figsize=(8, 6))
plt.plot([0, 1], [0, 1], 'k--', label='Perfect calibration')
plt.plot(prob_pred, prob_true, 's-', label='Model (uncalibrated)')

# Fix calibration with Platt scaling or isotonic regression
from sklearn.calibration import CalibratedClassifierCV
cal_model = CalibratedClassifierCV(base_model, method='isotonic', cv=5)
cal_model.fit(X_train, y_train)
y_prob_cal = cal_model.predict_proba(X_test)[:, 1]

prob_true_cal, prob_pred_cal = calibration_curve(y_test, y_prob_cal, n_bins=10)
plt.plot(prob_pred_cal, prob_true_cal, 'o-', label='Model (calibrated)')

plt.xlabel('Mean Predicted Probability')
plt.ylabel('Fraction of Positives')
plt.title('Calibration Curve')
plt.legend()
plt.savefig('calibration_curve.png', dpi=150)

print(f"Brier Score (uncalibrated): {brier_score_loss(y_test, y_prob):.4f}")
print(f"Brier Score (calibrated):   {brier_score_loss(y_test, y_prob_cal):.4f}")

Regression Metrics¶

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
import numpy as np

def regression_report(y_true, y_pred):
    mse = mean_squared_error(y_true, y_pred)
    rmse = np.sqrt(mse)
    mae = mean_absolute_error(y_true, y_pred)
    r2 = r2_score(y_true, y_pred)
    # MAPE: careful — divide by zero if any y_true == 0
    mape = np.mean(np.abs((y_true - y_pred) / (y_true + 1e-8))) * 100
    
    print(f"RMSE:  {rmse:.4f}  (same units as target)")
    print(f"MAE:   {mae:.4f}   (robust to outliers)")
    print(f"R²:    {r2:.4f}   (1.0 = perfect, 0 = mean predictor, <0 = worse than mean)")
    print(f"MAPE:  {mape:.2f}%  (percent error, avoid if target near 0)")
    
    # Residual analysis — non-random residuals = model misspecification
    residuals = y_true - y_pred
    print(f"\nResidual Analysis:")
    print(f"  Mean:    {residuals.mean():.4f}  (should be ~0)")
    print(f"  Std:     {residuals.std():.4f}")
    print(f"  Skew:    {pd.Series(residuals).skew():.4f}  (should be ~0)")

5. Hyperparameter Tuning¶

Grid Search vs Random Search vs Bayesian Optimization¶

from sklearn.model_selection import GridSearchCV, RandomizedSearchCV
from scipy.stats import randint, uniform

# Grid Search: exhaustive, works for small grids
param_grid = {
    'n_estimators': [100, 200, 500],
    'max_depth': [3, 5, 10, None],
    'min_samples_split': [2, 5, 10]
}
# Total: 3 × 4 × 3 = 36 configurations × 5-fold = 180 model fits

grid_cv = GridSearchCV(
    RandomForestClassifier(random_state=42),
    param_grid,
    cv=StratifiedKFold(5),
    scoring='roc_auc',
    n_jobs=-1,
    verbose=1
)
grid_cv.fit(X_train, y_train)
print(f"Best params: {grid_cv.best_params_}")
print(f"Best CV AUC: {grid_cv.best_score_:.4f}")

# Random Search: better for large spaces (Bergstra & Bengio 2012: 60 random > full grid)
param_dist = {
    'n_estimators': randint(50, 500),
    'max_depth': [3, 5, 10, None],
    'min_samples_split': randint(2, 20),
    'max_features': uniform(0.3, 0.7)
}

random_cv = RandomizedSearchCV(
    RandomForestClassifier(random_state=42),
    param_dist,
    n_iter=50,  # 50 random configurations
    cv=StratifiedKFold(5),
    scoring='roc_auc',
    n_jobs=-1,
    random_state=42
)
random_cv.fit(X_train, y_train)
print(f"Random Search Best AUC: {random_cv.best_score_:.4f}")

# Bayesian Optimization (when CV is expensive): use optuna
try:
    import optuna
    optuna.logging.set_verbosity(optuna.logging.WARNING)
    
    def objective(trial):
        params = {
            'n_estimators': trial.suggest_int('n_estimators', 50, 500),
            'max_depth': trial.suggest_int('max_depth', 2, 15),
            'min_samples_split': trial.suggest_int('min_samples_split', 2, 20),
            'max_features': trial.suggest_float('max_features', 0.3, 1.0)
        }
        model = RandomForestClassifier(**params, random_state=42, n_jobs=-1)
        score = cross_val_score(model, X_train, y_train, cv=3, scoring='roc_auc').mean()
        return score
    
    study = optuna.create_study(direction='maximize')
    study.optimize(objective, n_trials=50, show_progress_bar=True)
    print(f"Optuna Best AUC: {study.best_value:.4f}")
    print(f"Optuna Best Params: {study.best_params}")
except ImportError:
    print("pip install optuna for Bayesian optimization")

Speed Comparison for Reference¶

On a typical tabular dataset with n=10,000, p=50:

Grid search (36 configs × 5-fold): ~180 model fits, ~45 seconds
Random search (50 configs × 5-fold): ~250 model fits, ~62 seconds
Optuna (50 trials × 3-fold): ~150 model fits, ~38 seconds, typically better result

6. Overfitting/Underfitting Diagnostics Checklist¶

Run this checklist before declaring a model “done.”

def model_health_check(model, X_train, y_train, X_val, y_val, task='classification'):
    """
    Systematic overfitting/underfitting diagnostics.
    """
    from sklearn.metrics import roc_auc_score, r2_score
    
    if task == 'classification':
        train_score = roc_auc_score(y_train, model.predict_proba(X_train)[:, 1])
        val_score = roc_auc_score(y_val, model.predict_proba(X_val)[:, 1])
        metric_name = 'AUC'
    else:
        train_score = r2_score(y_train, model.predict(X_train))
        val_score = r2_score(y_val, model.predict(X_val))
        metric_name = 'R²'
    
    gap = train_score - val_score
    
    print(f"Train {metric_name}: {train_score:.4f}")
    print(f"Val   {metric_name}: {val_score:.4f}")
    print(f"Gap:                 {gap:.4f}")
    
    if train_score < 0.7 and val_score < 0.7:
        print("⚠️  UNDERFITTING: Both scores low. Increase model complexity.")
    elif gap > 0.05:
        print(f"⚠️  OVERFITTING: Gap {gap:.4f} > 0.05. Add regularization or more data.")
    elif val_score > 0.95:
        print("⚠️  SUSPICIOUSLY HIGH: Check for data leakage.")
    else:
        print("✅ Model appears healthy. Verify on holdout test set.")

7. Statistical Significance of Model Comparisons¶

⚠️ What Most People Get Wrong: Reporting “Model A got 0.84 AUC vs Model B’s 0.83 AUC” without a confidence interval. On a test set of n=500, that gap is not statistically meaningful. Use McNemar’s test for classification or the Diebold-Mariano test for regression.

from scipy.stats import wilcoxon
import numpy as np

def compare_models_statistically(model_a_scores, model_b_scores):
    """
    model_a_scores, model_b_scores: CV scores from cross_val_score()
    Uses Wilcoxon signed-rank test (non-parametric, appropriate for CV scores)
    """
    stat, p_value = wilcoxon(model_a_scores, model_b_scores)
    mean_diff = np.mean(model_a_scores - model_b_scores)
    
    print(f"Model A mean CV score: {model_a_scores.mean():.4f} ± {model_a_scores.std():.4f}")
    print(f"Model B mean CV score: {model_b_scores.mean():.4f} ± {model_b_scores.std():.4f}")
    print(f"Mean difference: {mean_diff:+.4f}")
    print(f"Wilcoxon p-value: {p_value:.4f}")
    
    if p_value < 0.05:
        better = "A" if mean_diff > 0 else "B"
        print(f"✅ Model {better} is statistically significantly better (p={p_value:.4f})")
    else:
        print(f"⚠️  Difference is NOT statistically significant (p={p_value:.4f})")
        print("    The models are equivalent given this evidence. Do not claim superiority.")

Practice Problems¶

Problem 1 — The Leaky Pipeline: Take the UCI Credit Card Default dataset. Build a pipeline that accidentally leaks the test set during preprocessing (e.g., fit a scaler on all data before splitting). Measure the inflated performance. Fix the leak. Quantify the actual vs. inflated AUC difference. Acceptance: you find a gap of at least 0.02-0.05 AUC and document why.

Problem 2 — Calibration Audit: Take any ensemble model (RandomForest, GBM) on a binary classification dataset with class imbalance (e.g., fraud detection dataset). Compute the calibration curve. Apply Platt scaling and isotonic regression. Compare Brier scores before/after. Write a one-paragraph conclusion on which method to prefer and why. Acceptance: Brier score improves by at least 10% with calibration.

Problem 3 — Bias-Variance Budget: On the California Housing dataset, train polynomial regression models of degrees 1 through 10. For each degree, compute train RMSE and 5-fold CV RMSE. Plot both curves. Identify the optimal degree. Write a two-sentence diagnosis for each end of the spectrum. Acceptance: clean plot with both curves, correct identification of underfitting and overfitting regions.

Problem 4 — Statistical Model Comparison: Train an SVM and a Random Forest on the same dataset using 10-fold CV. Use the Wilcoxon signed-rank test to determine if the better-performing model is statistically significantly better. Report the p-value and your conclusion in one sentence. Acceptance: correct test applied, p-value reported, conclusion stated.

Resources¶

Resource	Type	Why	Time
sklearn User Guide: Model Selection	Docs	Canonical reference, exhaustive	Reference
ISL Chapter 5 (Resampling)	Book	Rigorous CV derivation	3 hrs
“Calibration of Probabilities” (Niculescu-Mizil & Caruana 2005)	Paper	Why ensembles are miscalibrated	30 min
Bergstra & Bengio (2012): Random Search vs Grid Search	Paper	Proves random search dominates	20 min
StatQuest: ROC/AUC series	Video	Clearest explanation of AUC semantics	45 min

Return to README.md · Previous: 02_unsupervised_learning.md · Next: 04_feature_engineering.md