Keeping an XGBoost model honest: drift monitoring, champion/challenger, and retraining cadence

LEARN · GRADIENT BOOSTING WITH XGBOOST

An XGBoost model can remain perfectly healthy from an infrastructure perspective while becoming operationally wrong. The endpoint still returns 200 OK. Latency stays flat. CPU and memory look normal. The model file has not changed. Yet its inputs no longer mean what they meant during training, its predicted-score distribution has shifted, or the relationship between features and outcomes has decayed.

Keeping a model honest therefore requires more than periodic retraining. It requires three connected control loops:

Drift monitoring that watches inputs, predictions, data quality, and delayed performance.
A champion/challenger process that compares models on the same production traffic before promotion.
A retraining cadence derived from measured decay, label delay, and deployment lead time rather than an arbitrary monthly or quarterly schedule.

The order matters. When an alert fires, investigate the data pipeline first. Retrain only after establishing that the data is valid and the learned relationship has genuinely changed.

The four questions a production monitor must answer

“Model drift” is often used as a catch-all term, but production diagnosis becomes easier when the problem is split into four questions.

Did the input distribution change?

This is feature drift, sometimes called covariate drift. Examples include:

Order values becoming larger.
A new customer segment entering the system.
A categorical feature gaining previously unseen values.
A feature’s missing rate increasing.
A timestamp or currency conversion changing upstream.

Feature drift is observable without labels, which makes it useful when outcomes arrive days or weeks later.

It is not automatically evidence that the model is worse. A changed feature may be irrelevant to the model, or the model may still generalize well in the shifted region.

Did the model’s output distribution change?

Prediction drift asks whether predicted probabilities, margins, classes, rankings, or numerical outputs have moved.

This deserves its own monitor because a model is a nonlinear compression of the entire feature vector. Ten individual feature distributions can each look stable while correlations among them change. Those joint changes can move model scores even though every one-dimensional feature monitor remains green.

Prediction drift is an early-warning signal, not proof of performance loss. It can reveal:

Joint feature drift.
Changed feature interactions.
A preprocessing or model-version mismatch.
A shift in the proportion of traffic reaching the model.
A change in missing-value defaults.
A genuine change in the population being scored.

Did measured performance change?

Once labels arrive, calculate the metrics that describe the actual service objective.

For a binary classifier, that might include:

ROC AUC for ranking quality.
Log loss or Brier score for probability quality.
Recall at a fixed review capacity.
Precision at a business threshold.
Capture rate in the highest-scored 5% or 10%.
Metrics split by region, channel, customer type, or product family.

Current scikit-learn documentation exposes roc_auc_score for ranking performance and log_loss for probability loss. The latter should be tracked whenever predicted probabilities drive pricing, prioritization, expected-value calculations, or capacity planning—not merely a final class label.

Did the relationship between features and labels change?

Concept drift means that similar inputs now imply different outcomes.

For example, a return-risk model might have learned that high discounts on expensive orders predict returns. A new returns policy, merchandising strategy, or customer incentive can weaken or reverse that relationship.

Unlabeled feature monitoring cannot prove concept drift. The input distribution could remain almost unchanged while the conditional relationship between inputs and labels changes substantially. Delayed-label performance monitoring is therefore not optional.

Build a trustworthy reference window

A drift metric is only as meaningful as its reference data.

Do not automatically compare production traffic with the entire training set. The training set may contain oversampling, class weights, synthetic rows, historical periods, or preprocessing artifacts that never represented real serving traffic.

A good reference window is usually:

Drawn from a period when the model was known to perform acceptably.
Processed through the exact production feature pipeline.
Large enough to represent normal segments and seasonality.
Stored immutably with its schema and feature definitions.
Associated with a specific model version.
Separate from hyperparameter-tuning data where possible.

For seasonal systems, one global reference may be insufficient. Monday morning traffic should often be compared with previous Monday morning traffic, not Saturday night traffic.

Store at least:

Feature values or approved distribution summaries.
Prediction scores.
Missing rates.
Category frequencies.
Model version.
Feature-schema version.
Transformation-code version.
Event timestamps.
Segment keys.
Reference-window start and end dates.

The reference is part of the deployed model artifact. Changing it silently is equivalent to changing an alert definition silently.

Use current, supported APIs

As of August 2026, XGBoost’s stable release notes list version 3.3.0, released on June 17, 2026. Its supported Python interface includes XGBClassifier, and the official documentation supports saving models in JSON format.

SciPy’s current statistics API includes scipy.stats.ks_2samp for the two-sample Kolmogorov–Smirnov test.

Create a small project:

honest-xgb/
├── drift_metrics.py
├── run_demo.py
└── artifacts/

Install the dependencies:

