Stacking and Blending

The Production Scenario

The Kaggle competition leaderboard. Top position: 0.9412 AUC. Your single XGBoost model, tuned for three days: 0.9381. The gap feels insurmountable. You have tried every hyperparameter combination, every feature engineering trick in your notebook. You are 0.0031 AUC away from the gold zone and running out of ideas.

A teammate suggests stacking. You are skeptical - stacking sounds like it just averages models together, which you already tried by averaging your XGBoost with a LightGBM and got 0.9389 (a modest improvement). Your teammate shakes their head: "Averaging is not stacking. Stacking learns how to combine the models."

You spend four hours implementing proper stacked generalization with out-of-fold predictions. A logistic regression meta-learner sits on top of five base models - XGBoost, LightGBM, Random Forest, Extra Trees, and a neural network. You submit. The leaderboard updates: 0.9407. You gained 0.0026 AUC from the exact same base features, zero additional data, and no more hyperparameter tuning. The difference was the architecture of how predictions were combined.

This lesson explains exactly how that 0.0026 was earned, why averaging failed where stacking succeeded, what out-of-fold predictions prevent, and when stacking is - and is not - worth the engineering complexity.

Why Ensembles Work: Complementary Errors

The mathematical foundation for all ensemble methods comes back to the variance decomposition:

$\text{Var}(\bar{X}) = \rho \sigma^2 + \frac{(1-\rho)\sigma^2}{B}$

But there is a second, equally important reason ensembles work beyond variance reduction: bias reduction through complementary error patterns. When different models make different types of errors, an ensemble that learns to combine them can correct errors that no individual model can fix.

Scenario A: Three correlated models (all GBM variants)
  Sample 1:  Model1=WRONG, Model2=WRONG, Model3=WRONG  → Ensemble=WRONG
  Sample 2:  Model1=RIGHT, Model2=RIGHT, Model3=RIGHT  → Ensemble=RIGHT
  Correlation ρ ≈ 0.95 → minimal ensemble benefit

Scenario B: Three diverse models (GBM + linear + neural)
  Sample 1:  GBM=WRONG, Linear=RIGHT, NN=RIGHT  → Ensemble=RIGHT (majority)
  Sample 2:  GBM=RIGHT, Linear=WRONG, NN=RIGHT  → Ensemble=RIGHT (majority)
  Sample 3:  GBM=RIGHT, Linear=RIGHT, NN=WRONG  → Ensemble=RIGHT (majority)
  Correlation ρ ≈ 0.55 → substantial ensemble benefit

Stacking goes further: it LEARNS which model to trust for which samples.

Ensemble Methods Taxonomy

Method	Training	Models	How Combined	Key Mechanism
Bagging	Parallel	Same algorithm	Average / vote	Variance reduction via bootstrap
Boosting	Sequential	Same algorithm	Weighted sum	Bias reduction via residual correction
Voting	Parallel	Different algorithms	Fixed weights	Simple combination, no meta-learning
Stacking	Parallel then sequential	Different algorithms	Trained meta-model	Learns optimal combination per sample
Blending	Parallel	Different algorithms	Trained on holdout	Simplified stacking via holdout split

Voting: Simple Ensembles

Hard Voting

Each model casts a class vote. The majority wins.

$\hat{y} = \text{mode}\left(\hat{y}_1, \hat{y}_2, \ldots, \hat{y}_M\right)$

Soft Voting

Average the predicted probabilities, then take the argmax:

$\hat{p} = \frac{1}{M}\sum_{m=1}^{M} \hat{p}_m(y=1 \mid x)$

Soft voting almost always outperforms hard voting because it preserves confidence information. A model predicting 0.99 and a model predicting 0.51 both "vote yes" in hard voting - in soft voting, the confident model gets proportionally more effective weight.

Weighted Voting

$\hat{p} = \frac{\sum_{m=1}^{M} w_m \hat{p}_m}{\sum_{m=1}^{M} w_m}$

Weights $w_m$ are usually set by validation AUC or learned via a constrained optimization. This is exactly where stacking begins - instead of setting weights manually, we train a model to learn them.

:::note Weighted voting is a special case of stacking Weighted voting is stacking where the meta-learner is constrained to be a linear combination with non-negative coefficients summing to one. Stacking relaxes all of these constraints and lets the meta-learner be any model - including one that assigns different weights to different base models depending on the input. :::

Stacked Generalization: Architecture

Introduced by David Wolpert in 1992, stacked generalization formalizes the idea that a model can learn to combine other models' predictions.

Level 0: Base Learners (trained on training set)
─────────────────────────────────────────────────────────────────────
  Model 1: XGBoost                    ─┐
  Model 2: LightGBM                    │── Out-of-Fold (OOF) predictions
  Model 3: Random Forest               │   shape: (n_train, n_models)
  Model 4: Extra Trees                 │   These are meta-features
  Model 5: Neural Network (MLP)       ─┘

                          ↓
Level 1: Meta-Learner (trained on OOF predictions)
─────────────────────────────────────────────────────────────────────
  Meta-model: Logistic Regression (simple by design)
  Input:  OOF predictions from Level 0 (n_train, n_models)
  Output: Final prediction

The meta-learner is almost always kept simple - logistic regression, ridge regression, or a small linear model. Complex meta-learners overfit the meta-features, which is especially problematic because the meta-feature space is small (only n_models columns).

