StableGLM

rashomon-py

Does your conclusion survive every equally-good model?

Many parameter vectors fit a training set about as well as the one your solver returned. If some of them reverse your conclusion (a coefficient’s sign, a feature ranking, a decision for one row), the conclusion depends on which optimum the solver landed on, not on the data. rashomon-py checks this for a fitted scikit-learn logistic or linear regression.

from sklearn.linear_model import LogisticRegression
from rashomon import audit

model = LogisticRegression().fit(X, y)       # your existing model (Pipelines work too)
report = audit(model, X, y)                  # pandas in -> feature names out

print(report.summary())
report.plot()
report.flipped                               # bool mask: rows whose prediction can flip
report.coefficients                          # DataFrame: estimate, low, high, sign_stable

Output for the ten “mean” features of the Wisconsin breast-cancer data with the default LogisticRegression():

Stability audit: LogisticRegression (classification), n=569, 10 features
Equally-good models: training log-loss within 0.01704 of the optimum 0.1435
  (one standard error of the 5-fold cross-validated log-loss, whose mean is 0.1458; models closer than this cannot be told apart by CV)
Method: hit-and-run sampling of the exact set (with ellipsoid proposals); 2,000 models, min ESS = 218 -> reliable

Predictions (flip test exact; disagreement over the sampled models)
  Flip under some equally-good model:           17.4%  (99 of 569)
  Worst single-model disagreement with yours:    5.1%

Coefficients (exact range across all equally-good models)
  feature                  estimate        low       high  sign
  mean radius               -0.9976     -4.634      2.621  UNSTABLE
  mean texture               -1.399     -2.527    -0.4664  stable
  mean perimeter            -0.9124     -4.613      2.775  UNSTABLE
  mean area                  -1.298     -5.061      2.433  UNSTABLE
  ...
Sign stable: 1 of 10 features.

audit plot

The model is 95% accurate. Still, one diagnosis in six is reversed by some model that cross-validation cannot distinguish from it, and only the coefficient on texture keeps its sign across all of them. Radius, perimeter and area are near-duplicates, so the data fixes their combined effect but not how to split it between them. A claim like “tumour radius lowers the odds” would not survive.

Install

pip install rashomon-py          # Python 3.10+; depends on numpy, scipy, scikit-learn, pandas, matplotlib

How this differs from a confidence interval

A confidence interval or bootstrap asks how much the estimate would move under a new sample from the population. rashomon-py asks how many different models fit the data you have about equally well, and whether they agree with yours. The first is sampling uncertainty; the second is model multiplicity (Breiman’s “Rashomon effect”). A coefficient can have a narrow confidence interval and still change sign across equally-good models when features are collinear, because the loss is nearly flat along the collinear directions. See Bootstrap, Rashomon, and Bayesian intervals for a worked comparison.

What “equally good” means

The tolerance sets how much worse than optimal a model may be and still count. audit() accepts four forms and prints the one in use:

tolerance= Meaning When to use
"cv" (default) one standard error of the cross-validated loss; models closer than this cannot be told apart by cross-validation (the one-standard-error rule from glmnet) data-driven default
0.01 (any float in (0,1)) models at most 1% worse than optimal on the training loss a rule that is easy to state; 0.01 to 0.05 are common
"lr" / ("lr", 0.05) the models not rejected by a likelihood-ratio test at level α unpenalized fits
("profile", 0.05) χ²₁(0.95)/(2n): for an unpenalized fit the coefficient ranges are the 95% profile-likelihood confidence intervals statistical reporting; see Validation below
("absolute", 0.002) a loss gap in training-loss units reproducing a published setting

The CV default is permissive, so expect larger sets than with a 1% rule. Results at two or three tolerances say more than any single one; see Choosing the tolerance.

What the report contains