mkdir honest-xgb
cd honest-xgb
python -m venv .venv
source .venv/bin/activate
python -m pip install --upgrade pip
python -m pip install "xgboost==3.3.0" "scikit-learn==1.9.0" "scipy==1.18.0" pandas numpy
mkdir artifacts

On Windows PowerShell, activate the environment with:

.venv\Scripts\Activate.ps1

Population Stability Index: useful, but easy to misuse

Population Stability Index, or PSI, compares how much probability mass falls into fixed buckets in the reference and current samples.

In plain-text notation:

PSI = Σ (current share − reference share) × ln(current share / reference share)

A production-quality implementation needs several details that short examples often omit:

Derive numeric bucket boundaries from the reference only.
Reuse those exact boundaries for every live window.
Include underflow and overflow values.
Include missing values as an explicit bucket.
Prevent division by zero with a small epsilon.
Treat low-cardinality features as categorical rather than forcing quantile bins.
Keep sample size and window definition consistent.

If bucket boundaries are recalculated on current data, each dataset can be made to look artificially similar because both are divided into their own quantiles.

Practical PSI alert bands

There is no universal PSI law. The following bands are useful starting points, not scientific constants:

PSI below 0.10: usually ordinary movement.
PSI from 0.10 to 0.20: warning and investigation.
PSI above 0.20: strong shift.
Some teams use 0.25 rather than 0.20 as the critical threshold.

Prediction-score PSI often deserves tighter starting thresholds, such as:

Warning at 0.05.
Critical at 0.10.

That is because the score combines many features and interactions. A moderate score shift can affect approval rates, review queues, inventory allocation, or customer treatment even when individual inputs have not crossed their thresholds.

Do not copy these values directly into a pager. Backtest them against months of known-stable historical windows, known incidents, seasonal events, and benign campaigns. Choose thresholds that control both missed incidents and alert fatigue.

Kolmogorov–Smirnov: alert on effect size, not just the p-value

For continuous values, the two-sample KS statistic is:

KS = max |F_reference(x) − F_current(x)|

It measures the largest distance between the two empirical cumulative-distribution functions.

SciPy returns both a statistic and a p-value. In monitoring, the statistic is usually the more operationally useful quantity.

With hundreds of thousands of observations, a tiny and harmless shift can produce a very small p-value. A rule such as “page whenever p < 0.05” will eventually page for almost everything.

A more useful starting policy is:

Require a minimum sample size.
Alert on a meaningful KS statistic, such as 0.10 or 0.20.
Record the p-value as supporting context.
Require persistence across multiple windows.
Compare equivalent time periods and segments.
Calibrate feature-specific thresholds from stable history.

KS is designed for continuous distributions. For categorical variables, monitor category frequencies with categorical PSI, total-variation distance, Jensen–Shannon distance, or an appropriately designed chi-squared procedure.

Create drift_metrics.py:

from __future__ import annotations

from typing import Any

import numpy as np
import pandas as pd
from scipy.stats import ks_2samp


def _numeric_edges(reference: np.ndarray, bins: int) -> np.ndarray:
    finite = reference[np.isfinite(reference)]
    if finite.size < 2:
        raise ValueError("Reference data needs at least two finite values.")

    unique = np.unique(finite)

    if unique.size == 1:
        value = float(unique[0])
        delta = max(abs(value), 1.0) * 1e-6
        return np.array(
            [-np.inf, value - delta, value + delta, np.inf]
        )

    quantiles = np.linspace(0.0, 1.0, bins + 1)
    boundaries = np.unique(np.quantile(finite, quantiles))
    return np.concatenate(
        ([-np.inf], boundaries[1:-1], [np.inf])
    )


def _numeric_psi(
    reference: np.ndarray,
    current: np.ndarray,
    bins: int,
    epsilon: float,
) -> float:
    edges = _numeric_edges(reference, bins)

    reference_counts, _ = np.histogram(
        reference[np.isfinite(reference)],
        bins=edges,
    )
    current_counts, _ = np.histogram(
        current[np.isfinite(current)],
        bins=edges,
    )

    reference_counts = np.append(
        reference_counts,
        np.count_nonzero(~np.isfinite(reference)),
    )
    current_counts = np.append(
        current_counts,
        np.count_nonzero(~np.isfinite(current)),
    )

    reference_share = np.clip(
        reference_counts / len(reference),
        epsilon,
        None,
    )
    current_share = np.clip(
        current_counts / len(current),
        epsilon,
        None,
    )

    reference_share /= reference_share.sum()
    current_share /= current_share.sum()

    value = np.sum(
        (current_share - reference_share)
        * np.log(current_share / reference_share)
    )
    return float(value)