The Data Leakage Problem in Stacking

This is the most commonly misunderstood and most frequently violated aspect of stacking.

The wrong approach (introduces leakage):

1. Train base models on full training set
2. Generate base model predictions on same training set
3. Train meta-learner on those predictions + labels

PROBLEM: Base models predict on data they were trained on.
Their predictions are optimistically overfit - they look perfect
in training but will be much worse at test time.
The meta-learner learns to trust predictions that are
unrealistically good. At test time, when base models make
normal (worse) predictions, the meta-learner is miscalibrated.
Result: spectacular train accuracy, collapsed test accuracy.

The correct approach: out-of-fold predictions

1. Split training set into K folds (typically K=5)
2. For each fold k:
   - Train base models on remaining K-1 folds
   - Predict on fold k (these models NEVER saw fold k)
3. Concatenate fold predictions → OOF predictions (n_train, n_models)
4. Each OOF prediction is made by a model that never saw that sample
5. Train meta-learner on (OOF predictions, y_train)

This is honest: the meta-learner sees realistic base model predictions.

Visual representation of 5-fold OOF for one base model:

Train: [██████████████████████████████████████████████████████████████]
                                           ↓ Split into 5 folds
Fold 1: [████░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░]  train on ○, predict ●
Fold 2: [░░░░████░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░]  train on ○, predict ●
Fold 3: [░░░░░░░░████░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░]  train on ○, predict ●
Fold 4: [░░░░░░░░░░░░████░░░░░░░░░░░░░░░░░░░░░░░░░░░░]  train on ○, predict ●
Fold 5: [░░░░░░░░░░░░░░░░████░░░░░░░░░░░░░░░░░░░░░░░░]  train on ○, predict ●

OOF for this model = concatenate [fold1_pred, fold2_pred, ..., fold5_pred]
Each prediction: made by a model that never saw this sample during training.
Repeat for all M base models → OOF matrix: (n_train, M)

Full OOF Stacking Implementation

import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.model_selection import StratifiedKFold, train_test_split, cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.metrics import roc_auc_score
import xgboost as xgb
import lightgbm as lgb

# ── Dataset ───────────────────────────────────────────────────────────────────
X, y = make_classification(
    n_samples=12_000, n_features=25, n_informative=18,
    n_redundant=4, n_clusters_per_class=2, random_state=42,
)
X = pd.DataFrame(X, columns=[f"f{i}" for i in range(X.shape[1])])
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)

# ── Define base learners ──────────────────────────────────────────────────────
base_models = {
    "xgb": xgb.XGBClassifier(
        n_estimators=400, max_depth=5, learning_rate=0.05,
        subsample=0.8, colsample_bytree=0.8,
        eval_metric="auc", random_state=42, verbosity=0,
    ),
    "lgb": lgb.LGBMClassifier(
        n_estimators=400, max_depth=5, learning_rate=0.05,
        subsample=0.8, colsample_bytree=0.8,
        random_state=42, verbose=-1,
    ),
    "rf": RandomForestClassifier(
        n_estimators=300, max_depth=12, random_state=42, n_jobs=-1,
    ),
    "et": ExtraTreesClassifier(
        n_estimators=300, max_depth=12, random_state=42, n_jobs=-1,
    ),
    "mlp": MLPClassifier(
        hidden_layer_sizes=(256, 128, 64), max_iter=500,
        random_state=42, early_stopping=True, validation_fraction=0.1,
    ),
}


def generate_oof_predictions(
    base_models: dict,
    X: pd.DataFrame,
    y: np.ndarray,
    n_splits: int = 5,
    random_state: int = 42,
) -> np.ndarray:
    """
    Generate out-of-fold predictions for all base models.

    Every training sample's OOF prediction is made by a model
    that never saw that sample during training - no leakage.

    Returns
    -------
    oof_preds : np.ndarray, shape (n_samples, n_models)
    """
    skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=random_state)
    n_models  = len(base_models)
    oof_preds = np.zeros((len(X), n_models))

    for model_idx, (model_name, model) in enumerate(base_models.items()):
        print(f"  Generating OOF for: {model_name}")
        fold_aucs = []

        for fold_idx, (train_idx, val_idx) in enumerate(skf.split(X, y)):
            X_tr_fold = X.iloc[train_idx]
            X_val_fold = X.iloc[val_idx]
            y_tr_fold  = y[train_idx]
            y_val_fold = y[val_idx]

            # Clone model to avoid state leakage between folds
            import copy
            fold_model = copy.deepcopy(model)
            fold_model.fit(X_tr_fold, y_tr_fold)

            val_proba = fold_model.predict_proba(X_val_fold)[:, 1]
            oof_preds[val_idx, model_idx] = val_proba

            fold_aucs.append(roc_auc_score(y_val_fold, val_proba))

        cv_mean = np.mean(fold_aucs)
        cv_std  = np.std(fold_aucs)
        print(f"    CV AUC: {cv_mean:.4f} (+/- {cv_std:.4f})")

    return oof_preds


# ── Generate OOF predictions (meta-features for meta-learner training) ────────
print("=== Generating OOF predictions ===")
oof_preds = generate_oof_predictions(base_models, X_train, y_train, n_splits=5)

