Model Stability

Model stability is how much a model’s behavior changes when the thing it was fitted to changes a little: a different sample, a different split, a different week. The usual treatment lists four dimensions — performance, prediction, feature importance, temporal — and then proposes measuring one of them. They do not move together.

The stable model that keeps changing its mind

Cross-validation variance is the standard answer, and it measures the stability of a performance estimate. That is not the same as the stability of an answer given to a person.

The experiments in this entry run in one session. The first two share one routine, bootstrap_models: fit k logistic regressions, each on its own bootstrap resample of the same training set. Every difference between the runs then comes from resampling alone. Here the question is what the thirty retrained models look like on one fixed test set, first through the AUC of each run and then through the predicted probability each customer receives from each run.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score

rng = np.random.default_rng(4242)
d, n = 18, 3000
w = rng.normal(size=d) * 0.45


def make(m):
    X = rng.normal(size=(m, d))
    return X, (rng.random(m) < 1 / (1 + np.exp(-X @ w))).astype(int)


def bootstrap_models(X, y, k, C=1.0):
    """Fit k logistic regressions, each on its own bootstrap resample of (X, y)."""
    models = []
    for _ in range(k):
        idx = np.floor(rng.random(len(X)) * len(X)).astype(int)
        models.append(LogisticRegression(C=C, max_iter=5000).fit(X[idx], y[idx]))
    return models


Xtr, ytr = make(n)
Xte, yte = make(12000)
runs = np.array([m.predict_proba(Xte)[:, 1] for m in bootstrap_models(Xtr, ytr, 30)])
aucs = np.array([roc_auc_score(yte, p) for p in runs])

# Reading 1: the performance estimate across the 30 retrainings.
print("AUC   min %.4f   max %.4f   range %.4f   sd %.4f"
      % (aucs.min(), aucs.max(), aucs.max() - aucs.min(), aucs.std(ddof=1)))

# Reading 2: what each customer experiences on the same 30 runs.
per = runs.std(axis=0, ddof=1)
print("spread of a customer's predicted probability, median %.4f" % np.median(per))
print("90th percentile of that spread                       %.4f" % np.quantile(per, 0.9))
for cut in (0.3, 0.5, 0.7):
    above = (runs > cut).mean(axis=0)
    crossed = (above > 0) & (above < 1)
    print("share of customers who cross the %.1f line at least once  %.2f%%"
          % (cut, crossed.mean() * 100))

AUC   min 0.8597   max 0.8639   range 0.0043   sd 0.0012
spread of a customer's predicted probability, median 0.0317
90th percentile of that spread                       0.0530
share of customers who cross the 0.3 line at least once  15.71%
share of customers who cross the 0.5 line at least once  16.83%
share of customers who cross the 0.7 line at least once  14.45%

In this simulation, thirty bootstrap retrainings land inside 0.0043 of AUC, so a check that looks only at the aggregate spread of AUC would report a stable model. On the same thirty runs, 16.83% of customers cross the 0.5 line at least once: for one person in six, whether they were flagged depended on which retraining happened to be live that day, and they did nothing differently.

Both readings are stability. The first is what a modeling team reports and the second is what a customer, an auditor or an appeals process experiences. They can be measured on exactly the same runs, so there is no reason to report only the first.

Correlated-feature importance can vary more than predictions

Comparing feature rankings across retrainings is often proposed as a stability check. It measures something real, and the thing it measures is not the model’s behavior. Here the same bootstrap_models routine is applied twice, with two of the four features uncorrelated and then correlated at 0.97, and the code records each model’s coefficients next to its predictions.

rng = np.random.default_rng(7788)
n = 4000


def make(m, rho):
    z = rng.normal(size=(m, 4))
    a = z[:, 0]
    b = rho * a + np.sqrt(1 - rho ** 2) * z[:, 1]
    X = np.column_stack([a, b, z[:, 2], z[:, 3]])
    logit = 0.7 * a + 0.7 * b + 0.5 * z[:, 2]
    return X, (rng.random(m) < 1 / (1 + np.exp(-logit))).astype(int)


for rho in (0.0, 0.97):
    Xtr, ytr = make(n, rho)
    Xte, yte = make(8000, rho)
    models = bootstrap_models(Xtr, ytr, 30)
    coefs = np.array([m.coef_[0] for m in models])
    preds = np.array([m.predict_proba(Xte)[:, 1] for m in models])
    aucs = [roc_auc_score(yte, p) for p in preds]
    top = np.argmax(np.abs(coefs), axis=1)
    print("features 1 and 2 correlated at %.2f" % rho)
    print("   AUC sd across retrainings                 %.4f" % np.std(aucs, ddof=1))
    print("   median spread of a customer's prediction  %.4f"
          % np.median(preds.std(axis=0, ddof=1)))
    print("   coefficient sd, feature 1                 %.4f" % coefs[:, 0].std(ddof=1))
    print("   coefficient sd, feature 2                 %.4f" % coefs[:, 1].std(ddof=1))
    print("   how often feature 1 ranks most important  %.0f%%" % ((top == 0).mean() * 100))