def _categorical_psi(
    reference: pd.Series,
    current: pd.Series,
    epsilon: float,
) -> float:
    reference_values = reference.astype("object").where(
        reference.notna(),
        "__MISSING__",
    )
    current_values = current.astype("object").where(
        current.notna(),
        "__MISSING__",
    )

    categories = reference_values.unique().tolist()

    for category in current_values.unique():
        if category not in categories:
            categories.append(category)

    reference_share = (
        reference_values.value_counts(normalize=True)
        .reindex(categories, fill_value=0.0)
        .to_numpy(dtype=float)
    )
    current_share = (
        current_values.value_counts(normalize=True)
        .reindex(categories, fill_value=0.0)
        .to_numpy(dtype=float)
    )

    reference_share = np.clip(reference_share, epsilon, None)
    current_share = np.clip(current_share, epsilon, None)
    reference_share /= reference_share.sum()
    current_share /= current_share.sum()

    value = np.sum(
        (current_share - reference_share)
        * np.log(current_share / reference_share)
    )
    return float(value)


def population_stability_index(
    reference: pd.Series | np.ndarray,
    current: pd.Series | np.ndarray,
    *,
    categorical: bool = False,
    bins: int = 10,
    epsilon: float = 1e-6,
) -> float:
    if categorical:
        return _categorical_psi(
            pd.Series(reference),
            pd.Series(current),
            epsilon,
        )

    return _numeric_psi(
        np.asarray(reference, dtype=float),
        np.asarray(current, dtype=float),
        bins,
        epsilon,
    )


def ks_effect(
    reference: pd.Series | np.ndarray,
    current: pd.Series | np.ndarray,
) -> tuple[float, float]:
    reference_array = np.asarray(reference, dtype=float)
    current_array = np.asarray(current, dtype=float)

    reference_array = reference_array[
        np.isfinite(reference_array)
    ]
    current_array = current_array[
        np.isfinite(current_array)
    ]

    result = ks_2samp(
        reference_array,
        current_array,
        alternative="two-sided",
        method="auto",
    )
    return float(result.statistic), float(result.pvalue)


def alert_level(
    *,
    psi: float,
    ks: float | None,
    prediction: bool,
) -> str:
    ks_value = 0.0 if ks is None else ks

    if prediction:
        if psi >= 0.10 or ks_value >= 0.10:
            return "critical"
        if psi >= 0.05 or ks_value >= 0.05:
            return "warning"
        return "ok"

    if psi >= 0.20 or ks_value >= 0.20:
        return "critical"
    if psi >= 0.10 or ks_value >= 0.10:
        return "warning"
    return "ok"


def build_drift_report(
    reference: pd.DataFrame,
    current: pd.DataFrame,
    features: list[str],
    reference_score: np.ndarray,
    current_score: np.ndarray,
) -> pd.DataFrame:
    rows: list[dict[str, Any]] = []

    for feature in features:
        categorical = (
            reference[feature].nunique(dropna=True) <= 20
        )

        psi = population_stability_index(
            reference[feature],
            current[feature],
            categorical=categorical,
        )

        if categorical:
            ks = None
            p_value = None
        else:
            ks, p_value = ks_effect(
                reference[feature],
                current[feature],
            )

        rows.append(
            {
                "name": feature,
                "kind": "feature",
                "psi": psi,
                "ks": ks,
                "ks_p_value": p_value,
                "reference_missing_rate": float(
                    reference[feature].isna().mean()
                ),
                "current_missing_rate": float(
                    current[feature].isna().mean()
                ),
                "status": alert_level(
                    psi=psi,
                    ks=ks,
                    prediction=False,
                ),
            }
        )

    score_psi = population_stability_index(
        reference_score,
        current_score,
        categorical=False,
    )
    score_ks, score_p_value = ks_effect(
        reference_score,
        current_score,
    )

    rows.append(
        {
            "name": "predicted_return_probability",
            "kind": "prediction",
            "psi": score_psi,
            "ks": score_ks,
            "ks_p_value": score_p_value,
            "reference_missing_rate": 0.0,
            "current_missing_rate": 0.0,
            "status": alert_level(
                psi=score_psi,
                ks=score_ks,
                prediction=True,
            ),
        }
    )

    severity = {
        "critical": 0,
        "warning": 1,
        "ok": 2,
    }

    report = pd.DataFrame(rows)
    report["_severity"] = report["status"].map(severity)

    return (
        report.sort_values(
            ["_severity", "kind", "psi"],
            ascending=[True, True, False],
        )
        .drop(columns="_severity")
        .reset_index(drop=True)
    )

A runnable e-commerce example

The following example predicts whether an order will be returned within 30 days. That delay is realistic enough to demonstrate why unlabeled drift signals and delayed performance metrics must coexist.

The demo creates three situations:

A silent seconds-to-milliseconds feature bug.
A joint-distribution change that leaves individual feature marginals stable but moves predictions.
A genuine concept change that justifies training and evaluating a challenger.