# ── Train meta-learner on OOF predictions ─────────────────────────────────────
meta_learner = LogisticRegression(C=0.1, random_state=42, max_iter=1000)
meta_learner.fit(oof_preds, y_train)

meta_cv = cross_val_score(
    meta_learner, oof_preds, y_train, cv=5, scoring="roc_auc"
)
print(f"\nMeta-learner CV AUC: {meta_cv.mean():.4f} (+/- {meta_cv.std():.4f})")
print(f"Meta-learner weights: {dict(zip(base_models.keys(), meta_learner.coef_[0].round(3)))}")

# ── Generate test predictions: retrain each base model on full training set ───
print("\n=== Retraining base models on full training set ===")
test_meta_features = np.zeros((len(X_test), len(base_models)))
for model_idx, (model_name, model) in enumerate(base_models.items()):
    print(f"  Training: {model_name}")
    model.fit(X_train, y_train)
    test_meta_features[:, model_idx] = model.predict_proba(X_test)[:, 1]

# ── Final stacked prediction ──────────────────────────────────────────────────
stacked_proba = meta_learner.predict_proba(test_meta_features)[:, 1]
stacked_auc   = roc_auc_score(y_test, stacked_proba)

# ── Compare against individual models and simple average ─────────────────────
individual_aucs = {
    name: roc_auc_score(y_test, test_meta_features[:, i])
    for i, name in enumerate(base_models.keys())
}
simple_avg_auc = roc_auc_score(y_test, test_meta_features.mean(axis=1))

print("\n=== Results Summary ===")
for name, auc in individual_aucs.items():
    print(f"  {name:4s}: {auc:.4f}")
print(f"  Simple average: {simple_avg_auc:.4f}")
print(f"  Stacked:        {stacked_auc:.4f}  "
      f"(+{stacked_auc - max(individual_aucs.values()):.4f} vs best individual)")

sklearn StackingClassifier

sklearn's StackingClassifier handles the OOF logic automatically:

from sklearn.ensemble import StackingClassifier
import copy

# Build estimators list (fresh instances - must not be pre-fitted)
estimators = [
    ("xgb", xgb.XGBClassifier(n_estimators=300, max_depth=5, learning_rate=0.05,
                               random_state=42, verbosity=0)),
    ("lgb", lgb.LGBMClassifier(n_estimators=300, max_depth=5, learning_rate=0.05,
                                random_state=42, verbose=-1)),
    ("rf",  RandomForestClassifier(n_estimators=300, max_depth=10, random_state=42, n_jobs=-1)),
    ("et",  ExtraTreesClassifier(n_estimators=300, max_depth=10, random_state=42, n_jobs=-1)),
]

sklearn_stack = StackingClassifier(
    estimators=estimators,
    final_estimator=LogisticRegression(C=0.1, max_iter=1000),
    cv=5,                         # StratifiedKFold by default
    stack_method="predict_proba", # pass probabilities to meta-learner
    n_jobs=1,                     # -1 can cause issues with some estimators
    passthrough=False,            # True = also pass raw features to meta-learner
)
sklearn_stack.fit(X_train, y_train)
sklearn_auc = roc_auc_score(y_test, sklearn_stack.predict_proba(X_test)[:, 1])
print(f"sklearn StackingClassifier AUC: {sklearn_auc:.4f}")

:::tip The passthrough trick Setting passthrough=True gives the meta-learner access to the original features alongside the base model predictions. This sometimes helps when base models have similar predictions but disagree in specific feature regions - the meta-learner can use raw features to resolve disagreements. Try it when your stacked ensemble is not improving over simple averaging. :::

Blending: The Simpler Alternative

Blending replaces cross-validation OOF with a single holdout split:

Blending procedure:
─────────────────────────────────────────────────────────
Training data  →  70-80% for base model training
               →  20-30% holdout (blend set)

1. Train each base model on the 70-80% split only
2. Predict holdout set → meta-features: (holdout_n, M)
3. Train meta-learner on (meta-features, holdout_y)
4. At test time: base models predict → meta-learner predicts final

Note: base models see less training data than in stacking.
Meta-learner trains on fewer samples than in stacking.

from sklearn.model_selection import train_test_split

# ── Blending implementation ───────────────────────────────────────────────────
X_blend_train, X_blend_val, y_blend_train, y_blend_val = train_test_split(
    X_train, y_train, test_size=0.25, stratify=y_train, random_state=42
)

# Train base models on 75% of training data
blend_test_features = np.zeros((len(X_test), len(base_models)))
blend_val_features  = np.zeros((len(X_blend_val), len(base_models)))

for model_idx, (model_name, model) in enumerate(base_models.items()):
    import copy
    model_fresh = copy.deepcopy(model)
    model_fresh.fit(X_blend_train, y_blend_train)
    blend_val_features[:, model_idx]  = model_fresh.predict_proba(X_blend_val)[:, 1]
    blend_test_features[:, model_idx] = model_fresh.predict_proba(X_test)[:, 1]

# Train meta-learner on holdout predictions
blend_meta = LogisticRegression(C=0.1, random_state=42, max_iter=1000)
blend_meta.fit(blend_val_features, y_blend_val)

