API reference

Estimators

FlagGAMClassifier / FlagGAMRegressor — the sklearn-compatible entry points.

estimator

FlagGAMClassifier and FlagGAMRegressor: sklearn-compatible estimators for rule-basis GAMs.

FlagGAMClassifier

FlagGAMClassifier(task='auto', quantile_low=(0.05, 0.45), quantile_high=(0.55, 0.95), quantile_step=0.05, min_support='auto', fdr_alpha=0.05, effect_size='risk_difference', representation='full', feature_weighting=None, head='additive', flexible_estimator=None, C=1.0, missing='no_evidence', monotonic_constraints=None, categorical_features=None, random_state=None)

Bases: ClassifierMixin, _BaseFlagGAM

Rule-basis generalized additive model classifier (sklearn contract).

Source code in src/flaggam/estimator.py

def __init__(
    self,
    task="auto",
    quantile_low=(0.05, 0.45),
    quantile_high=(0.55, 0.95),
    quantile_step=0.05,
    min_support="auto",
    fdr_alpha=0.05,
    effect_size="risk_difference",
    representation="full",
    feature_weighting=None,
    head="additive",
    flexible_estimator=None,
    C=1.0,
    missing="no_evidence",
    monotonic_constraints=None,
    categorical_features=None,
    random_state=None,
) -> None:
    self.task = task
    self.quantile_low = quantile_low
    self.quantile_high = quantile_high
    self.quantile_step = quantile_step
    self.min_support = min_support
    self.fdr_alpha = fdr_alpha
    self.effect_size = effect_size
    self.representation = representation
    self.feature_weighting = feature_weighting
    self.head = head
    self.flexible_estimator = flexible_estimator
    self.C = C
    self.missing = missing
    self.monotonic_constraints = monotonic_constraints
    self.categorical_features = categorical_features
    self.random_state = random_state

fit

fit(X, y)

Discover flag bases, fit head on Z(X).

Source code in src/flaggam/estimator.py

def fit(self, X, y):
    """Discover flag bases, fit head on Z(X)."""
    if y is None:
        raise ValueError("requires y to be passed, but the target y is None")
    if self.task not in {"auto", "binary", "multiclass"}:
        raise ValueError(
            f"task={self.task!r} is not recognised; "
            "valid values are 'auto', 'binary', 'multiclass'"
        )
    if self.representation not in {"full", "compact"}:
        raise ValueError(
            f"representation={self.representation!r} is not recognised; "
            "valid values are 'full', 'compact'"
        )
    if self.missing not in {"no_evidence", "indicator"}:
        raise ValueError(
            f"missing={self.missing!r} is not recognised; "
            "valid values are 'no_evidence', 'indicator'"
        )
    df = self._to_frame(X, reset=True)
    y = np.asarray(y)
    if len(df) != len(y):
        raise ValueError(
            f"Found input variables with inconsistent numbers of samples: "
            f"[{len(df)}, {len(y)}]"
        )
    check_classification_targets(y)
    self.label_encoder_ = LabelEncoder().fit(y)
    self.classes_ = self.label_encoder_.classes_
    y_enc = self.label_encoder_.transform(y)

    task = self.task
    if task == "auto":
        task = "binary" if len(self.classes_) == 2 else "multiclass"

    if self.monotonic_constraints is not None:
        _validate_constraints(self.monotonic_constraints, list(df.columns))
        if self.head != "additive":
            raise ValueError("monotonic_constraints require head='additive'")
        if self.representation == "compact":
            raise ValueError(
                "monotonic_constraints are incompatible with representation='compact'"
            )
        if task == "multiclass":
            raise NotImplementedError(
                "monotonic constraints support binary classification only"
            )
    if self.head == "flexible" and self.flexible_estimator is None:
        raise ValueError("head='flexible' requires flexible_estimator")

    self.core_ = self._build_core(task).fit(df, y_enc)
    Z = self.core_.transform(df)

    self.feature_weights_ = None
    H = Z
    if self.representation == "compact":
        if self.feature_weighting == "auto":
            self.feature_weights_ = feature_weights(
                df, y_enc, task, self.core_.numerical_features_
            )
        H = compact_scores(
            Z,
            self.core_.bases_,
            np.arange(len(self.classes_)),
            self.feature_weights_,
        )

    if self.head == "additive":
        if self.monotonic_constraints is not None:
            from .monotonic import MonotonicAdditiveHead, bounds_for_bases

            C = 1.0 if isinstance(self.C, (list, tuple)) else self.C
            bounds = bounds_for_bases(self.core_.bases_, dict(self.monotonic_constraints))
            self.head_ = MonotonicAdditiveHead(task, bounds, C=C)
        else:
            self.head_ = AdditiveHead(task, C=self.C, random_state=self.random_state)
    else:
        self.head_ = FlexibleHead(self.flexible_estimator, task, random_state=self.random_state)
    self.head_.fit(H, y_enc)
    return self

FlagGAMRegressor

FlagGAMRegressor(quantile_low=(0.05, 0.45), quantile_high=(0.55, 0.95), quantile_step=0.05, min_support='auto', fdr_alpha=0.05, head='additive', flexible_estimator=None, alpha=1.0, missing='no_evidence', monotonic_constraints=None, categorical_features=None, random_state=None)

Bases: RegressorMixin, _BaseFlagGAM

Rule-basis generalized additive model regressor (sklearn contract).

Source code in src/flaggam/estimator.py

def __init__(
    self,
    quantile_low=(0.05, 0.45),
    quantile_high=(0.55, 0.95),
    quantile_step=0.05,
    min_support="auto",
    fdr_alpha=0.05,
    head="additive",
    flexible_estimator=None,
    alpha=1.0,
    missing="no_evidence",
    monotonic_constraints=None,
    categorical_features=None,
    random_state=None,
) -> None:
    self.quantile_low = quantile_low
    self.quantile_high = quantile_high
    self.quantile_step = quantile_step
    self.min_support = min_support
    self.fdr_alpha = fdr_alpha
    self.head = head
    self.flexible_estimator = flexible_estimator
    self.alpha = alpha
    self.missing = missing
    self.monotonic_constraints = monotonic_constraints
    self.categorical_features = categorical_features
    self.random_state = random_state

fit

fit(X, y)

Discover flag bases, fit regression head on Z(X).

Source code in src/flaggam/estimator.py

def fit(self, X, y):
    """Discover flag bases, fit regression head on Z(X)."""
    if y is None:
        raise ValueError("requires y to be passed, but the target y is None")
    if self.missing not in {"no_evidence", "indicator"}:
        raise ValueError(
            f"missing={self.missing!r} is not recognised; "
            "valid values are 'no_evidence', 'indicator'"
        )
    df = self._to_frame(X, reset=True)
    y = np.asarray(y)
    if y.dtype.kind == "c":
        raise ValueError(f"Complex data not supported\n{y}\n")
    if y.ndim == 2 and y.shape[1] == 1:
        warnings.warn(
            "A column-vector y was passed when a 1d array was expected. "
            "Please change the shape of y to (n_samples,), for example using ravel().",
            DataConversionWarning,
            stacklevel=2,
        )
    y = y.astype(float).ravel()
    if len(df) != len(y):
        raise ValueError(
            f"Found input variables with inconsistent numbers of samples: "
            f"[{len(df)}, {len(y)}]"
        )
    if self.monotonic_constraints is not None:
        _validate_constraints(self.monotonic_constraints, list(df.columns))
        if self.head != "additive":
            raise ValueError("monotonic_constraints require head='additive'")
    if self.head == "flexible" and self.flexible_estimator is None:
        raise ValueError("head='flexible' requires flexible_estimator")
    self.core_ = self._build_core("regression").fit(df, y)
    Z = self.core_.transform(df)
    if self.head == "additive":
        if self.monotonic_constraints is not None:
            from .monotonic import MonotonicAdditiveHead, bounds_for_bases

            alpha = 1.0 if isinstance(self.alpha, (list, tuple)) else self.alpha
            bounds = bounds_for_bases(self.core_.bases_, dict(self.monotonic_constraints))
            self.head_ = MonotonicAdditiveHead("regression", bounds, alpha=alpha)
        else:
            self.head_ = AdditiveHead(
                "regression", alpha=self.alpha, random_state=self.random_state
            )
    else:
        self.head_ = FlexibleHead(
            self.flexible_estimator, "regression", random_state=self.random_state
        )
    self.head_.fit(Z, y)
    return self