Create run_demo.py:

from __future__ import annotations

import json
import os
from datetime import datetime, timedelta, timezone
from pathlib import Path
from typing import Any

import numpy as np
import pandas as pd
from scipy.special import expit
from sklearn.metrics import (
    brier_score_loss,
    log_loss,
    roc_auc_score,
)
from xgboost import XGBClassifier

from drift_metrics import build_drift_report


FEATURES = [
    "order_value",
    "items_in_cart",
    "account_age_days",
    "seconds_since_last_order",
    "discount_fraction",
    "mobile_order",
    "prior_returns_90d",
]

ARTIFACTS = Path("artifacts")
ARTIFACTS.mkdir(exist_ok=True)


def make_orders(
    n: int,
    seed: int,
    *,
    concept_shift: bool = False,
    joint_shift: bool = False,
) -> tuple[pd.DataFrame, np.ndarray]:
    rng = np.random.default_rng(seed)

    order_value = rng.lognormal(
        mean=4.25,
        sigma=0.65,
        size=n,
    )
    items = np.clip(
        np.rint(
            order_value / 35
            + rng.normal(1.5, 1.2, size=n)
        ),
        1,
        20,
    ).astype(int)
    account_age = rng.gamma(
        shape=2.2,
        scale=180.0,
        size=n,
    )
    seconds_since = rng.lognormal(
        mean=7.0,
        sigma=0.9,
        size=n,
    )
    discount = rng.beta(
        2.0,
        8.0,
        size=n,
    )
    mobile = rng.binomial(
        1,
        0.62,
        size=n,
    )
    prior_returns = np.clip(
        rng.poisson(0.7, size=n),
        0,
        8,
    )

    if joint_shift:
        order_rank = np.argsort(order_value)
        sorted_discount = np.sort(discount)[::-1]

        shifted_discount = np.empty_like(discount)
        shifted_discount[order_rank] = sorted_discount
        discount = shifted_discount

    high_value_high_discount = (
        (order_value > 110)
        & (discount > 0.25)
    ).astype(float)

    if not concept_shift:
        logit = (
            -3.7
            + 0.008 * (order_value - 70)
            + 0.16 * (items - 3)
            - 0.0012 * account_age
            - 0.00008 * seconds_since
            + 1.9 * discount
            + 0.45 * mobile
            + 0.62 * prior_returns
            + 1.10 * high_value_high_discount
        )
    else:
        logit = (
            -3.3
            + 0.003 * (order_value - 70)
            + 0.08 * (items - 3)
            - 0.0005 * account_age
            - 0.00002 * seconds_since
            + 0.8 * discount
            + 0.05 * mobile
            + 0.82 * prior_returns
            - 0.35 * high_value_high_discount
        )

    probability = expit(logit)
    labels = rng.binomial(1, probability)

    frame = pd.DataFrame(
        {
            "order_value": order_value,
            "items_in_cart": items,
            "account_age_days": account_age,
            "seconds_since_last_order": seconds_since,
            "discount_fraction": discount,
            "mobile_order": mobile,
            "prior_returns_90d": prior_returns,
        }
    )

    return frame, labels


def train_model(
    X: pd.DataFrame,
    y: np.ndarray,
    *,
    seed: int,
) -> XGBClassifier:
    model = XGBClassifier(
        objective="binary:logistic",
        eval_metric="logloss",
        tree_method="hist",
        n_estimators=300,
        max_depth=4,
        learning_rate=0.05,
        min_child_weight=4,
        subsample=0.85,
        colsample_bytree=0.85,
        reg_lambda=2.0,
        random_state=seed,
        n_jobs=-1,
    )
    model.fit(X[FEATURES], y)
    return model


def performance_metrics(
    y_true: np.ndarray,
    probability: np.ndarray,
) -> dict[str, float]:
    cutoff = np.quantile(probability, 0.90)
    selected = probability >= cutoff
    positive_count = max(int(y_true.sum()), 1)

    return {
        "roc_auc": float(
            roc_auc_score(y_true, probability)
        ),
        "log_loss": float(
            log_loss(y_true, probability)
        ),
        "brier": float(
            brier_score_loss(y_true, probability)
        ),
        "top_decile_capture": float(
            y_true[selected].sum() / positive_count
        ),
        "selected_rate": float(selected.mean()),
    }