blend_auc = roc_auc_score(y_test, blend_meta.predict_proba(blend_test_features)[:, 1])
print(f"Blending AUC: {blend_auc:.4f}  (vs Stacking AUC: {stacked_auc:.4f})")

Blending vs Stacking:

Aspect	Blending	Stacking
Implementation complexity	Simple	Moderate
Data efficiency	Wastes holdout split for base model training	Uses all training data for OOF
Meta-learner training data	20-25% of train (smaller)	100% of train via OOF (larger)
Risk of leakage	Lower (single holdout, simpler workflow)	Negligible if OOF implemented correctly
Typical AUC vs stacking	0.001-0.003 AUC lower	Higher (more data)
Best use case	Quick prototyping, very large datasets	Competition ML, careful production builds

Snapshot Ensembling: Multiple Checkpoints from One Training Run

Snapshot ensembling (Huang et al., 2017) generates diverse ensemble members from a single training run by saving model checkpoints at different points in the learning rate cycle. Because different checkpoints correspond to different local minima, they make complementary errors.

import torch
import torch.nn as nn
import numpy as np
from torch.optim.lr_scheduler import CyclicLR

class SimpleNet(nn.Module):
    def __init__(self, input_dim: int, hidden_dims: list, output_dim: int):
        super().__init__()
        layers = []
        prev_dim = input_dim
        for h in hidden_dims:
            layers.extend([nn.Linear(prev_dim, h), nn.ReLU(), nn.Dropout(0.3)])
            prev_dim = h
        layers.append(nn.Linear(prev_dim, output_dim))
        self.net = nn.Sequential(*layers)

    def forward(self, x):
        return self.net(x)


def snapshot_ensemble_train(
    X_train: np.ndarray,
    y_train: np.ndarray,
    n_cycles: int = 5,
    epochs_per_cycle: int = 20,
    lr_max: float = 0.1,
    lr_min: float = 0.001,
) -> list:
    """
    Train with cyclic LR - save snapshot at end of each cycle.
    Returns list of trained models (snapshots).
    """
    import torch
    from torch.utils.data import TensorDataset, DataLoader

    device = "cuda" if torch.cuda.is_available() else "cpu"
    X_t = torch.FloatTensor(X_train).to(device)
    y_t = torch.LongTensor(y_train).to(device)
    dataset = TensorDataset(X_t, y_t)
    loader = DataLoader(dataset, batch_size=256, shuffle=True)

    model = SimpleNet(X_train.shape[1], [256, 128, 64], 2).to(device)
    optimizer = torch.optim.SGD(model.parameters(), lr=lr_max, momentum=0.9)
    scheduler = CyclicLR(
        optimizer, base_lr=lr_min, max_lr=lr_max,
        step_size_up=len(loader) * epochs_per_cycle // 2,
        cycle_momentum=True,
    )
    criterion = nn.CrossEntropyLoss()

    snapshots = []
    epoch = 0
    for cycle in range(n_cycles):
        for _ in range(epochs_per_cycle):
            model.train()
            for X_batch, y_batch in loader:
                optimizer.zero_grad()
                loss = criterion(model(X_batch), y_batch)
                loss.backward()
                optimizer.step()
                scheduler.step()
            epoch += 1

        # Save snapshot at end of each cycle (model is at a local minimum)
        import copy
        snapshots.append(copy.deepcopy(model))
        print(f"Cycle {cycle+1}/{n_cycles}: snapshot saved (epoch {epoch})")

    return snapshots


def snapshot_predict_proba(snapshots: list, X: np.ndarray) -> np.ndarray:
    """Average predictions from all snapshots."""
    device = "cuda" if torch.cuda.is_available() else "cpu"
    X_t = torch.FloatTensor(X).to(device)
    all_probs = []
    for model in snapshots:
        model.eval()
        with torch.no_grad():
            logits = model(X_t)
            probs = torch.softmax(logits, dim=1)[:, 1].cpu().numpy()
        all_probs.append(probs)
    return np.column_stack(all_probs).mean(axis=1)  # average across snapshots

Temporal Ensembling for Sequential Data

Temporal ensembling maintains exponential moving averages of model predictions as a form of ensemble that leverages predictions from past epochs.

class TemporalEnsembleTrainer:
    """
    Temporal ensembling: maintain an EMA of predictions across epochs.
    Useful for semi-supervised learning and consistency regularization.
    """

    def __init__(self, alpha: float = 0.6, n_samples: int = 1000):
        """
        alpha: EMA coefficient (higher = more weight on historical predictions)
        """
        self.alpha = alpha
        self.ensemble_preds = np.zeros(n_samples)  # Z in original paper

    def update_ensemble(self, current_preds: np.ndarray, epoch: int) -> np.ndarray:
        """
        Update EMA ensemble predictions.
        Z_t = alpha * Z_{t-1} + (1 - alpha) * current_preds
        """
        self.ensemble_preds = self.alpha * self.ensemble_preds + (1 - self.alpha) * current_preds

        # Bias correction: account for startup with zero initialization
        bias_correction = 1.0 - self.alpha ** (epoch + 1)
        return self.ensemble_preds / bias_correction

Meta-Feature Engineering

The meta-learner's input can be enriched beyond just base model probabilities:

