Scikit-Learn Advanced Model Selection

Real-World Scenario

A competition entry achieves 94% accuracy on the leaderboard but scores 87% on the private test set. The reason: the team made 50 submissions, selecting the model with the highest public leaderboard score. This is test set overfitting — standard practice in competitions but a critical mistake in production. Nested cross-validation would have given an honest estimate.

Validation Curves

from sklearn.datasets import make_classification
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import validation_curve, StratifiedKFold
import numpy as np

X, y = make_classification(n_samples=2000, n_features=20, n_informative=10, random_state=42)

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# How does performance change as we vary max_depth?
param_range = [1, 2, 3, 4, 5, 7, 10, 15]

train_scores, val_scores = validation_curve(
    GradientBoostingClassifier(n_estimators=100, random_state=42),
    X, y,
    param_name="max_depth",
    param_range=param_range,
    cv=cv,
    scoring="roc_auc",
    n_jobs=-1,
)

print(f"{'depth':>6} {'train AUC':>12} {'val AUC':>12} {'gap':>8}")
print("-" * 42)
for depth, tr, val in zip(
    param_range,
    train_scores.mean(axis=1),
    val_scores.mean(axis=1),
):
    gap    = tr - val
    status = "⚠ overfit" if gap > 0.05 else "✓"
    print(f"{depth:>6} {tr:>12.4f} {val:>12.4f} {gap:>8.4f}  {status}")

# Find optimal depth
best_depth = param_range[val_scores.mean(axis=1).argmax()]
print(f"\nBest depth: {best_depth}")

Stratified Splitting for Multi-Label and Regression

from sklearn.model_selection import (
    StratifiedKFold, StratifiedGroupKFold,
    TimeSeriesSplit, GroupShuffleSplit
)
from sklearn.datasets import make_classification
import numpy as np

# ── Standard stratified split ─────────────────────────────────────
X, y = make_classification(n_samples=1000, n_features=20,
                            weights=[0.9, 0.1], random_state=42)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

print("Stratified 5-fold class balance per fold:")
for fold, (train_idx, val_idx) in enumerate(cv.split(X, y)):
    pos_rate = y[val_idx].mean()
    print(f"  Fold {fold+1}: {len(val_idx)} samples, {pos_rate:.1%} positive")

# ── Group K-Fold: keep groups in the same fold ───────────────────
# Use case: patient data — all visits from one patient must be in the same fold
n = 1000
groups = np.repeat(np.arange(100), 10)  # 100 patients, 10 visits each
X = np.random.randn(n, 5)
y = np.random.randint(0, 2, n)

from sklearn.model_selection import GroupKFold
group_cv = GroupKFold(n_splits=5)

print("\nGroup K-Fold (no patient leaks across folds):")
for fold, (train_idx, val_idx) in enumerate(group_cv.split(X, y, groups)):
    train_groups = set(groups[train_idx])
    val_groups   = set(groups[val_idx])
    leak = train_groups & val_groups
    print(f"  Fold {fold+1}: {len(val_groups)} patients in val, overlap={len(leak)}")

# ── Stratified Group K-Fold: preserve both group integrity and class balance ─
from sklearn.model_selection import StratifiedGroupKFold
sgkf = StratifiedGroupKFold(n_splits=5)

print("\nStratified Group K-Fold:")
for fold, (train_idx, val_idx) in enumerate(sgkf.split(X, y, groups)):
    pos_rate = y[val_idx].mean()
    val_groups = set(groups[val_idx])
    print(f"  Fold {fold+1}: {pos_rate:.1%} positive, {len(val_groups)} groups")

Model Comparison with Paired Tests

from sklearn.datasets import make_classification
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.ensemble import (
    RandomForestClassifier, GradientBoostingClassifier, HistGradientBoostingClassifier
)
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from scipy import stats
import numpy as np

X, y = make_classification(n_samples=3000, n_features=20, n_informative=12, random_state=42)
cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)

models = {
    "LR":  Pipeline([("scaler", StandardScaler()), ("clf", LogisticRegression(max_iter=1000))]),
    "RF":  RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1),
    "GB":  GradientBoostingClassifier(n_estimators=200, random_state=42),
    "HGB": HistGradientBoostingClassifier(max_iter=200, random_state=42),
}

scores = {}
for name, model in models.items():
    s = cross_val_score(model, X, y, cv=cv, scoring="roc_auc", n_jobs=-1)
    scores[name] = s
    print(f"{name:5s}: {s.mean():.4f} ± {s.std():.4f}")