def paired_bootstrap_log_loss_delta(
    y_true: np.ndarray,
    champion_probability: np.ndarray,
    challenger_probability: np.ndarray,
    *,
    iterations: int = 1000,
    seed: int = 123,
) -> tuple[float, float, float]:
    epsilon = np.finfo(float).eps

    champion = np.clip(
        champion_probability,
        epsilon,
        1.0 - epsilon,
    )
    challenger = np.clip(
        challenger_probability,
        epsilon,
        1.0 - epsilon,
    )

    champion_loss = -(
        y_true * np.log(champion)
        + (1 - y_true) * np.log(1.0 - champion)
    )
    challenger_loss = -(
        y_true * np.log(challenger)
        + (1 - y_true) * np.log(1.0 - challenger)
    )

    per_row_delta = challenger_loss - champion_loss

    rng = np.random.default_rng(seed)
    bootstrap_means = np.empty(iterations)
    n = len(y_true)

    for iteration in range(iterations):
        indices = rng.integers(0, n, size=n)
        bootstrap_means[iteration] = np.mean(
            per_row_delta[indices]
        )

    point = float(per_row_delta.mean())
    low, high = np.quantile(
        bootstrap_means,
        [0.025, 0.975],
    )

    return point, float(low), float(high)


def atomic_write_json(
    path: Path,
    payload: dict[str, Any],
) -> None:
    temporary = path.with_suffix(
        path.suffix + ".tmp"
    )
    temporary.write_text(
        json.dumps(payload, indent=2),
        encoding="utf-8",
    )
    os.replace(temporary, path)


def forecast_retraining_trigger(
    history: pd.DataFrame,
    *,
    metric_floor: float,
    preparation_days: int,
    safety_buffer_days: int,
) -> dict[str, Any]:
    if len(history) < 6:
        raise ValueError(
            "Use at least six labeled windows."
        )

    ordered = history.sort_values(
        "window_end"
    ).copy()

    origin = ordered["window_end"].iloc[0]
    age_days = (
        ordered["window_end"] - origin
    ).dt.days.to_numpy(dtype=float)
    metric = ordered["roc_auc"].to_numpy(
        dtype=float
    )

    slope, intercept = np.polyfit(
        age_days,
        metric,
        deg=1,
    )

    if slope >= 0:
        return {
            "slope_per_day": float(slope),
            "forecast_breach_date": None,
            "recommended_start_date": None,
            "reason": "No declining trend.",
        }

    breach_age = (
        metric_floor - intercept
    ) / slope

    breach_date = (
        origin
        + pd.to_timedelta(
            breach_age,
            unit="D",
        )
    )

    total_lead = timedelta(
        days=(
            preparation_days
            + safety_buffer_days
        )
    )
    recommended_start = (
        breach_date.to_pydatetime()
        - total_lead
    )

    latest_date = (
        ordered["window_end"]
        .iloc[-1]
        .to_pydatetime()
    )

    return {
        "slope_per_day": float(slope),
        "forecast_breach_date": (
            breach_date.isoformat()
        ),
        "recommended_start_date": (
            recommended_start.isoformat()
        ),
        "already_due": (
            recommended_start <= latest_date
        ),
    }


def print_report(
    title: str,
    report: pd.DataFrame,
) -> None:
    print(f"\n{title}")
    print(
        report.head(5).to_string(
            index=False
        )
    )