def build_augmented_meta_features(
    oof_preds: np.ndarray,
    model_names: list,
) -> pd.DataFrame:
    """
    Augment OOF predictions with derived meta-features that help
    the meta-learner understand ensemble behavior.
    """
    meta_df = pd.DataFrame(oof_preds, columns=model_names)

    # Hard predictions from each model
    for name in model_names:
        meta_df[f"{name}_hard"] = (meta_df[name] > 0.5).astype(int)

    # Inter-model disagreement: high disagreement = hard examples
    meta_df["variance"] = oof_preds.var(axis=1)
    meta_df["range"]    = oof_preds.max(axis=1) - oof_preds.min(axis=1)
    meta_df["std_dev"]  = oof_preds.std(axis=1)

    # Model agreement count (how many predict positive)
    meta_df["n_positive_votes"] = (oof_preds > 0.5).sum(axis=1)

    # Rank-based features (position in probability ranking)
    for name in model_names:
        meta_df[f"{name}_rank"] = meta_df[name].rank(pct=True)

    return meta_df


augmented_oof = build_augmented_meta_features(oof_preds, list(base_models.keys()))
print(f"OOF meta-features: {oof_preds.shape[1]} → {augmented_oof.shape[1]} columns")

# Train meta-learner on augmented features
from sklearn.linear_model import RidgeClassifier
meta_ridge = LogisticRegression(C=0.05, random_state=42, max_iter=2000)
meta_ridge.fit(augmented_oof, y_train)

# Generate augmented test meta-features
augmented_test = build_augmented_meta_features(test_meta_features, list(base_models.keys()))
augmented_auc = roc_auc_score(y_test, meta_ridge.predict_proba(augmented_test)[:, 1])
print(f"Augmented meta-features AUC: {augmented_auc:.4f}")

Model Calibration Before Stacking

A critical detail: if base models output raw log-odds or poorly calibrated probabilities, the meta-learner may weight them incorrectly. Always calibrate base models before stacking.

from sklearn.calibration import CalibratedClassifierCV

# Check calibration of each base model
from sklearn.calibration import calibration_curve
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, len(base_models), figsize=(5*len(base_models), 4))
for ax, (name, probs) in zip(axes, zip(base_models.keys(),
                                         [test_meta_features[:, i]
                                          for i in range(len(base_models))])):
    prob_true, prob_pred = calibration_curve(y_test, probs, n_bins=10)
    ax.plot(prob_pred, prob_true, "o-", label=name, color="#2563eb")
    ax.plot([0, 1], [0, 1], "k--", alpha=0.5, label="Perfect")
    ax.set_title(f"{name} Calibration")
    ax.set_xlabel("Mean predicted probability")
    ax.set_ylabel("Fraction of positives")
    ax.legend(fontsize=8)

plt.suptitle("Base Model Calibration (for stacking quality assessment)")
plt.tight_layout()
plt.savefig("calibration_check.png", dpi=150)
plt.close()

# Calibrate an individual base model using isotonic regression
rf_calibrated = CalibratedClassifierCV(
    RandomForestClassifier(n_estimators=200, random_state=42),
    cv=5, method="isotonic",
)
rf_calibrated.fit(X_train, y_train)

Diversity Strategy for Effective Stacking

GOOD diversity sources:
────────────────────────────────────────────────────────────────
✓ Algorithm diversity: trees + linear + neural + kernel
✓ Feature diversity: different feature subsets or transformations
✓ Hyperparameter diversity: deep trees vs shallow trees
✓ Data diversity: different bootstrap samples or augmentations
✓ Training objective diversity: different loss functions
✓ Architecture diversity: tree depth, hidden layers, regularization