Core

Rule discovery and Z(X) basis-matrix construction.

core

Flag Core Module: rule discovery and Z(X) construction (training data only).

FlagCoreModule

FlagCoreModule(task: str, quantile_low: tuple[float, float] = (0.05, 0.45), quantile_high: tuple[float, float] = (0.55, 0.95), quantile_step: float = 0.05, min_support: int | str = 'auto', fdr_alpha: float = 0.05, effect_size: str = 'risk_difference', missing: str = 'no_evidence')

Discovers per-feature flag bases and assembles the sparse Z(X) matrix.

Source code in src/flaggam/core.py

def __init__(
    self,
    task: str,
    quantile_low: tuple[float, float] = (0.05, 0.45),
    quantile_high: tuple[float, float] = (0.55, 0.95),
    quantile_step: float = 0.05,
    min_support: int | str = "auto",
    fdr_alpha: float = 0.05,
    effect_size: str = "risk_difference",
    missing: str = "no_evidence",
) -> None:
    self.task = task
    self.quantile_low = quantile_low
    self.quantile_high = quantile_high
    self.quantile_step = quantile_step
    self.min_support = min_support
    self.fdr_alpha = fdr_alpha
    self.effect_size = effect_size
    self.missing = missing

Screening

UFA screening statistics used to select candidate rules.

screening

Screening statistics for rule discovery.

These are screening tools, not confirmatory inference (Zhao & Welsch, arXiv:2605.31189). Two-proportion test falls back to Fisher's exact test when any expected cell count is below 5 (documented decision, spec §6.3).

compute_min_support

compute_min_support(n_train: int) -> int

Minimum tail/level support: min(200, max(20, ceil(0.02 * n_train))).

Source code in src/flaggam/screening.py

def compute_min_support(n_train: int) -> int:
    """Minimum tail/level support: min(200, max(20, ceil(0.02 * n_train)))."""
    return min(200, max(20, math.ceil(0.02 * n_train)))

bh_adjust

bh_adjust(p_values: ndarray) -> ndarray