def main() -> None:
    train_X, train_y = make_orders(
        40_000,
        seed=1,
    )
    reference_X, reference_y = make_orders(
        15_000,
        seed=2,
    )

    champion = train_model(
        train_X,
        train_y,
        seed=42,
    )
    champion_path = (
        ARTIFACTS / "champion.json"
    )
    champion.save_model(champion_path)

    reference_score = champion.predict_proba(
        reference_X[FEATURES]
    )[:, 1]

    print("\nReference performance")
    print(
        json.dumps(
            performance_metrics(
                reference_y,
                reference_score,
            ),
            indent=2,
        )
    )

    buggy_live_X, _ = make_orders(
        10_000,
        seed=5,
    )
    buggy_live_X[
        "seconds_since_last_order"
    ] *= 1000

    buggy_score = champion.predict_proba(
        buggy_live_X[FEATURES]
    )[:, 1]

    print_report(
        "Silent unit-change incident",
        build_drift_report(
            reference_X,
            buggy_live_X,
            FEATURES,
            reference_score,
            buggy_score,
        ),
    )

    corrected_live_X = buggy_live_X.copy()
    corrected_live_X[
        "seconds_since_last_order"
    ] /= 1000

    corrected_score = champion.predict_proba(
        corrected_live_X[FEATURES]
    )[:, 1]

    print_report(
        "After correcting the upstream unit",
        build_drift_report(
            reference_X,
            corrected_live_X,
            FEATURES,
            reference_score,
            corrected_score,
        ),
    )

    recent_X, recent_y = make_orders(
        40_000,
        seed=3,
        concept_shift=True,
        joint_shift=True,
    )
    shadow_X, shadow_y = make_orders(
        15_000,
        seed=4,
        concept_shift=True,
        joint_shift=True,
    )

    champion_shadow_score = (
        champion.predict_proba(
            shadow_X[FEATURES]
        )[:, 1]
    )

    print_report(
        "Stable marginals, shifted predictions",
        build_drift_report(
            reference_X,
            shadow_X,
            FEATURES,
            reference_score,
            champion_shadow_score,
        ),
    )

    challenger = train_model(
        recent_X,
        recent_y,
        seed=43,
    )
    challenger_path = (
        ARTIFACTS / "challenger.json"
    )
    challenger.save_model(challenger_path)

    challenger_shadow_score = (
        challenger.predict_proba(
            shadow_X[FEATURES]
        )[:, 1]
    )

    champion_metrics = performance_metrics(
        shadow_y,
        champion_shadow_score,
    )
    challenger_metrics = performance_metrics(
        shadow_y,
        challenger_shadow_score,
    )

    delta, ci_low, ci_high = (
        paired_bootstrap_log_loss_delta(
            shadow_y,
            champion_shadow_score,
            challenger_shadow_score,
        )
    )

    comparison = {
        "champion": champion_metrics,
        "challenger": challenger_metrics,
        "challenger_minus_champion_log_loss": (
            delta
        ),
        "log_loss_delta_95_percent_interval": [
            ci_low,
            ci_high,
        ],
    }

    print(
        "\nSame-traffic champion/challenger comparison"
    )
    print(
        json.dumps(
            comparison,
            indent=2,
        )
    )

    promote = (
        challenger_metrics["log_loss"]
        <= champion_metrics["log_loss"] - 0.01
        and ci_high < 0.0
        and challenger_metrics["roc_auc"]
        >= champion_metrics["roc_auc"] - 0.002
        and challenger_metrics[
            "top_decile_capture"
        ]
        >= champion_metrics[
            "top_decile_capture"
        ]
    )

    shadow_log = shadow_X.copy()
    shadow_log.insert(
        0,
        "request_id",
        np.arange(len(shadow_log)),
    )
    shadow_log["champion_score"] = (
        champion_shadow_score
    )
    shadow_log["challenger_score"] = (
        challenger_shadow_score
    )
    shadow_log["label"] = shadow_y
    shadow_log.to_csv(
        ARTIFACTS / "shadow_scores.csv",
        index=False,
    )

    if promote:
        registry = {
            "champion": {
                "version": (
                    "returns-model-2026-08-03"
                ),
                "model_path": str(
                    challenger_path
                ),
                "promoted_at": (
                    datetime.now(
                        timezone.utc
                    ).isoformat()
                ),
            },
            "previous_champion": {
                "version": (
                    "returns-model-2026-05-01"
                ),
                "model_path": str(
                    champion_path
                ),
            },
            "evidence": comparison,
        }

        atomic_write_json(
            ARTIFACTS / "registry.json",
            registry,
        )
        print(
            "\nPromotion decision: "
            "promote challenger."
        )
    else:
        print(
            "\nPromotion decision: "
            "keep champion."
        )

    history = pd.DataFrame(
        {
            "window_end": pd.date_range(
                "2026-05-17",
                periods=8,
                freq="W-SUN",
            ),
            "roc_auc": [
                0.809,
                0.807,
                0.804,
                0.800,
                0.796,
                0.792,
                0.786,
                0.780,
            ],
        }
    )

    cadence = forecast_retraining_trigger(
        history,
        metric_floor=0.770,
        preparation_days=14,
        safety_buffer_days=7,
    )

    print(
        "\nDecay-derived retraining trigger"
    )
    print(
        json.dumps(
            cadence,
            indent=2,
        )
    )


if __name__ == "__main__":
    main()

Run it:

python run_demo.py

The exact floating-point values can vary slightly by platform, but the important shape of the result should be similar to:

Silent unit-change incident

seconds_since_last_order
PSI: approximately 12
KS: approximately 1.0
Status: critical

predicted_return_probability
PSI: greater than 1
Status: critical

After correcting the upstream unit
All major feature and prediction alerts return to normal.

Stable marginals, shifted predictions
Individual feature PSI values remain close to zero.
Prediction PSI crosses the tighter warning threshold.

Same-traffic champion/challenger comparison
The challenger improves log loss, ROC AUC, and top-decile capture.
The paired bootstrap interval for the log-loss difference remains below zero.

Why prediction drift can fire when features look stable

The demo deliberately changes the association between order value and discount while preserving the one-dimensional discount distribution.

Imagine these four feature pairs:

Reference:
high order value + high discount
low order value  + low discount

Current:
high order value + low discount
low order value  + high discount

The list of order values is unchanged. The list of discounts is unchanged. Therefore:

Order-value PSI can remain near zero.
Discount PSI can remain near zero.
Their individual KS statistics can remain near zero.

But the joint distribution changed.

A tree ensemble can split first on order value, then on discount. The frequency with which rows reach particular leaf combinations changes, causing the score distribution to move.