features 1 and 2 correlated at 0.00
   AUC sd across retrainings                 0.0004
   median spread of a customer's prediction  0.0137
   coefficient sd, feature 1                 0.0464
   coefficient sd, feature 2                 0.0448
   how often feature 1 ranks most important  67%
features 1 and 2 correlated at 0.97
   AUC sd across retrainings                 0.0003
   median spread of a customer's prediction  0.0138
   coefficient sd, feature 1                 0.1359
   coefficient sd, feature 2                 0.1502
   how often feature 1 ranks most important  43%

Correlating two features leaves the predictions and the AUC where they were and roughly triples the coefficient spread. Which of the two the importance table calls strongest becomes a coin flip, 43% against 57%, where before it was 67% against 33%.

Two correlated features can trade the same signal back and forth without changing a single output. Importance stability is a property of the parameterization, and treating it as a proxy for model stability leads to retraining a model that is behaving perfectly — or to a governance conversation about a feature that never mattered.

What actually buys stability

The standard remedies are regularization and ensembling. Both are worth measuring rather than assuming, because the variation being complained about has a specific source. This third experiment starts from a new data generator with 40 features and builds the bootstrap refit into evaluate instead of calling bootstrap_models, because it needs every customer’s rank and flag, not coefficients. Each variant is refit on 25 bootstrap resamples and scored on one fixed test set with the rule “flag the top 20% of scores”, so the table reports mean AUC, the spread of each customer’s rank, and the share of customers who move in or out of the flagged group.

rng = np.random.default_rng(31337)
d = 40
w = rng.normal(size=d) * 0.4


def make(m):
    X = rng.normal(size=(m, d))
    return X, (rng.random(m) < 1 / (1 + np.exp(-X @ w))).astype(int)


Xte, yte = make(10000)


def fit_score(X, y, C):
    return LogisticRegression(C=C, max_iter=5000).fit(X, y).predict_proba(Xte)[:, 1]


def bagged(X, y, k):
    out = np.zeros(len(Xte))
    for _ in range(k):
        j = np.floor(rng.random(len(X)) * len(X)).astype(int)
        out += fit_score(X[j], y[j], 100.0)
    return out / k


def evaluate(pool, builder, reps=25):
    """Refit on bootstrap resamples of the pool; return mean AUC, rank spread, flip share."""
    Xp, yp = pool
    ranks, aucs, flags = [], [], []
    for _ in range(reps):
        idx = np.floor(rng.random(len(Xp)) * len(Xp)).astype(int)
        p = builder(Xp[idx], yp[idx])
        aucs.append(roc_auc_score(yte, p))
        r = np.argsort(np.argsort(p)) / len(p)
        ranks.append(r)
        flags.append(r >= 0.80)
    ranks, flags = np.array(ranks), np.array(flags)
    flips = (flags.mean(0) > 0) & (flags.mean(0) < 1)
    return np.mean(aucs), np.median(ranks.std(axis=0, ddof=1)), flips.mean() * 100


plain = lambda X, y: fit_score(X, y, 100.0)
small = make(700)
results = {
    "700 rows, plain fit": evaluate(small, plain),
    "700 rows, C = 0.01": evaluate(small, lambda X, y: fit_score(X, y, 0.01)),
    "700 rows, averaged over 20 fits": evaluate(small, lambda X, y: bagged(X, y, 20)),
    "2800 rows, plain fit": evaluate(make(2800), plain),
    "11200 rows, plain fit": evaluate(make(11200), plain),
}

print("ceiling AUC %.4f" % roc_auc_score(yte, Xte @ w))
print("%-34s %9s %12s %16s" % ("", "mean AUC", "rank spread", "flips in or out"))
for label, (auc, spread, flip) in results.items():
    print("%-34s %9.4f %12.4f %15.1f%%" % (label, auc, spread, flip))

ceiling AUC 0.8810
                                    mean AUC  rank spread  flips in or out
700 rows, plain fit                   0.8435       0.0929            34.5%
700 rows, C = 0.01                    0.8451       0.0846            31.9%
700 rows, averaged over 20 fits       0.8454       0.0914            33.9%
2800 rows, plain fit                  0.8736       0.0400            14.9%
11200 rows, plain fit                 0.8791       0.0209             7.8%

Regularization and bootstrap averaging make relatively small changes for the stable logistic estimator in this simulation; sixteen times the data reduces spread by 77%. This is not an impossibility result for bagging. Resampling adds no new observations, but averaging unstable fitted estimators can reduce variation across training samples. Measure the effect for the actual learner and prediction metric.

Instability can reflect limited information, an unstable learner, weak identification or changing data. Additional representative labels, regularization, model changes and ensembling are candidate remedies; select them using repeated-fit evaluation rather than assuming one lever is always effective.

Performance degradation over time can result from concept drift, population mix, pipeline faults, label changes or measurement problems. Diagnose those causes separately from sensitivity to training samples; drift guardrails provide evidence rather than a unique diagnosis.

References

• scikit-learn documentation, sklearn.linear_model.LogisticRegression.
• Gama, J., Žliobaitė, I., Bifet, A., Pechenizkiy, M., & Bouchachia, A. (2014). A survey on concept drift adaptation. ACM Computing Surveys, 46(4), 1–37.