BAD diversity (doesn't help much):
────────────────────────────────────────────────────────────────
✗ 5 XGBoost models with slightly different hyperparameters
✗ 3 LightGBM models with different seeds (correlation ρ ≈ 0.97)
✗ Multiple versions of the same feature engineering pipeline
✗ Ensembling models that all predict very high or very low probabilities

When to Stop Adding Models

def evaluate_incremental_stack_value(
    X_test: pd.DataFrame,
    y_test: np.ndarray,
    existing_meta_features: np.ndarray,
    new_model_proba: np.ndarray,
    existing_meta: LogisticRegression,
    model_names: list,
) -> dict:
    """
    Compute the AUC gain from adding a new model to the stack.
    If gain < 0.001, the model is not worth the added complexity.
    """
    # Stack without new model
    baseline_auc = roc_auc_score(
        y_test,
        existing_meta.predict_proba(existing_meta_features)[:, 1]
    )

    # Stack with new model
    augmented_features = np.column_stack([existing_meta_features, new_model_proba])
    new_meta = LogisticRegression(C=0.1, max_iter=1000, random_state=42)
    # Note: in practice, regenerate OOF - this is a simplified check
    new_meta.fit(augmented_features, y_test)  # simplified for demo
    new_auc = roc_auc_score(y_test, new_meta.predict_proba(augmented_features)[:, 1])

    gain = new_auc - baseline_auc
    verdict = "ADD" if gain > 0.001 else "SKIP"
    print(f"Incremental gain from new model: {gain:+.4f} → {verdict}")
    return {"baseline_auc": baseline_auc, "new_auc": new_auc, "gain": gain}

Production Retraining Pipeline

def retrain_stack_pipeline(
    X_new: pd.DataFrame,
    y_new: np.ndarray,
    base_models: dict,
    n_splits: int = 5,
) -> tuple:
    """
    Full stack retrain pipeline.
    Run on a schedule or when drift is detected.
    Cost: O(K * M) + M model trains (K=folds, M=base models)
    """
    print("=== Stack Retraining Pipeline ===")

    # Step 1: Generate fresh OOF predictions with new data
    print("Step 1: Generating OOF predictions...")
    new_oof = generate_oof_predictions(base_models, X_new, y_new, n_splits=n_splits)

    # Step 2: Retrain meta-learner on new OOF
    print("Step 2: Retraining meta-learner...")
    meta_learner = LogisticRegression(C=0.1, max_iter=1000, random_state=42)
    meta_learner.fit(new_oof, y_new)

    # Step 3: Retrain base models on full new dataset for inference
    print("Step 3: Retraining base models for inference...")
    for name, model in base_models.items():
        model.fit(X_new, y_new)
        print(f"  Retrained: {name}")

    # Cost summary
    total_trains = n_splits * len(base_models) + len(base_models)
    print(f"\nTotal model trains: {total_trains} ({n_splits}K OOF + {len(base_models)} final)")

    return base_models, meta_learner

Comparison: Bagging vs Boosting vs Stacking

Dimension	Voting	Bagging	Boosting	Stacking
Algorithm diversity	Yes	No	No	Yes
Combination method	Fixed weights	Average	Weighted sum	Trained meta-model
Reduces bias	No	No	Yes	Yes (via diversity)
Reduces variance	Partly	Yes	No	Yes
Training complexity	Low	Medium	Medium	High
Inference complexity	Medium	Medium	Medium	High
Risk of overfitting	Low	Low	Medium	Medium-High (if OOF violated)
Typical AUC gain vs best single model	+0.001-0.003	+0.003-0.007	+0.005-0.015	+0.003-0.010
When to use	Quick win	High variance	High bias	Maximum performance

Production Engineering Notes

Latency

A stack of five models runs inference on all five in sequence (or parallel). Five models at 5 ms each = 5 ms parallel vs 25 ms sequential. In real-time systems, model parallelism is essential:

import concurrent.futures

def parallel_base_predictions(base_models_fitted, X_input):
    """Run all base model predictions in parallel."""
    with concurrent.futures.ThreadPoolExecutor(max_workers=len(base_models_fitted)) as executor:
        futures = {
            name: executor.submit(model.predict_proba, X_input)
            for name, model in base_models_fitted.items()
        }
        results = {name: fut.result()[:, 1] for name, fut in futures.items()}
    return np.column_stack(list(results.values()))


# Benchmark sequential vs parallel
import time
X_batch = X_test.iloc[:100]

start = time.perf_counter()
sequential_preds = np.column_stack([m.predict_proba(X_batch)[:, 1] for m in base_models.values()])
sequential_ms = (time.perf_counter() - start) * 1000

start = time.perf_counter()
parallel_preds = parallel_base_predictions(base_models, X_batch)
parallel_ms = (time.perf_counter() - start) * 1000

print(f"Sequential: {sequential_ms:.1f}ms  Parallel: {parallel_ms:.1f}ms  "
      f"Speedup: {sequential_ms/parallel_ms:.1f}x")

Maintenance Surface

Six models (5 base + 1 meta) means six models to retrain, monitor, and version. When one base model drifts, retraining it invalidates the meta-learner's OOF predictions - you must regenerate OOF and retrain the meta-learner as well. Factor this into your retraining schedule.

When stacking is worth it in production:

• High-stakes domain where every 0.001 AUC matters (credit, healthcare, fraud)
• Batch prediction pipeline where latency is not a constraint
• The performance gap exceeds the SLA cost of added latency

When a single model is better in production:

• Real-time serving with tight SLA (less than 20 ms)
• Continuous or online learning pipelines (stacks are hard to update incrementally)
• Small team - six models means six failure modes to debug at 2 AM

:::danger Never skip OOF in production stacking The most common stacking implementation mistake is generating base model predictions on the same data they were trained on, then training the meta-learner on those overfit predictions. The resulting stack looks fantastic in training evaluation and collapses at test time. The fix is non-negotiable: always use K-fold OOF predictions for meta-learner training. There is no shortcut. If you cannot afford the compute for full K-fold OOF, use blending (single holdout) instead - it is honest, if less data-efficient. :::

YouTube Resources

Video	Channel	Why Watch It
Stacking Ensemble - Complete Tutorial	StatQuest with Josh Starmer	Best conceptual walkthrough with OOF predictions
Ensemble Methods in Machine Learning	Krish Naik	All ensemble types with sklearn implementations
Kaggle Competition Ensembling Tricks	Abhishek Thakur	Competition-proven ensembling strategies
sklearn StackingClassifier Tutorial	Data Science Dojo	End-to-end sklearn stacking with real dataset
Snapshot Ensembling Paper Explained	Yannic Kilcher	Snapshot ensembling and cyclic learning rates

Interview Questions and Answers

Q1: Why must stacking use out-of-fold predictions rather than base model predictions on the training set?

When a base model predicts on the same data it was trained on, it produces optimistic, overfit predictions - near-perfect for training examples due to memorization. If the meta-learner trains on these overfit predictions, it learns to trust predictions that are unrealistically good. At test time, base models produce normal (more realistic, worse) predictions. The meta-learner, calibrated to overfit predictions, is now miscalibrated and performs poorly. Out-of-fold (OOF) predictions solve this by ensuring every training sample's prediction comes from a model that never saw that sample. This means the meta-learner sees realistic base model predictions during training - the same quality it will see at test time. OOF predictions are the stacking equivalent of a held-out validation set, but generalized to use all training data efficiently.

Q2: Explain the difference between hard voting, soft voting, weighted voting, and stacking. Why does each improve on the previous?

Hard voting aggregates class labels by majority vote - ignores prediction confidence, so a model predicting 0.99 and one predicting 0.51 have equal say. Soft voting averages predicted probabilities before thresholding - preserves confidence, giving confident models more effective weight. Weighted voting assigns learned or heuristic weights per model, so stronger models have more influence - better than equal weighting when models differ in quality. Stacking trains a meta-learner on base model predictions - instead of fixed weights, it learns sample-adaptive combination. For a specific input where the linear model is reliably better than the tree model, stacking can down-weight the tree model's prediction for that type of input. This is something fixed-weight methods cannot do. Stacking's meta-learner discovers which models are more reliable for which regions of the input space.

Q3: What is model diversity and why does it matter more than individual model accuracy for ensembles?

Adding a sixth XGBoost variant with slightly different hyperparameters to a stack that already has five XGBoost variants provides near-zero benefit - the six models are highly correlated (ρ ≈ 0.97), so the variance formula $\text{Var}(\bar{X}) = \rho\sigma^2 + (1-\rho)\sigma^2/B$ shows that the variance floor $\rho\sigma^2$ remains nearly unchanged. Adding a fundamentally different model - a linear model, a k-NN, or a neural network - introduces complementary errors (low ρ), which the variance formula translates directly into variance reduction. In practice: an ensemble of one XGBoost (0.90 AUC), one logistic regression (0.83 AUC), and one MLP (0.86 AUC) often outperforms an ensemble of three slightly different XGBoost models all at 0.90 AUC. Diversity matters more than raw accuracy of individual models.

Q4: What is snapshot ensembling and how does it provide diversity from a single training run?

Snapshot ensembling (Huang et al., 2017) uses cyclical learning rates to drive the model to different local minima during training. In each cycle, the learning rate rises from a low value to a high value, allowing the optimizer to escape local minima, then falls back down, converging to a new local minimum. Saving the model at the end of each cycle captures a model at a different local minimum - these different convergence points represent different learned representations with complementary strengths. The ensemble is formed by averaging predictions from all snapshots. The key insight: local minima in neural network loss surfaces tend to have similar loss values but different functional behavior on specific examples - exactly the complementarity needed for effective ensembling. Snapshot ensembling gets diversity at the cost of training time (M cycles instead of 1) rather than at the cost of model storage or inference time.

Q5: You have a stacked ensemble with AUC 0.93 on CV but 0.88 on the test set. What went wrong and how do you fix it?

A 0.05 AUC gap between CV and test strongly suggests data leakage in the stacking procedure. The most likely cause: base model predictions were generated on training data without OOF, allowing the meta-learner to learn overfit predictions. Fix: (1) verify that OOF is implemented correctly - every training sample's prediction must come from a model trained without that sample; (2) check for target leakage in features - any feature computed using future information or the target variable itself; (3) check whether the CV folds match the data's temporal structure - for time series data, standard K-fold leaks future into past; use time-based splits instead; (4) check if feature engineering was done on the full training set before splitting (e.g., mean-encoding with global statistics) - this leaks target information into features; feature engineering must be done within each fold. The 0.05 gap is too large to be random variation - it almost certainly indicates structural leakage somewhere in the pipeline.

Q6: In a production system, when would you use blending instead of stacking?

Blending uses a single holdout split instead of K-fold OOF for generating meta-learner training data. Use blending when: (1) dataset is very large (50M+ rows) where K-fold OOF with 5 base models means 5×5=25 training runs - computationally prohibitive; (2) time constraints prevent full K-fold iteration; (3) data has a temporal structure where K-fold leaks future information - blending with a time-based holdout is cleaner; (4) quick prototyping to validate whether ensembling will help before investing in full stacking infrastructure. The cost: base models see 20-30% less training data than in stacking (since holdout is excluded), and the meta-learner trains on 20-30% of samples instead of 100% via OOF. This typically costs 0.001-0.003 AUC versus full stacking on medium-sized datasets. For very large datasets, the data efficiency advantage of stacking is negligible, making blending the practical choice.

Incremental Stack Value: Know When to Stop Adding Models

A key production skill is knowing when adding another base model to a stack stops paying off:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.linear_model import LogisticRegression, RidgeClassifier
from sklearn.svm import SVC
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.metrics import roc_auc_score
from sklearn.base import clone

def evaluate_stack_value(X: np.ndarray, y: np.ndarray,
                           base_models: list, meta_model,
                           model_names: list,
                           cv_folds: int = 5):
    """
    Evaluate the marginal value of adding each model to the stack.
    Reports AUC at each step: single model → 2-model stack → ... → full stack.
    """
    cv = StratifiedKFold(n_splits=cv_folds, shuffle=True, random_state=42)

    # Step 1: Generate OOF predictions for all base models
    print("Generating OOF predictions...")
    oof_predictions = np.zeros((len(y), len(base_models)))

    for m_idx, model in enumerate(base_models):
        for fold, (train_idx, val_idx) in enumerate(cv.split(X, y)):
            m_clone = clone(model)
            m_clone.fit(X[train_idx], y[train_idx])
            if hasattr(m_clone, 'predict_proba'):
                oof_predictions[val_idx, m_idx] = m_clone.predict_proba(X[val_idx])[:, 1]
            else:
                oof_predictions[val_idx, m_idx] = m_clone.decision_function(X[val_idx])
        print(f"  {model_names[m_idx]}: done (AUC={roc_auc_score(y, oof_predictions[:,m_idx]):.4f})")

    # Step 2: Evaluate meta-learner at each stack size
    results = []
    for k in range(1, len(base_models) + 1):
        meta_preds = []
        for train_idx, val_idx in cv.split(X, y):
            meta_train = oof_predictions[train_idx, :k]
            meta_val   = oof_predictions[val_idx, :k]

            meta_clone = clone(meta_model)
            meta_clone.fit(meta_train, y[train_idx])

            if hasattr(meta_clone, 'predict_proba'):
                preds = meta_clone.predict_proba(meta_val)[:, 1]
            else:
                preds = meta_clone.decision_function(meta_val)

            meta_preds.extend(zip(val_idx, preds))

        meta_preds.sort(key=lambda x: x[0])
        final_preds = np.array([p for _, p in meta_preds])
        auc = roc_auc_score(y, final_preds)
        results.append({'k': k, 'models': model_names[:k], 'auc': auc})
        print(f"  Stack ({'+'.join(model_names[:k])}): AUC={auc:.4f}")

    # Plot marginal value curve
    ks   = [r['k'] for r in results]
    aucs = [r['auc'] for r in results]

    plt.figure(figsize=(10, 5))
    plt.plot(ks, aucs, 'o-', linewidth=2, markersize=8, color='#3b82f6')
    plt.xticks(ks, ['\n'.join(model_names[:k]) for k in ks], fontsize=8)
    plt.xlabel("Stack composition (models added left to right)")
    plt.ylabel("Cross-validated AUC")
    plt.title("Incremental Stack Value\n"
               "(diminishing returns → stop when gain < 0.001)")
    plt.grid(True, alpha=0.3)

    # Annotate marginal gains
    for i in range(1, len(aucs)):
        gain = aucs[i] - aucs[i-1]
        plt.annotate(f"+{gain:.4f}", xy=(ks[i], aucs[i]), xytext=(0, 10),
                     textcoords='offset points', ha='center', fontsize=8,
                     color='#16a34a' if gain > 0.001 else '#dc2626')

    plt.tight_layout()
    plt.show()

    return results


# Usage
from sklearn.datasets import make_classification

X, y = make_classification(n_samples=5000, n_features=20,
                             n_informative=12, n_redundant=4, random_state=42)

base_models = [
    RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1),
    GradientBoostingClassifier(n_estimators=200, learning_rate=0.05, random_state=42),
    LogisticRegression(C=1.0, max_iter=500, random_state=42),
    SVC(probability=True, kernel='rbf', C=1.0, random_state=42),
]
model_names = ['RF', 'GBM', 'LR', 'SVM']
meta_model  = LogisticRegression(C=0.1, max_iter=500, random_state=42)