Benjamini-Hochberg adjusted p-values (within one feature's candidates).

Source code in src/flaggam/screening.py

def bh_adjust(p_values: np.ndarray) -> np.ndarray:
    """Benjamini-Hochberg adjusted p-values (within one feature's candidates)."""
    p = np.asarray(p_values, dtype=float)
    m = p.size
    order = np.argsort(p)
    ranked = p[order] * m / np.arange(1, m + 1)
    ranked = np.minimum.accumulate(ranked[::-1])[::-1]
    adj = np.empty(m, dtype=float)
    adj[order] = np.clip(ranked, 0.0, 1.0)
    return adj

two_proportion_test

two_proportion_test(k_tail: int, n_tail: int, k_base: int, n_base: int) -> float

Two-sided two-proportion z-test; Fisher's exact if any expected count < 5.

Source code in src/flaggam/screening.py

def two_proportion_test(k_tail: int, n_tail: int, k_base: int, n_base: int) -> float:
    """Two-sided two-proportion z-test; Fisher's exact if any expected count < 5."""
    if n_tail == 0 or n_base == 0:
        return 1.0
    if _expected_counts(k_tail, n_tail, k_base, n_base) < 5.0:
        table = [[k_tail, n_tail - k_tail], [k_base, n_base - k_base]]
        return float(stats.fisher_exact(table, alternative="two-sided")[1])
    p1, p0 = k_tail / n_tail, k_base / n_base
    pooled = (k_tail + k_base) / (n_tail + n_base)
    se = math.sqrt(pooled * (1.0 - pooled) * (1.0 / n_tail + 1.0 / n_base))
    if se == 0.0:
        return 1.0
    z = (p1 - p0) / se
    return float(2.0 * stats.norm.sf(abs(z)))

chi_square_test

chi_square_test(y_tail: ndarray, y_base: ndarray) -> float

Pearson chi-square p-value on the 2 x K (region x class) table.

Source code in src/flaggam/screening.py

def chi_square_test(y_tail: np.ndarray, y_base: np.ndarray) -> float:
    """Pearson chi-square p-value on the 2 x K (region x class) table."""
    classes = np.union1d(np.unique(y_tail), np.unique(y_base))
    table = np.array(
        [[np.sum(y_tail == c) for c in classes], [np.sum(y_base == c) for c in classes]]
    )
    keep = table.sum(axis=0) > 0
    table = table[:, keep]
    if table.shape[1] < 2 or table.sum(axis=1).min() == 0:
        return 1.0
    return float(stats.chi2_contingency(table, correction=False)[1])

welch_t_test

welch_t_test(a: ndarray, b: ndarray) -> float

Two-sided Welch t-test p-value; 1.0 when a side is degenerate.

Source code in src/flaggam/screening.py

def welch_t_test(a: np.ndarray, b: np.ndarray) -> float:
    """Two-sided Welch t-test p-value; 1.0 when a side is degenerate."""
    if len(a) < 2 or len(b) < 2:
        return 1.0
    res = stats.ttest_ind(a, b, equal_var=False)
    return float(res.pvalue) if np.isfinite(res.pvalue) else 1.0

risk_difference

risk_difference(k_tail: int, n_tail: int, k_base: int, n_base: int) -> float

Absolute difference in positive-class rate, tail vs baseline (spec §6.1).

Source code in src/flaggam/screening.py

def risk_difference(k_tail: int, n_tail: int, k_base: int, n_base: int) -> float:
    """Absolute difference in positive-class rate, tail vs baseline (spec §6.1)."""
    return abs(k_tail / n_tail - k_base / n_base)

log_odds_ratio

log_odds_ratio(k_tail: int, n_tail: int, k_base: int, n_base: int) -> float

Absolute log odds ratio with 0.5 continuity correction.

Source code in src/flaggam/screening.py

def log_odds_ratio(k_tail: int, n_tail: int, k_base: int, n_base: int) -> float:
    """Absolute log odds ratio with 0.5 continuity correction."""
    a, b = k_tail + 0.5, n_tail - k_tail + 0.5
    c, d = k_base + 0.5, n_base - k_base + 0.5
    return abs(math.log((a * d) / (b * c)))

standardized_mean_difference

standardized_mean_difference(a: ndarray, b: ndarray) -> float

Absolute standardized mean difference with pooled (average) variance.

Source code in src/flaggam/screening.py

def standardized_mean_difference(a: np.ndarray, b: np.ndarray) -> float:
    """Absolute standardized mean difference with pooled (average) variance."""
    var = (np.var(a, ddof=1) + np.var(b, ddof=1)) / 2.0
    if var == 0.0:
        return 0.0
    return float(abs(np.mean(a) - np.mean(b)) / np.sqrt(var))

Bases

Basis objects — one column of Z(X) each, with screening metadata.

bases

Basis objects: one column of Z(X) each, with screening metadata.

Missing semantics ("no_evidence", spec §7): NaN/None input never triggers a basis; transform returns 0.0 there. For TrendBasis a missing value maps to 0.0, i.e. the feature mean (documented decision; paper is silent).

Basis dataclass

Basis(feature: str, support: int, effect_size: float, p_value: float, p_adj: float, enriched_class: Any = None)

One univariate basis function z_ir(x_i) plus its discovery metadata.

ThresholdBasis dataclass

ThresholdBasis(feature: str, support: int, effect_size: float, p_value: float, p_adj: float, enriched_class: Any = None, cutoff: float = 0.0, side: str = 'low')

Bases: Basis

Tail flag 1{x <= c} (side='low') or 1{x >= c} (side='high').

CategoryBasis dataclass

CategoryBasis(feature: str, support: int, effect_size: float, p_value: float, p_adj: float, enriched_class: Any = None, level: Any = None)

Bases: Basis

Level flag 1{x == v}; also the regression step basis.

HingeBasis dataclass

HingeBasis(feature: str, support: int, effect_size: float, p_value: float, p_adj: float, enriched_class: Any = None, cutoff: float = 0.0, side: str = 'low')

Bases: Basis

Tail-deviation hinge (x - c)+ (side='high') or (c - x)+ (side='low').

TrendBasis dataclass

TrendBasis(feature: str, support: int, effect_size: float, p_value: float, p_adj: float, enriched_class: Any = None, mean: float = 0.0)

Bases: Basis

Centered baseline trend x - mean (regression numerical features).

MissingIndicatorBasis dataclass

MissingIndicatorBasis(feature: str, support: int, effect_size: float, p_value: float, p_adj: float, enriched_class: Any = None)

Bases: Basis

Explicit flag 1{x is missing} (missing='indicator' mode only).

Missing

Missing-indicator discovery.

missing

Missing-indicator discovery. 'no_evidence' semantics live in bases.transform.

discover_missing_indicators

discover_missing_indicators(X: DataFrame, y: ndarray, task: str, min_support: int, fdr_alpha: float) -> list[MissingIndicatorBasis]

One candidate per feature; BH across features (one test each).

Source code in src/flaggam/missing.py

def discover_missing_indicators(
    X: pd.DataFrame, y: np.ndarray, task: str, min_support: int, fdr_alpha: float
) -> list[MissingIndicatorBasis]:
    """One candidate per feature; BH across features (one test each)."""
    candidates: list[tuple[str, int, float]] = []
    for col in X.columns:
        miss = X[col].isna().to_numpy()
        n_miss = int(miss.sum())
        if n_miss < min_support or (len(miss) - n_miss) < min_support:
            continue
        if task == "binary":
            k1, k0 = int(y[miss].sum()), int(y[~miss].sum())
            p = two_proportion_test(k1, n_miss, k0, len(miss) - n_miss)
        else:
            p = chi_square_test(y[miss], y[~miss])
        candidates.append((col, n_miss, p))
    if not candidates:
        return []
    p_adj = bh_adjust(np.array([c[2] for c in candidates]))
    out = []
    for (col, n_miss, p), pa in zip(candidates, p_adj, strict=False):
        if pa <= fdr_alpha:
            out.append(
                MissingIndicatorBasis(
                    feature=col,
                    support=n_miss,
                    effect_size=float("nan"),
                    p_value=p,
                    p_adj=float(pa),
                    enriched_class=None,
                )
            )
    return out

Heads

Prediction heads (additive / flexible) fit on Z(X).

heads

Prediction heads fit on Z(X) only. Z is used unstandardized (see DECISIONS.md).

AdditiveHead

AdditiveHead(task: str, C: float | list[float] = 1.0, alpha: float | list[float] = 1.0, random_state: int | None = None)

L2 logistic/softmax (classification) or ridge (regression) on Z(X).

Source code in src/flaggam/heads.py

def __init__(
    self,
    task: str,
    C: float | list[float] = 1.0,
    alpha: float | list[float] = 1.0,
    random_state: int | None = None,
) -> None:
    self.task = task
    self.C = C
    self.alpha = alpha
    self.random_state = random_state

FlexibleHead

FlexibleHead(estimator, task: str, random_state: int | None = None)

User-supplied tree-ensemble estimator fit on Z(X) only (no raw features).

Source code in src/flaggam/heads.py

def __init__(self, estimator, task: str, random_state: int | None = None) -> None:
    self.estimator = estimator
    self.task = task
    self.random_state = random_state

Weighting

Feature-weight statistics and the compact score representation.

weighting

Feature-weight statistics and the compact score representation (spec §6.6-6.7).

Compact scores are a classification-only ablation/option; benchmark default is the full Z(X). A flag contributes to the score of its enriched class.

point_biserial

point_biserial(x: ndarray, y: ndarray) -> float

Absolute Pearson correlation between numerical x and binary y.

Source code in src/flaggam/weighting.py

def point_biserial(x: np.ndarray, y: np.ndarray) -> float:
    """Absolute Pearson correlation between numerical x and binary y."""
    if np.std(x) == 0.0 or np.std(y) == 0.0:
        return 0.0
    return float(abs(np.corrcoef(x, y)[0, 1]))

correlation_ratio

correlation_ratio(x: ndarray, y: ndarray) -> float

Correlation ratio eta of numerical x across the groups defined by y.

Source code in src/flaggam/weighting.py

def correlation_ratio(x: np.ndarray, y: np.ndarray) -> float:
    """Correlation ratio eta of numerical x across the groups defined by y."""
    grand = x.mean()
    ss_total = float(((x - grand) ** 2).sum())
    if ss_total == 0.0:
        return 0.0
    ss_between = sum(len(x[y == c]) * (x[y == c].mean() - grand) ** 2 for c in np.unique(y))
    return float(np.sqrt(ss_between / ss_total))

cramers_v

cramers_v(x: ndarray, y: ndarray) -> float

Cramer's V between categorical x and target labels y.

Source code in src/flaggam/weighting.py

def cramers_v(x: np.ndarray, y: np.ndarray) -> float:
    """Cramer's V between categorical x and target labels y."""
    table = pd.crosstab(pd.Series(x), pd.Series(y)).to_numpy(dtype=float)
    n = table.sum()
    if n == 0 or min(table.shape) < 2:
        return 0.0
    expected = np.outer(table.sum(1), table.sum(0)) / n
    with np.errstate(divide="ignore", invalid="ignore"):
        chi2 = np.nansum((table - expected) ** 2 / expected)
    k = min(table.shape) - 1
    return float(np.sqrt(chi2 / (n * k)))

feature_weights

feature_weights(X: DataFrame, y: ndarray, task: str, numerical: list[str]) -> dict[str, float]

Per-feature association weight on non-missing training rows.

Source code in src/flaggam/weighting.py

def feature_weights(
    X: pd.DataFrame, y: np.ndarray, task: str, numerical: list[str]
) -> dict[str, float]:
    """Per-feature association weight on non-missing training rows."""
    weights: dict[str, float] = {}
    for col in X.columns:
        obs = X[col].notna().to_numpy()
        if obs.sum() < 2:
            weights[col] = 0.0
            continue
        xv, yv = X[col].to_numpy()[obs], y[obs]
        if col in numerical:
            xv = xv.astype(float)
            weights[col] = point_biserial(xv, yv) if task == "binary" else correlation_ratio(xv, yv)
        else:
            weights[col] = cramers_v(xv, yv)
    return weights

compact_scores

compact_scores(Z: spmatrix, bases: list[Basis], classes: ndarray, weights: dict[str, float] | None) -> ndarray

Per-class (optionally feature-weighted) sums of triggered flags: (n, K).

Source code in src/flaggam/weighting.py

def compact_scores(
    Z: sparse.spmatrix,
    bases: list[Basis],
    classes: np.ndarray,
    weights: dict[str, float] | None,
) -> np.ndarray:
    """Per-class (optionally feature-weighted) sums of triggered flags: (n, K)."""
    Zd = np.asarray(Z.todense(), dtype=float)
    out = np.zeros((Zd.shape[0], len(classes)))
    class_index = {c: k for k, c in enumerate(classes)}
    for j, b in enumerate(bases):
        if b.kind not in _FLAG_KINDS or b.enriched_class not in class_index:
            continue
        w = 1.0 if weights is None else weights.get(b.feature, 0.0)
        out[:, class_index[b.enriched_class]] += w * Zd[:, j]
    return out

Inspection

Rule export (export_rules) and per-row reason codes (explain).

inspection

Rule export and per-row reason codes for fitted FlagGAM estimators.

Calibration

PD calibration diagnostics and recalibration. An original addition — see Calibration.

calibration

PD calibration diagnostics and recalibration for FlagGAM.

This module is an ORIGINAL ADDITION and is not part of Zhao & Welsch (arXiv:2605.31189): the paper evaluates only ranking metrics and never assesses probability calibration. Design notes in docs/DECISIONS.md entry 19.

CalibratedFlagGAM

CalibratedFlagGAM(estimator: Any, method: str = 'platt', cv: int | str = 5, target_rate: float | None = None)

Recalibration wrapper: platt / isotonic / base_rate over a FlagGAM classifier.

cv=k cross-fits k clones to obtain leak-free out-of-fold predictions for the single pooled calibrator, then refits the estimator on all data; cv="prefit" treats estimator as already fitted and (X, y) in fit() as pure calibration data. Binary classification only.

Source code in src/flaggam/calibration.py

def __init__(
    self,
    estimator: Any,
    method: str = "platt",
    cv: int | str = 5,
    target_rate: float | None = None,
) -> None:
    self.estimator = estimator
    self.method = method
    self.cv = cv
    self.target_rate = target_rate

reliability_curve

reliability_curve(y_true: Any, y_prob: Any, n_bins: int = 10, strategy: str = 'uniform') -> DataFrame

Binned reliability table (empty bins dropped).

Source code in src/flaggam/calibration.py

def reliability_curve(
    y_true: Any, y_prob: Any, n_bins: int = 10, strategy: str = "uniform"
) -> pd.DataFrame:
    """Binned reliability table (empty bins dropped)."""
    y, p = _validate_binary(y_true, y_prob)
    edges = _bin_edges(p, n_bins, strategy)
    idx = np.clip(np.digitize(p, edges[1:-1], right=False), 0, len(edges) - 2)
    rows = []
    for b in range(len(edges) - 1):
        mask = idx == b
        if not mask.any():
            continue
        rows.append(
            {
                "bin_lower": edges[b],
                "bin_upper": edges[b + 1],
                "mean_predicted": float(p[mask].mean()),
                "fraction_positive": float(y[mask].mean()),
                "count": int(mask.sum()),
            }
        )
    return pd.DataFrame(rows)

Monotonic

Exact monotonicity constraints for the additive head. An original addition — see Monotonicity.

monotonic

Exact monotonicity constraints for the FlagGAM additive head.

This module is an ORIGINAL ADDITION and is not part of Zhao & Welsch (arXiv:2605.31189). Because each numerical feature contributes monotone step/ramp bases (tail flags, hinges, trend), sign constraints on their coefficients yield EXACT monotonicity of the additive contribution. Sign table and design notes in docs/DECISIONS.md entry 20.

MonotonicAdditiveHead

MonotonicAdditiveHead(task: str, bounds: list[tuple[float | None, float | None]], C: float = 1.0, alpha: float = 1.0)

L-BFGS-B box-constrained L2 logistic (binary) or ridge (regression) head.

Drop-in for AdditiveHead when monotonic_constraints is active. Unlike AdditiveHead, C/alpha are single floats only: CV tuning of the constrained head is out of scope (spec routes list-valued C/alpha to the spec default of 1.0 before construction, see estimator.py).

Source code in src/flaggam/monotonic.py

def __init__(
    self,
    task: str,
    bounds: list[tuple[float | None, float | None]],
    C: float = 1.0,
    alpha: float = 1.0,
) -> None:
    self.task = task
    self.bounds = bounds
    self.C = C
    self.alpha = alpha

bounds_for_bases

bounds_for_bases(bases: list, constraints: dict) -> list[tuple[float | None, float | None]]

Per-basis-column box bounds implementing spec §8.2 sign constraints.

Source code in src/flaggam/monotonic.py

def bounds_for_bases(bases: list, constraints: dict) -> list[tuple[float | None, float | None]]:
    """Per-basis-column box bounds implementing spec §8.2 sign constraints."""
    out: list[tuple[float | None, float | None]] = []
    for b in bases:
        direction = constraints.get(b.feature, 0)
        rule = _PLUS_ONE.get(b.kind)
        if direction == 0 or rule is None:  # unconstrained, categorical, or missing
            if direction != 0 and b.kind not in ("category", "missing_indicator"):
                logger.warning(
                    "basis kind %r on constrained feature %r is not in the "
                    "monotonicity sign table; leaving its coefficient unconstrained — "
                    "monotonicity is no longer guaranteed",
                    b.kind,
                    b.feature,
                )
            out.append(_FREE)
        elif direction == 1:
            out.append(rule)
        else:  # -1: mirror the +1 bound by swapping the tuple
            out.append((rule[1], rule[0]))
    return out

Fairness

Group metrics and rule-level proxy audit. An original addition — see Fairness.

fairness

Fairness diagnostics and rule-level proxy audit for FlagGAM.

This module is an ORIGINAL ADDITION and is not part of Zhao & Welsch (arXiv:2605.31189); it operationalizes the paper's own Impact-Statement warning that selected rules may encode bias or proxies for protected attributes. Thresholds and binarization notes in docs/DECISIONS.md entry 21.

ProxyAudit

ProxyAudit(estimator: Any)

Rank fitted rule bases by association with a protected attribute.

Association is computed on the BINARIZED basis indicator z > 0: exact for threshold/category/missing_indicator flag bases, and a documented approximation (fires vs. does-not-fire) for the continuous hinge/trend bases. See docs/DECISIONS.md entry 21.

Source code in src/flaggam/fairness.py

def __init__(self, estimator: Any) -> None:
    self.estimator = estimator

report

report(X: Any, A: Any, threshold: float = 0.2) -> DataFrame

One row per basis: feature, rule, kind, association, method, flagged.

Source code in src/flaggam/fairness.py

def report(self, X: Any, A: Any, threshold: float = 0.2) -> pd.DataFrame:
    """One row per basis: feature, rule, kind, association, method, flagged."""
    Zb, bases = self._indicators(X)
    a = np.asarray(A).ravel()
    numeric_a = np.issubdtype(a.dtype, np.number)
    rows = []
    for j, basis in enumerate(bases):
        z = Zb[:, j]
        if numeric_a:
            assoc = 0.0 if len(np.unique(z)) < 2 else abs(stats.pointbiserialr(z, a)[0])
            method = "point_biserial"
        else:
            assoc = _cramers_v(z, a)
            method = "cramers_v"
        rows.append(
            {
                "feature": basis.feature,
                "rule": basis.name,
                "kind": basis.kind,
                "association": assoc,
                "method": method,
                "flagged": assoc > threshold,
            }
        )
    return (
        pd.DataFrame(rows)
        .sort_values("association", ascending=False, kind="stable")
        .reset_index(drop=True)
    )

drop_proxies

drop_proxies(X: Any, y: Any, A: Any, threshold: float = 0.2) -> tuple[Any, DataFrame]

Refit the head without flagged bases; return (new_estimator, trade-off row).

Supports only fitted binary classifiers with representation="full" and head="additive": compact-score columns don't map 1:1 to bases and flexible heads can't be refit column-wise; multiclass is out of PD scope. Estimators fitted with monotonic_constraints are also rejected, since the head refit would discard the constraints (see docs/DECISIONS.md entry 21). y must be numeric 0/1, matching the group_metrics contract used for the before/after trade-off.

Source code in src/flaggam/fairness.py

def drop_proxies(
    self, X: Any, y: Any, A: Any, threshold: float = 0.2
) -> tuple[Any, pd.DataFrame]:
    """Refit the head without flagged bases; return (new_estimator, trade-off row).

    Supports only fitted binary classifiers with `representation="full"`
    and `head="additive"`: compact-score columns don't map 1:1 to bases
    and flexible heads can't be refit column-wise; multiclass is out of
    PD scope. Estimators fitted with `monotonic_constraints` are also
    rejected, since the head refit would discard the constraints
    (see docs/DECISIONS.md entry 21). `y` must be numeric 0/1, matching
    the `group_metrics` contract used for the before/after trade-off.
    """
    est = self.estimator
    n_classes = len(getattr(est, "classes_", []))
    if not (n_classes == 2 and est.representation == "full" and est.head == "additive"):
        raise ValueError(
            "drop_proxies requires a fitted binary classifier with "
            "representation='full' and head='additive'"
        )
    if est.monotonic_constraints is not None:
        raise ValueError(
            "drop_proxies does not support monotonic-constrained estimators "
            "(the head refit would discard the constraints)"
        )
    report = self.report(X, A, threshold=threshold)
    flagged_rules = set(report.loc[report.flagged, "rule"])
    y_arr = np.asarray(y).ravel()
    p_before = est.predict_proba(X)[:, 1]

    new_est = copy.deepcopy(est)
    new_est.core_.bases_ = [b for b in est.core_.bases_ if b.name not in flagged_rules]
    df = new_est._to_frame(X, reset=False)
    Z = new_est.core_.transform(df)
    y_enc = new_est.label_encoder_.transform(y_arr)
    new_est.head_ = AdditiveHead(
        "binary", C=new_est.C, random_state=new_est.random_state
    ).fit(Z, y_enc)
    p_after = new_est.predict_proba(X)[:, 1]

    trade = pd.DataFrame(
        [
            {
                "n_dropped": int(len(est.core_.bases_) - len(new_est.core_.bases_)),
                "auroc_before": float(roc_auc_score(y_arr, p_before)),
                "auroc_after": float(roc_auc_score(y_arr, p_after)),
                "dp_diff_before": group_metrics(y_arr, p_before, A)["gaps"][
                    "demographic_parity_diff"
                ],
                "dp_diff_after": group_metrics(y_arr, p_after, A)["gaps"][
                    "demographic_parity_diff"
                ],
            }
        ]
    )
    if not flagged_rules:
        logger.info("drop_proxies: nothing flagged at threshold %.3f", threshold)
    return new_est, trade

group_metrics

group_metrics(y_true: Any, y_prob: Any, A: Any, threshold: float = 0.5, n_bins: int = 10) -> dict[str, Any]

Per-group PD metrics and max-minus-min gaps for protected attribute A.

Source code in src/flaggam/fairness.py

def group_metrics(
    y_true: Any, y_prob: Any, A: Any, threshold: float = 0.5, n_bins: int = 10
) -> dict[str, Any]:
    """Per-group PD metrics and max-minus-min gaps for protected attribute A."""
    y = np.asarray(y_true).ravel()
    p = np.asarray(y_prob, dtype=float).ravel()
    a = np.asarray(A).ravel()
    if set(np.unique(y)) - {0, 1}:
        raise ValueError("group_metrics supports binary y_true only")
    rows = {}
    for g in pd.unique(a):
        m = a == g
        yg, pg = y[m], p[m]
        sel = (pg >= threshold).astype(int)
        pos = yg == 1
        single_class = len(np.unique(yg)) < 2
        rows[g] = {
            "n": int(m.sum()),
            "base_rate": float(yg.mean()),
            "mean_predicted": float(pg.mean()),
            "selection_rate": float(sel.mean()),
            "tpr": float(sel[pos].mean()) if pos.any() else np.nan,
            "auroc": np.nan if single_class else float(roc_auc_score(yg, pg)),
            "ece": expected_calibration_error(yg, pg, n_bins=n_bins),
        }
    by_group = pd.DataFrame.from_dict(rows, orient="index")

    def _gap(col: str) -> float:
        vals = by_group[col].dropna()
        return float(vals.max() - vals.min()) if len(vals) > 1 else np.nan

    gaps = {
        "demographic_parity_diff": _gap("selection_rate"),
        "equal_opportunity_diff": _gap("tpr"),
        "auroc_gap": _gap("auroc"),
    }
    return {"by_group": by_group, "gaps": gaps}

Datasets

Benchmark dataset loaders with local caching.

datasets

Benchmark dataset loaders with local caching and license notes.

Each loader returns (X, y): X a DataFrame whose categorical features are pd.Categorical dtype, y a Series (binary targets are int 0/1 with the "event"/positive class = 1). Data is fetched at runtime and cached under data_dir(); raw files are never committed to the repository. Verify each dataset's license on its source page before redistributing anything.

data_dir

data_dir() -> Path

Cache directory: $FLAGGAM_DATA_DIR or ~/.cache/flaggam (created).

Source code in src/flaggam/datasets.py

def data_dir() -> Path:
    """Cache directory: $FLAGGAM_DATA_DIR or ~/.cache/flaggam (created)."""
    d = Path(os.environ.get("FLAGGAM_DATA_DIR", Path.home() / ".cache" / "flaggam"))
    d.mkdir(parents=True, exist_ok=True)
    return d

load_breast_cancer

load_breast_cancer() -> tuple[DataFrame, Series]

Wisconsin Diagnostic Breast Cancer (569 x 30, binary; positive=malignant).

Source: https://archive.ics.uci.edu/dataset/17 (bundled via scikit-learn). License: CC BY 4.0 per UCI page — verify before redistribution.

Source code in src/flaggam/datasets.py

def load_breast_cancer() -> tuple[pd.DataFrame, pd.Series]:
    """Wisconsin Diagnostic Breast Cancer (569 x 30, binary; positive=malignant).

    Source: https://archive.ics.uci.edu/dataset/17 (bundled via scikit-learn).
    License: CC BY 4.0 per UCI page — verify before redistribution.
    """
    from sklearn.datasets import load_breast_cancer as _load

    b = _load(as_frame=True)
    X = b.data.astype(float)
    y = pd.Series((b.target == 0).astype(int), name="malignant")  # sklearn: 0=malignant
    return X, y

load_california

load_california() -> tuple[DataFrame, Series]

California Housing (20640 x 8, regression on median house value).

Source: https://scikit-learn.org/stable/datasets/real_world.html#california-housing-dataset License: public domain (US Census derived) — verify before redistribution.

Source code in src/flaggam/datasets.py

def load_california() -> tuple[pd.DataFrame, pd.Series]:
    """California Housing (20640 x 8, regression on median house value).

    Source: https://scikit-learn.org/stable/datasets/real_world.html#california-housing-dataset
    License: public domain (US Census derived) — verify before redistribution.
    """
    from sklearn.datasets import fetch_california_housing

    b = fetch_california_housing(as_frame=True, data_home=str(data_dir()))
    return b.data.astype(float), b.target.rename("median_house_value")

load_pima

load_pima() -> tuple[DataFrame, Series]

Pima Indians Diabetes (768 x 8, binary; positive=diabetic).

Source: https://www.openml.org/d/37 (Pima Indians Diabetes, dataset id=37). License: CC0 per OpenML page — verify before redistribution. Observed variant: OpenML dataset 37 (pima_diabetes), 768 rows, 8 features. Short attribute names renamed: plas→glucose, pres→blood_pressure, skin→skin_thickness, insu→insulin, mass→bmi, pedi→diabetes_pedigree.

DECISIONS 17: Physiologically impossible zeros (glucose, blood_pressure, skin_thickness, insulin, bmi) are replaced with NaN per clinical convention and FlagGAM's native-missing design.

Source code in src/flaggam/datasets.py

def load_pima() -> tuple[pd.DataFrame, pd.Series]:
    """Pima Indians Diabetes (768 x 8, binary; positive=diabetic).

    Source: https://www.openml.org/d/37 (Pima Indians Diabetes, dataset id=37).
    License: CC0 per OpenML page — verify before redistribution.
    Observed variant: OpenML dataset 37 (pima_diabetes), 768 rows, 8 features.
    Short attribute names renamed: plas→glucose, pres→blood_pressure,
    skin→skin_thickness, insu→insulin, mass→bmi, pedi→diabetes_pedigree.

    DECISIONS 17: Physiologically impossible zeros (glucose, blood_pressure,
    skin_thickness, insulin, bmi) are replaced with NaN per clinical convention
    and FlagGAM's native-missing design.
    """

    def _fetch() -> pd.DataFrame:
        import openml

        dataset = openml.datasets.get_dataset(37)
        X_raw, y_raw, _, _ = dataset.get_data(target="class")
        assert isinstance(X_raw, pd.DataFrame)
        _rename = {
            "plas": "glucose",
            "pres": "blood_pressure",
            "skin": "skin_thickness",
            "insu": "insulin",
            "mass": "bmi",
            "pedi": "diabetes_pedigree",
        }
        X_raw = X_raw.rename(columns={k: v for k, v in _rename.items() if k in X_raw.columns})
        X_raw["_target"] = y_raw
        return X_raw

    df = _cached("pima", _fetch)
    y = (df["_target"] == "tested_positive").astype(int)
    y.name = "diabetic"
    X = df.drop("_target", axis=1).copy()

    _zero_cols = ["glucose", "blood_pressure", "skin_thickness", "insulin", "bmi"]
    for col in _zero_cols:
        if col in X.columns:
            X[col] = X[col].replace(0.0, float("nan"))

    return X, y

load_heart

load_heart() -> tuple[DataFrame, Series]

Heart Disease Cleveland (303 x 13, binary; positive=disease present).

Source: https://archive.ics.uci.edu/dataset/45 (Heart Disease, Cleveland). License: CC BY 4.0 — verify before redistribution. Target: num > 0 → 1. Observed variant: UCI id=45, Cleveland 303 rows, 13 features (all numeric). ca/thal have missing values (kept as NaN). Categorical features (originally coded as integers): cp, restecg, slope, thal, sex, fbs, exang.

Source code in src/flaggam/datasets.py

def load_heart() -> tuple[pd.DataFrame, pd.Series]:
    """Heart Disease Cleveland (303 x 13, binary; positive=disease present).

    Source: https://archive.ics.uci.edu/dataset/45 (Heart Disease, Cleveland).
    License: CC BY 4.0 — verify before redistribution. Target: num > 0 → 1.
    Observed variant: UCI id=45, Cleveland 303 rows, 13 features (all numeric).
    ca/thal have missing values (kept as NaN).
    Categorical features (originally coded as integers): cp, restecg, slope,
    thal, sex, fbs, exang.
    """

    def _fetch() -> pd.DataFrame:
        from ucimlrepo import fetch_ucirepo

        dataset = fetch_ucirepo(id=45)
        X = dataset.data.features.copy()
        y = dataset.data.targets.copy()
        X["_target"] = y.iloc[:, 0]
        return X

    df = _cached("heart", _fetch)
    y = (df["_target"] > 0).astype(int)
    y.name = "heart_disease"
    X = df.drop("_target", axis=1).copy()

    _cat_cols = ["cp", "restecg", "slope", "thal", "sex", "fbs", "exang"]
    for col in _cat_cols:
        if col in X.columns:
            X[col] = pd.Categorical(X[col])

    return X, y

load_german_credit

load_german_credit() -> tuple[DataFrame, Series]

Statlog German Credit (1000 x 20, binary; positive=bad credit).

Source: https://archive.ics.uci.edu/dataset/144 (Statlog German Credit). License: CC BY 4.0 — verify before redistribution. Target: bad credit = 1. Observed variant: UCI id=144, 1000 rows, 20 features (Attribute* columns). Target: class==2 → bad credit=1 (per UCI docs, 1=good, 2=bad).

Source code in src/flaggam/datasets.py

def load_german_credit() -> tuple[pd.DataFrame, pd.Series]:
    """Statlog German Credit (1000 x 20, binary; positive=bad credit).

    Source: https://archive.ics.uci.edu/dataset/144 (Statlog German Credit).
    License: CC BY 4.0 — verify before redistribution. Target: bad credit = 1.
    Observed variant: UCI id=144, 1000 rows, 20 features (Attribute* columns).
    Target: class==2 → bad credit=1 (per UCI docs, 1=good, 2=bad).
    """

    def _fetch() -> pd.DataFrame:
        from ucimlrepo import fetch_ucirepo

        dataset = fetch_ucirepo(id=144)
        X = dataset.data.features.copy()
        y = dataset.data.targets.copy()
        X["_target"] = y.iloc[:, 0]
        return X

    df = _cached("german_credit", _fetch)
    y = (df["_target"] == 2).astype(int)
    y.name = "bad_credit"
    X = df.drop("_target", axis=1).copy()

    for col in X.select_dtypes(include=["object", "string"]).columns:
        X[col] = pd.Categorical(X[col])

    return X, y

load_adult

load_adult() -> tuple[DataFrame, Series]

Adult / Census Income (~48842 rows, binary; positive=income >50K).

Source: https://archive.ics.uci.edu/dataset/2 (Adult / Census Income). License: CC BY 4.0 — verify before redistribution. Target: income >50K = 1. Observed variant: UCI id=2, 48842 rows, 14 features. Whitespace stripped; '?' replaced with NaN; string cols → pd.Categorical.

Source code in src/flaggam/datasets.py

def load_adult() -> tuple[pd.DataFrame, pd.Series]:
    """Adult / Census Income (~48842 rows, binary; positive=income >50K).

    Source: https://archive.ics.uci.edu/dataset/2 (Adult / Census Income).
    License: CC BY 4.0 — verify before redistribution. Target: income >50K = 1.
    Observed variant: UCI id=2, 48842 rows, 14 features.
    Whitespace stripped; '?' replaced with NaN; string cols → pd.Categorical.
    """

    def _fetch() -> pd.DataFrame:
        from ucimlrepo import fetch_ucirepo

        dataset = fetch_ucirepo(id=2)
        X = dataset.data.features.copy()
        y = dataset.data.targets.copy()
        for col in X.select_dtypes(include=["object", "string"]).columns:
            X[col] = X[col].str.strip()
        y_col = y.iloc[:, 0]
        X["_target"] = y_col.str.strip() if hasattr(y_col, "str") else y_col
        return X

    df = _cached("adult", _fetch)
    df = df.dropna(subset=["_target"])
    y_raw = df["_target"]  # already stripped pre-cache
    X = df.drop("_target", axis=1).copy()

    for col in X.select_dtypes(include=["object", "string"]).columns:
        X[col] = X[col].replace("?", float("nan"))
        X[col] = pd.Categorical(X[col])

    y = y_raw.str.contains(">50K", na=False).astype(int)
    y.name = "high_income"

    return X, y

load_bank_marketing

load_bank_marketing() -> tuple[DataFrame, Series]

Bank Marketing — bank-additional-full variant (41188 rows, binary; positive=subscribed).

Source: https://www.openml.org/d/42813 (OpenML id=42813, bank-additional-full variant). License: CC BY 4.0 — verify before redistribution. Observed variant: bank-additional-full (spec §9), 41188 rows, 21 original features including socioeconomic columns emp.var.rate, cons.price.idx, cons.conf.idx, euribor3m, nr.employed. After dropping 'duration' → 19 feature columns + y. The post-call 'duration' column is always dropped (UCI recommendation, spec §9). 'unknown' is kept as a regular category level (FlagGAM/UFA design: treat as explicit category, not NaN, to preserve structure of missingness).

Source code in src/flaggam/datasets.py

def load_bank_marketing() -> tuple[pd.DataFrame, pd.Series]:
    """Bank Marketing — bank-additional-full variant (41188 rows, binary; positive=subscribed).

    Source: https://www.openml.org/d/42813 (OpenML id=42813, bank-additional-full variant).
    License: CC BY 4.0 — verify before redistribution.
    Observed variant: bank-additional-full (spec §9), 41188 rows, 21 original features
    including socioeconomic columns emp.var.rate, cons.price.idx, cons.conf.idx,
    euribor3m, nr.employed. After dropping 'duration' → 19 feature columns + y.
    The post-call 'duration' column is always dropped (UCI recommendation, spec §9).
    'unknown' is kept as a regular category level (FlagGAM/UFA design: treat
    as explicit category, not NaN, to preserve structure of missingness).
    """

    # invalidate stale bank-full cache (UCI id=222, 45211 rows)
    _stale_path = data_dir() / "bank_marketing.parquet"
    if _stale_path.exists():
        try:
            _stale_df = pd.read_parquet(_stale_path)
            if len(_stale_df) != 41188:
                _stale_path.unlink()
                logger.info(
                    "deleted stale bank_marketing.parquet (%d rows, expected 41188)",
                    len(_stale_df),
                )
        except Exception:
            _stale_path.unlink()
            logger.info("deleted unreadable bank_marketing.parquet")

    def _fetch() -> pd.DataFrame:
        import openml

        # OpenML id=42813 is 'bankmarketing' — bank-additional-full (41188×20+target).
        # id=44234 is bank-full (45211 rows); id=1461 is also bank-full.
        dataset = openml.datasets.get_dataset(42813)
        X_raw, y_raw, _, _ = dataset.get_data(target=dataset.default_target_attribute or "y")
        assert isinstance(X_raw, pd.DataFrame)
        if X_raw.shape[0] != 41188 or "euribor3m" not in X_raw.columns:
            dl = openml.datasets.list_datasets(
                data_name="bankmarketing", output_format="dataframe"
            )
            did = int(dl[dl["NumberOfInstances"] == 41188].index[0])
            dataset = openml.datasets.get_dataset(did)
            X_raw, y_raw, _, _ = dataset.get_data(target=dataset.default_target_attribute or "y")
            assert isinstance(X_raw, pd.DataFrame)
        assert "euribor3m" in X_raw.columns, "bank-additional-full variant not found"
        X_raw.drop(columns=["duration"], errors="ignore", inplace=True)
        X_raw["_target"] = y_raw
        return X_raw

    df = _cached("bank_marketing", _fetch)
    y = (df["_target"] == "yes").astype(int)
    y.name = "subscribed"
    X = df.drop("_target", axis=1).copy()

    for col in X.select_dtypes(include=["object", "string"]).columns:
        X[col] = pd.Categorical(X[col])

    return X, y

load_ames

load_ames() -> tuple[DataFrame, Series]

Ames Housing (2930 rows, regression; target = log(Sale_Price)).

Source: https://www.openml.org/search?type=data&q=ames+housing (Ames Housing). License: public domain (De Cock, 2011) — verify before redistribution. Observed variant: OpenML 'ames_housing' version 1 (id=43926), 2930 rows, 80 features. OpenML returns features already as ordered CategoricalDtype. Target is log(Sale_Price); RMSE is reported on the log scale.

Source code in src/flaggam/datasets.py

def load_ames() -> tuple[pd.DataFrame, pd.Series]:
    """Ames Housing (2930 rows, regression; target = log(Sale_Price)).

    Source: https://www.openml.org/search?type=data&q=ames+housing (Ames Housing).
    License: public domain (De Cock, 2011) — verify before redistribution.
    Observed variant: OpenML 'ames_housing' version 1 (id=43926), 2930 rows,
    80 features. OpenML returns features already as ordered CategoricalDtype.
    Target is log(Sale_Price); RMSE is reported on the log scale.
    """

    def _fetch() -> pd.DataFrame:
        import openml

        try:
            dataset = openml.datasets.get_dataset("ames_housing", version=1)
        except (openml.exceptions.OpenMLServerException, ValueError, KeyError):
            dl = openml.datasets.list_datasets(data_name="ames_housing", output_format="dataframe")
            did = int(dl.index[0]) if hasattr(dl, "index") else next(iter(dl))
            dataset = openml.datasets.get_dataset(did)

        target_col = dataset.default_target_attribute or "Sale_Price"
        X_raw, y_raw, _, _ = dataset.get_data(target=target_col)
        assert isinstance(X_raw, pd.DataFrame)
        X_raw["_target"] = y_raw
        return X_raw

    df = _cached("ames", _fetch)
    y = np.log(df["_target"].astype(float)).rename("log_sale_price")
    X = df.drop("_target", axis=1).copy()

    for col in X.select_dtypes(include=["object", "string"]).columns:
        X[col] = pd.Categorical(X[col])

    return X, y

load_wine_white

load_wine_white() -> tuple[DataFrame, Series]

Wine Quality white subset (4898 x 11, regression on quality score).

Source: https://archive.ics.uci.edu/dataset/186 (Wine Quality, white subset). License: CC BY 4.0 — verify before redistribution. Observed variant: UCI id=186 (combined red+white, 6497 rows); 'color' column found in ds.data.original — filtered to white wines → 4898 rows, 11 features. All features are float; no categorical columns.

Source code in src/flaggam/datasets.py

def load_wine_white() -> tuple[pd.DataFrame, pd.Series]:
    """Wine Quality white subset (4898 x 11, regression on quality score).

    Source: https://archive.ics.uci.edu/dataset/186 (Wine Quality, white subset).
    License: CC BY 4.0 — verify before redistribution.
    Observed variant: UCI id=186 (combined red+white, 6497 rows); 'color' column
    found in ds.data.original — filtered to white wines → 4898 rows, 11 features.
    All features are float; no categorical columns.
    """

    def _fetch() -> pd.DataFrame:
        from ucimlrepo import fetch_ucirepo

        dataset = fetch_ucirepo(id=186)
        orig = dataset.data.original.copy()
        if "color" in orig.columns:
            orig = orig[orig["color"] == "white"].copy()
            orig.drop(columns=["color"], inplace=True)
        # orig now has 11 feature cols + quality (target)
        return orig

    df = _cached("wine_white", _fetch)
    y = df["quality"].astype(float)
    y.name = "quality"
    X = df.drop("quality", axis=1).copy()

    return X, y

Plots

Matplotlib visualization helpers (optional viz extra).

plots

Matplotlib visualization helpers for fitted FlagGAM estimators and diagnostics.

This module is an ORIGINAL ADDITION and is not part of Zhao & Welsch (arXiv:2605.31189): the paper specifies no plotting API. matplotlib is an OPTIONAL dependency (pip install flaggam[viz]); importing this module never imports matplotlib eagerly. Every function calls _plt() first, which raises a helpful ImportError if matplotlib is not installed, so import flaggam and import flaggam.plots both work without matplotlib present.

plot_shape

plot_shape(estimator: Any, feature: str, ax: Axes | None = None, grid_points: int = 200) -> Axes

Plot the fitted additive contribution for one feature.

Numeric features get a step curve of contribution vs. value over the range spanned by that feature's discovered cutoffs (padded 10%), plus a rug of the cutoffs themselves. Categorical features get one bar per discovered level, height equal to that level's coefficient. Requires representation='full' (the compact head's coefficients are per-class scores, not per-basis weights).

Source code in src/flaggam/plots.py

def plot_shape(
    estimator: Any, feature: str, ax: "Axes | None" = None, grid_points: int = 200
) -> "Axes":
    """Plot the fitted additive contribution for one feature.

    Numeric features get a step curve of contribution vs. value over the
    range spanned by that feature's discovered cutoffs (padded 10%), plus a
    rug of the cutoffs themselves. Categorical features get one bar per
    discovered level, height equal to that level's coefficient. Requires
    `representation='full'` (the compact head's coefficients are per-class
    scores, not per-basis weights).
    """
    plt = _plt()
    if getattr(estimator, "representation", "full") == "compact":
        raise ValueError(
            "plot_shape requires representation='full'; under 'compact' the head "
            "coefficients are per-class scores, not per-rule weights"
        )
    if hasattr(estimator, "classes_") and len(estimator.classes_) > 2:
        raise ValueError("plot_shape supports binary classification and regression only")
    if not isinstance(estimator.head_, _ADDITIVE_HEADS):
        raise ValueError("plot_shape requires the additive head")
    bases = estimator.core_.bases_
    coef = np.ravel(estimator.head_.coef_)
    feat_bases = [(j, b) for j, b in enumerate(bases) if b.feature == feature]
    if not feat_bases:
        available = sorted({b.feature for b in bases})
        raise ValueError(
            f"no bases discovered for feature {feature!r}; features with bases: {available}"
        )
    y_label = "contribution (log-odds)" if hasattr(estimator, "classes_") else "contribution"

    if ax is None:
        _, ax = plt.subplots()

    if feature in estimator.core_.categorical_features_:
        labels = [str(getattr(b, "level", b.name)) for _, b in feat_bases]
        heights = [coef[j] for j, _ in feat_bases]
        ax.bar(labels, heights)
    else:
        anchors = [
            v
            for _, b in feat_bases
            for v in (getattr(b, "cutoff", None), getattr(b, "mean", None))
            if v is not None
        ]
        if not anchors:
            raise ValueError(
                f"feature {feature!r} has no threshold/hinge/trend bases to plot a shape for"
            )
        lo, hi = min(anchors), max(anchors)
        pad = 0.1 * (hi - lo) if hi > lo else max(abs(lo), 1.0)
        grid = np.linspace(lo - pad, hi + pad, grid_points)
        contribution = np.zeros_like(grid)
        for j, b in feat_bases:
            contribution += coef[j] * b.transform(grid)
        ax.step(grid, contribution, where="post")
        rug_y = ax.get_ylim()[0]
        ax.plot(anchors, [rug_y] * len(anchors), "|", color="black", markersize=12)

    ax.set_xlabel(feature)
    ax.set_ylabel(y_label)
    ax.set_title(f"{feature} shape function")
    return ax

plot_rule_importance

plot_rule_importance(estimator: Any, top_n: int = 20, ax: Axes | None = None) -> Axes

Horizontal bar chart of the top top_n rules by |weight| from export_rules().

Source code in src/flaggam/plots.py

def plot_rule_importance(estimator: Any, top_n: int = 20, ax: "Axes | None" = None) -> "Axes":
    """Horizontal bar chart of the top `top_n` rules by |weight| from `export_rules()`."""
    plt = _plt()
    if not isinstance(estimator.head_, _ADDITIVE_HEADS):
        raise ValueError("plot_rule_importance requires the additive head")
    rules = estimator.export_rules()
    top = rules.reindex(rules["weight"].abs().sort_values(ascending=False).index).head(top_n)
    if ax is None:
        _, ax = plt.subplots()
    ax.barh(top["rule"].iloc[::-1], top["weight"].iloc[::-1])
    ax.set_xlabel("weight")
    ax.set_title("Rule importance")
    return ax

plot_waterfall

plot_waterfall(estimator: Any, x_row: Any, ax: Axes | None = None, max_rules: int = 15) -> Axes

Cumulative horizontal bars from intercept to total score for one row.

Rules are sorted by |contribution| and collapsed beyond max_rules into a single "(other rules)" bucket; the final bar marks the total score.

Source code in src/flaggam/plots.py

def plot_waterfall(
    estimator: Any, x_row: Any, ax: "Axes | None" = None, max_rules: int = 15
) -> "Axes":
    """Cumulative horizontal bars from intercept to total score for one row.

    Rules are sorted by |contribution| and collapsed beyond `max_rules` into
    a single "(other rules)" bucket; the final bar marks the total score.
    """
    plt = _plt()
    exp = estimator.explain(x_row)
    row0 = exp[exp["row"] == 0]
    intercept = float(row0.loc[row0["feature"] == "<intercept>", "contribution"].sum())
    rules = row0[row0["feature"] != "<intercept>"]
    rules = rules.reindex(rules["contribution"].abs().sort_values(ascending=False).index)
    if len(rules) > max_rules:
        head, tail = rules.iloc[:max_rules], rules.iloc[max_rules:]
        names = [*head["rule"], "(other rules)"]
        contributions = [*head["contribution"].astype(float), float(tail["contribution"].sum())]
    else:
        names = list(rules["rule"])
        contributions = list(rules["contribution"].astype(float))

    labels = ["<intercept>", *names, "total"]
    values = [intercept, *contributions]
    running = 0.0
    lefts: list[float] = []
    widths: list[float] = []
    for v in values:
        lefts.append(min(running, running + v))
        widths.append(abs(v))
        running += v
    total = running
    lefts.append(min(0.0, total))
    widths.append(abs(total))

    if ax is None:
        _, ax = plt.subplots()
    y_pos = np.arange(len(labels))
    colors = (
        ["tab:blue"] + ["tab:green" if v >= 0 else "tab:red" for v in values[1:]] + ["black"]
    )
    ax.barh(y_pos, widths, left=lefts, color=colors)
    ax.set_yticks(y_pos)
    ax.set_yticklabels(labels)
    ax.invert_yaxis()
    ax.set_xlabel("contribution")
    ax.set_title(f"Waterfall (total = {total:.3f})")
    return ax

plot_reliability

plot_reliability(y_true: Any, y_prob: Any, n_bins: int = 10, strategy: str = 'uniform', ax: Axes | None = None) -> Axes

Reliability diagram: mean predicted vs. observed rate, with per-bin counts.

Source code in src/flaggam/plots.py

def plot_reliability(
    y_true: Any,
    y_prob: Any,
    n_bins: int = 10,
    strategy: str = "uniform",
    ax: "Axes | None" = None,
) -> "Axes":
    """Reliability diagram: mean predicted vs. observed rate, with per-bin counts."""
    plt = _plt()
    table = reliability_curve(y_true, y_prob, n_bins=n_bins, strategy=strategy)
    if ax is None:
        _, ax = plt.subplots()
    ax.plot([0.0, 1.0], [0.0, 1.0], linestyle="--", color="gray", label="perfectly calibrated")
    ax.plot(table["mean_predicted"], table["fraction_positive"], marker="o", label="model")
    count_ax = ax.twinx()
    count_ax.bar(
        table["mean_predicted"], table["count"], width=1.0 / n_bins, alpha=0.2, color="gray"
    )
    count_ax.set_ylabel("count")
    ax.set_xlabel("mean predicted probability")
    ax.set_ylabel("fraction of positives")
    ax.set_title("Reliability diagram")
    ax.legend()
    return ax

plot_proxy_association

plot_proxy_association(report: DataFrame, top_n: int = 20, ax: Axes | None = None) -> Axes

Horizontal bars of rule-level proxy association, flagged rules highlighted.

Source code in src/flaggam/plots.py

def plot_proxy_association(
    report: pd.DataFrame, top_n: int = 20, ax: "Axes | None" = None
) -> "Axes":
    """Horizontal bars of rule-level proxy association, flagged rules highlighted."""
    plt = _plt()
    top = report.reindex(report["association"].sort_values(ascending=False).index).head(top_n)
    if ax is None:
        _, ax = plt.subplots()
    colors = ["tab:red" if f else "tab:blue" for f in top["flagged"]]
    ax.barh(top["rule"].iloc[::-1], top["association"].iloc[::-1], color=list(reversed(colors)))
    if top["flagged"].any():
        # Exact audit threshold isn't carried in `report`; approximate the
        # boundary with the smallest association among the flagged rows.
        line_x = float(top.loc[top["flagged"], "association"].min())
        ax.axvline(line_x, linestyle="--", color="black")
    ax.set_xlabel("association")
    ax.set_title("Proxy association by rule")
    return ax

plot_group_metrics

plot_group_metrics(metrics: dict[str, Any], ax: Axes | None = None) -> Axes

Grouped bar chart of selection_rate/tpr/auroc per protected-attribute group.

Source code in src/flaggam/plots.py

def plot_group_metrics(metrics: dict[str, Any], ax: "Axes | None" = None) -> "Axes":
    """Grouped bar chart of selection_rate/tpr/auroc per protected-attribute group."""
    plt = _plt()
    by = metrics["by_group"][["selection_rate", "tpr", "auroc"]]
    if ax is None:
        _, ax = plt.subplots()
    groups = list(by.index)
    cols = list(by.columns)
    n_metrics = len(cols)
    x = np.arange(len(groups))
    width = 0.8 / n_metrics
    for j, col in enumerate(cols):
        ax.bar(x + j * width, by[col].to_numpy(dtype=float), width, label=col)
    ax.set_xticks(x + width * (n_metrics - 1) / 2)
    ax.set_xticklabels([str(g) for g in groups])
    gap_str = ", ".join(
        f"{k}={v:.3f}" if not np.isnan(v) else f"{k}=nan" for k, v in metrics["gaps"].items()
    )
    ax.set_title(f"Group metrics (gaps: {gap_str})")
    ax.legend()
    return ax

Explorer

Self-contained interactive HTML rules explorer.

explorer

Self-contained interactive HTML rules explorer for fitted FlagGAM estimators.

This module is an ORIGINAL ADDITION (not part of Zhao & Welsch, arXiv:2605.31189): the paper specifies no explorer/export API. export_rules_html renders a single HTML document with inlined CSS and vanilla JS and NO external resources (no CDN, fonts, or images), so the result is fully self-contained: it works offline and can be embedded in an iframe (e.g. the docs site) with no network access.

export_rules_html

export_rules_html(estimator: Any, path: str | Path | None = None, title: str = 'FlagGAM rules explorer') -> str

Render a self-contained interactive HTML explorer of the fitted rules.

Requires an estimator fitted with representation='full' and the additive head (binary classification or regression only; the compact head's coefficients are per-class scores, not per-rule weights). Returns the complete HTML document as a string; if path is given, also writes it (UTF-8) and still returns the string.

Source code in src/flaggam/explorer.py

def export_rules_html(
    estimator: Any, path: "str | Path | None" = None, title: str = "FlagGAM rules explorer"
) -> str:
    """Render a self-contained interactive HTML explorer of the fitted rules.

    Requires an estimator fitted with `representation='full'` and the additive
    head (binary classification or regression only; the compact head's
    coefficients are per-class scores, not per-rule weights). Returns the
    complete HTML document as a string; if `path` is given, also writes it
    (UTF-8) and still returns the string.
    """
    check_is_fitted(estimator, "core_")
    if getattr(estimator, "representation", "full") == "compact":
        raise ValueError(
            "export_rules_html requires representation='full'; under 'compact' the head "
            "coefficients are per-class scores, not per-rule weights"
        )
    if hasattr(estimator, "classes_") and len(estimator.classes_) > 2:
        raise ValueError("export_rules_html supports binary classification and regression only")
    if not isinstance(estimator.head_, _ADDITIVE_HEADS):
        raise ValueError("export_rules_html requires the additive head")

    payload = _build_payload(estimator, title)
    data = json.dumps(payload, allow_nan=False).replace("</", "<\\/")
    doc = _TEMPLATE.replace("%%DATA%%", data).replace("%%TITLE%%", html.escape(title, quote=True))

    if path is not None:
        Path(path).write_text(doc, encoding="utf-8")
    return doc