report. Meaning Term in the literature
flip_rate, flipped share (and mask) of rows whose predicted label changes under some equally-good model ambiguity (Marx, Calmon & Ustun 2020)
max_disagreement the largest share of rows on which one equally-good model disagrees with yours discrepancy (Marx et al. 2020)
coefficients each coefficient’s range across all equally-good models, and whether its sign is stable variable importance cloud (Dong & Rudin 2020); hacking intervals (Coker, Rudin & King 2021)
prediction_ranges per-row range of predicted probability or fitted value prediction bands over the Rashomon set
predict_ranges(X_new) the same for new rows, through your pipeline  
tolerance the loss gap defining the set ε in the ε-Rashomon set (Fisher, Rudin & Dominici 2019)
rashomon_set, samples the underlying RashomonSet and the sampled parameter vectors  

Coefficient ranges and the flip test are exact: each coefficient range is the solution of a convex program over the true set, and each row not already flipped by a sampled model is settled by one more convex program (automatic when the problem is small enough: n · d² ≤ 5·10⁶ for ranges, n_undecided · n · d² ≤ 2·10¹⁰ for flips; force with exact_ranges=True / exact_flips=True). The disagreement figure and the prediction ranges come from models sampled from the exact set, so they are lower bounds that tighten as n_samples grows. The report states the effective sample size and whether it is reliable.

Validation

For an unpenalized logistic regression, the range of a coefficient over the Rashomon set with tolerance χ²₁(0.95)/(2n) is by definition its 95% profile-likelihood confidence interval, so the exact-range machinery can be checked against a standard statistical result. On the UCLA graduate-admissions data (admit ~ gre + gpa + rank, n = 400), coef_extremes() reproduces all twelve bounds of R’s confint() output to within 1e-4 (R’s own interpolation precision) and the MLE to 5e-7; the test also agrees with an independent profile root-finder to 2e-6. A second test grids a two-dimensional Rashomon set and checks the exact ranges, the exact flip test and the sampler’s moments against brute force. See tests/test_validation.py. The evaluation compares exact, sampled and ellipsoid quantities on four real datasets.

Supported models

LogisticRegression / LogisticRegressionCV (binary; L2 or unpenalized; with or without class_weight), Ridge / RidgeCV, LinearRegression, and a Pipeline whose last step is one of these. If the model was fitted with sample_weight, pass the same weights to audit(..., sample_weight=w). The regularization strength, the unpenalized intercept and the row weights are converted exactly, and the audit checks that your fitted coefficients sit at the optimum of the reconstructed objective. A mismatch is reported.

Not supported: multiclass (planned), L1 and elastic-net penalties (the level set of a non-smooth objective is not a convex set of the same kind), trees, neural networks. The method needs a convex, twice-differentiable loss, so other L2-penalized GLMs (Poisson, multinomial) are possible extensions.

Coefficient ranges and flips are exact, so dimension affects their cost, not their validity (about a minute at n = 5,000 and d = 101). The sampled quantities depend on the hit-and-run chain mixing, which it does well up to a few dozen features and poorly beyond about 60; the report’s effective sample size says which case applies. Sampled extremes understate: with 2,000 draws they cover about half the exact coefficient range at d = 30 and a fifth at d = 100, which is why the exact computations are the default. See the evaluation.

Expert API

RashomonSet gives direct access to ε calibration (percent_loss, LR_alpha, absolute), the membership oracle, hit-and-run sampling (with optional ellipsoid proposals), the Hessian-ellipsoid approximation, functional_range / coef_extremes, model class reliance, Shapley-VIC, and bootstrap and Bayesian comparisons. RashomonSet.from_sklearn(model, X, y, ...) builds one from a fitted model. RashomonSet(C=...) is not scikit-learn’s C: it is 1/λ for the mean-loss objective, and scikit-learn’s C equals C/n. audit and from_sklearn do the conversion.

from rashomon import RashomonSet
rs = RashomonSet.from_sklearn(model, X, y, epsilon=0.02, random_state=0)
rs.functional_range(x_row)          # exact prediction range for one row, on the logit scale
rs.sample_hitandrun(2000, ellipsoid_mix=0.5)

Documentation

References

License

MIT

This site is open source. Improve this page.