Prediction monitoring is therefore a low-cost detector for some multivariate changes. It does not tell you which relationship changed, but it tells you that the model’s view of live traffic is no longer the same.

When prediction drift fires without obvious feature drift, investigate:

Pairwise correlations.
Segment proportions.
Interaction features.
Category combinations.
Missingness patterns conditional on other features.
Feature freshness.
Whether a request router changed which traffic reaches the model.
Whether preprocessing and model versions are correctly paired.

Do not conclude that the model needs retraining merely because prediction PSI crossed 0.05.

What to do when drift fires

The safest incident-response order is:

1. Preserve evidence

Record:

Alerting window.
Model version.
Feature-schema version.
Pipeline deployment version.
Raw-source version.
Reference version.
A sample of request identifiers.
Segment-level drift metrics.
Current decision rates.

Avoid replacing the reference or mutating logs while the incident is being investigated.

2. Check data contracts

Verify:

Column names and data types.
Units.
Currency.
Time zones.
Allowed ranges.
Category vocabularies.
Nullability.
Default values.
Array or vector dimensions.
Timestamp freshness.
Duplicate rates.

A feature called duration is insufficiently specified. Prefer a contract such as duration_seconds, accompanied by a documented valid range and transformation owner.

3. Inspect pipeline changes

Ask whether an upstream service recently changed:

Serialization formats.
Database fields.
Aggregation windows.
Join keys.
Deduplication rules.
Time-zone conversions.
Currency conversions.
Missing-value handling.
Batch schedules.
Source tables.

Compare raw source values, transformed values, and model-ready values for the same request.

4. Check training-serving consistency

Recompute a sample of production features through the offline training pipeline and compare it with online values.

This can reveal:

Different aggregation logic.
Different category encodings.
Stale online features.
Point-in-time leakage in training data.
Different default values.
A bug deployed only to one path.

5. Check expected events and seasonality

A promotion, holiday, product launch, policy change, or market expansion can create legitimate drift.

Legitimate does not mean harmless. The model may still require intervention, but the remediation differs from fixing corrupted data.

6. Evaluate labeled performance

When labels are available, compare performance by event cohort—not by label-arrival date.

An order created in May but labeled in June belongs to the May prediction cohort. Otherwise, changing label delays can make performance charts misleading.

7. Retrain only after the input is trusted

Retraining on corrupted milliseconds-as-seconds data teaches the model to accommodate a broken contract.

That can turn a temporary pipeline bug into a permanent model dependency.

Cherry on the cake: the model was fine, but the unit was wrong

A Zalando Payments presentation on live monitoring described a particularly instructive failure. A fraud model had been trained with a time-to-order feature measured in seconds. In live serving, another microservice sent the feature in milliseconds.

The feature mean moved from approximately 200 to 200,000. The model service itself remained operational, and conventional CPU, memory, and latency monitoring did not identify the semantic corruption. The presentation reported that the resulting predictions were corrupt and used the incident to motivate live comparison of feature and prediction distributions.

This is why “investigate upstream first” is not a slogan. It is a safeguard against solving the wrong problem.

A range contract could have caught the issue before inference:

def validate_seconds_since_last_order(
    value: float,
) -> None:
    if not np.isfinite(value):
        raise ValueError(
            "seconds_since_last_order must be finite"
        )

    if value < 0 or value > 31_536_000:
        raise ValueError(
            "seconds_since_last_order is outside "
            "the supported range of 0 to one year"
        )

A stronger contract also carries an explicit unit:

{
  "name": "seconds_since_last_order",
  "data_type": "float64",
  "unit": "seconds",
  "minimum": 0,
  "maximum": 31536000,
  "nullable": false,
  "owner": "orders-feature-platform",
  "schema_version": "4.2.0"
}

Champion/challenger means the same traffic, not merely the same test file

An offline comparison is necessary, but production promotion should include shadow scoring.

In a shadow deployment:

The champion continues to drive decisions.
The challenger receives the same model-ready feature vector.
Both predictions are logged under the same request identifier.
Only the champion’s result affects the customer or downstream system.
Delayed labels are joined to both scores later.

This design controls for traffic differences. If one model receives Monday traffic and the other receives Tuesday traffic, model quality is confounded with population differences.

Each shadow log should include:

Request or entity identifier.
Event timestamp.
Champion version.
Challenger version.
Feature-schema hash.
Champion score.
Challenger score.
Champion decision.
Challenger counterfactual decision.
Decision threshold.
Segment keys.
Label and label timestamp when available.
End-to-end inference latency for both models.

Promotion criteria should be written before the experiment

Do not inspect ten metrics and promote based on whichever one looks favorable.

A defensible promotion policy might require:

Challenger log loss improves by at least 0.01.
A paired 95% bootstrap interval for challenger-minus-champion log loss remains below zero.
ROC AUC is not worse by more than 0.002.
Capture rate at the available review capacity does not decline.
No important segment violates its guardrail.
Calibration remains acceptable.
Latency and memory stay within service limits.
The challenger passes schema and data-quality checks.
The rollback artifact has been tested.

The runnable example uses paired bootstrap resampling because both models score the same rows. Pairing preserves the row-level relationship and produces a more informative comparison than independently resampling each model’s results.

Promotion should be atomic and reversible

The example writes a temporary registry file and replaces the active registry atomically.

A production registry should also retain:

Previous champion.
Training-data cutoff.
Training-code commit.
Feature-pipeline version.
Dependency lock file.
Approval record.
Evaluation report.
Promotion timestamp.
Rollback instructions.

Do not overwrite the previous model artifact. A rollback should be a pointer change, not an emergency retraining job.

Derive retraining cadence from measured decay

“Retrain every month” is easy to schedule and difficult to justify.

A better cadence begins with a business performance floor.

Suppose the model’s ROC AUC is currently 0.81, and the minimum acceptable level is 0.77. Labeled weekly cohorts show a downward trend.

Estimate the time at which the metric will cross the floor, then subtract:

Time required to collect enough recent labels.
Training and tuning time.
Validation time.
Shadow-evaluation time.
Approval and deployment time.
A safety buffer.

In plain text:

Recommended retraining start = forecast floor-breach date − preparation lead time − safety buffer

If the model is forecast to breach its floor in 35 days, while training, validation, and shadowing require 21 days, waiting for a monthly calendar job may already be too late.

A practical decay-derived process

Choose the performance floor from a business or risk requirement.
Measure performance on fixed event-time cohorts.
Attach confidence intervals and sample counts.
Exclude confirmed data incidents from the decay fit.
Fit a trend over enough recent labeled windows.
Forecast the breach date.
Subtract operational lead time.
Recalculate whenever a new labeled cohort matures.

The included function uses a transparent linear trend. Production systems may need:

Robust regression.
Exponential decay.
Piecewise trends.
Seasonal components.
Bayesian change-point models.
Separate trends by important segment.

The objective is not to predict the crossing date to the hour. It is to replace an unsupported calendar interval with an evidence-based operating horizon.

Keep a maximum-age policy as a backstop

A measured-decay policy and a calendar policy can coexist.

For example:

Start early when forecast decay demands it.
Require a review after a maximum model age even when metrics appear stable.
Retrain immediately only when concept decay is confirmed and the release can be validated safely.

The maximum-age rule protects against weak labels, blind spots, and low statistical power. It should not be the primary reason for every retraining run.

Alert design that works in production

A single noisy window should rarely wake someone at 03:00.

A mature alert policy includes:

A minimum row count.
A minimum number of non-missing observations.
Persistence, such as two of three consecutive windows.
Severity based on both drift magnitude and affected feature importance.
Separate thresholds for features and predictions.
Seasonally matched references.
Segment-level monitoring for high-risk populations.
Suppression during known migrations, with an expiry time.
Links to dashboards, samples, owners, and runbooks.

Consider two channels:

A ticket or chat warning for moderate drift.
A page only for severe drift combined with operational impact, data-contract failure, or decision-rate movement.

A critical unit violation should page quickly. A small PSI increase in a low-impact feature can wait for business hours.

Do not let delayed labels create a blind period

When labels take 30 days, there are effectively two monitoring speeds.

The fast loop uses:

Schema checks.
Unit and range checks.
Missingness.
Feature PSI.
Feature KS.
Prediction PSI.
Prediction KS.
Decision rates.
Traffic and segment mix.

The slow loop uses:

ROC AUC.
Log loss.
Calibration.
Precision and recall.
Business outcomes.
Segment fairness and reliability metrics.
Champion/challenger comparisons.

The fast loop reduces time to detection. The slow loop determines whether the model remains useful.

Neither replaces the other.

Production hardening checklist

Before calling the system complete, verify that it can answer these questions:

Which exact reference window produced this alert?
Were the same transformations used in reference and serving?
Did the feature’s unit or semantic definition change?
Is the shift global or isolated to a segment?
Did the predicted-score distribution move?
Did the decision rate move?
Are labels complete for the cohort being evaluated?
Is performance loss statistically and operationally meaningful?
Did the challenger score the same requests as the champion?
Can the system roll back without rebuilding anything?
Is the next retraining date based on measured decay and actual lead time?
Who owns the first upstream investigation?

Start by running the example against your own model’s reference and recent production windows. Replace the generic thresholds with values backtested on your historical traffic, add unit-bearing feature contracts, and put the next candidate model into shadow mode before it is allowed to make a single live decision.