results = evaluate_stack_value(X, y, base_models, meta_model, model_names)

## Common Stacking Mistakes

```python
# ── MISTAKE 1: Training base models on full training set ──────────────────────
# WRONG: base models trained on X_train, then predict on X_train for meta-learner
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
import numpy as np

rf = RandomForestClassifier().fit(X_train, y_train)
# This is WRONG - RF has seen all of X_train - it overfits to produce near-perfect
# in-sample predictions that fool the meta-learner:
meta_train_features_WRONG = rf.predict_proba(X_train)[:, 1]

# CORRECT: OOF predictions - RF was trained without the predicted samples
from sklearn.model_selection import StratifiedKFold, cross_val_predict
meta_train_features_CORRECT = cross_val_predict(rf, X_train, y_train,
                                                  cv=5, method='predict_proba')[:, 1]
# Now each sample's prediction comes from a model that didn't see it

# ── MISTAKE 2: Leaking test set into blending holdout ─────────────────────────
# WRONG: using the test set as the blending holdout
# This evaluates on the same data used to train the meta-learner → optimistic AUC

# CORRECT: 3-way split for blending
# train_blend → train base models
# holdout     → generate predictions to train meta-learner (blending)
# test        → final evaluation only (never used until the very end)

# ── MISTAKE 3: Not including original features in meta-learner input ───────────
# SOMETIMES WRONG: meta-learner only gets base model predictions
# Often BETTER: also include original features (or SHAP values from base models)
# because the meta-learner can correct systematic biases that all base models share

# Extended meta-features: base model predictions + original features
meta_features_extended = np.hstack([
    meta_train_features_CORRECT.reshape(-1, 1),
    X_train  # original features - gives meta-learner more context
])

:::tip 🎮 Interactive Playground

Visualize this concept: Try the Ensemble Methods demo on the EngineersOfAI Playground - no code required.

:::