# Paired t-test to check if difference between models is statistically significant
best_model = max(scores, key=lambda k: scores[k].mean())
print(f"\nBest model: {best_model}")
print("\nPaired t-test vs best model (5% significance):")
for name, s in scores.items():
    if name == best_model:
        continue
    t_stat, p_val = stats.ttest_rel(scores[best_model], s)
    sig = "★ significant" if p_val < 0.05 else "not significant"
    print(f"  {best_model} vs {name:5s}: p={p_val:.4f}  {sig}")

Threshold-Free Evaluation

from sklearn.datasets import make_classification
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_predict
from sklearn.metrics import (
    roc_auc_score, average_precision_score,
    roc_curve, precision_recall_curve,
)
import numpy as np

X, y = make_classification(n_samples=3000, n_features=20, weights=[0.85, 0.15], random_state=42)
model = HistGradientBoostingClassifier(max_iter=200, random_state=42)
cv    = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# Get out-of-fold probability predictions (avoids train/test contamination)
oof_probs = cross_val_predict(model, X, y, cv=cv, method="predict_proba")[:, 1]

# Threshold-free metrics
auc_roc = roc_auc_score(y, oof_probs)
avg_prec = average_precision_score(y, oof_probs)

print(f"OOF AUC-ROC:          {auc_roc:.4f}")
print(f"OOF Average Precision: {avg_prec:.4f}")

# ROC curve analysis
fpr, tpr, thresholds = roc_curve(y, oof_probs)
# Find threshold for desired TPR ≥ 0.85
idx_85tpr = np.argmax(tpr >= 0.85)
print(f"\nAt TPR ≥ 85%: threshold={thresholds[idx_85tpr]:.3f}  "
      f"FPR={fpr[idx_85tpr]:.3f}")

# PR curve analysis
prec, rec, pr_thresholds = precision_recall_curve(y, oof_probs)
# Find threshold for desired precision ≥ 0.80
idx_80prec = np.argmax(prec >= 0.80)
if idx_80prec < len(pr_thresholds):
    print(f"At Precision ≥ 80%: threshold={pr_thresholds[idx_80prec]:.3f}  "
          f"Recall={rec[idx_80prec]:.3f}")

# The J statistic (Youden's J): maximize TPR - FPR tradeoff
j_scores  = tpr - fpr
best_j    = np.argmax(j_scores)
print(f"\nYouden's J optimal threshold: {thresholds[best_j]:.3f}  "
      f"TPR={tpr[best_j]:.3f}  FPR={fpr[best_j]:.3f}")

Detecting and Preventing Data Leakage

import pandas as pd
import numpy as np
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.datasets import make_classification

X, y = make_classification(n_samples=2000, n_features=50, n_informative=10, random_state=42)
cv   = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# ── WRONG: fit PCA on all data before cross-validation ──────────────
pca          = PCA(n_components=20)
X_pca_leaky  = pca.fit_transform(X)   # uses validation fold info — LEAKAGE

leaky_scores = cross_val_score(
    LogisticRegression(max_iter=1000),
    X_pca_leaky, y, cv=cv, scoring="roc_auc"
)
print(f"Leaky PCA + LR:   {leaky_scores.mean():.4f} ± {leaky_scores.std():.4f}")

# ── CORRECT: PCA inside pipeline ────────────────────────────────────
correct_pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("pca",    PCA(n_components=20)),
    ("model",  LogisticRegression(max_iter=1000)),
])

correct_scores = cross_val_score(correct_pipe, X, y, cv=cv, scoring="roc_auc")
print(f"Pipeline PCA + LR: {correct_scores.mean():.4f} ± {correct_scores.std():.4f}")

# The gap between leaky and correct is the "optimism bias"
bias = leaky_scores.mean() - correct_scores.mean()
print(f"Optimism bias from leakage: {bias:.4f}")

# ── Other common leakage patterns ───────────────────────────────────
LEAKAGE_CHECKLIST = """
Common data leakage sources:
1. Fitting scalers/encoders on full dataset before CV
2. Target encoding without out-of-fold estimation
3. Using future data as features in time-series (look-ahead bias)
4. Duplicate rows split across train and test
5. Group leakage: related samples (e.g. patient visits) in both folds
6. Feature created from target statistics including the test row itself
7. Pre-filtering data using the target variable before splitting

Prevention:
- Always use sklearn Pipelines for CV
- Use cross_val_predict for target-based features
- Use TimeSeriesSplit for temporal data
- Use GroupKFold when samples are grouped
- Check for duplicates before splitting
"""
print(LEAKAGE_CHECKLIST)