Average vs max rank on tied data

When the target y contains ties (e.g. a discretised or rounded target), rank_method="average" and rank_method="max" rank the tied target values differently and can therefore select different features. This example compares the two methods’ test R² on the seeds where they disagree; see Average vs. Max Ranking on Tied Targets for a discussion of why neither method is uniformly better.

By default only the fixed SEEDS_OF_INTEREST list is evaluated, keeping the documentation build fast. Set SEARCH_SEEDS = True to rescan the first SEARCH_SEEDS_LIMIT seeds and, if necessary, print a refreshed list in Python syntax to paste back.

Seeds compared: 34
  average better test R² : 18
  max better test R²     : 16
  equal test R²          : 0
  mean R² difference (average - max): +0.0314

import matplotlib.pyplot as plt
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split

from pyFOCI import FOCISelector

# -- config ----------------------------------------------------------------
N_SAMPLES = 350
N_FEATURES = 20
N_INFORMATIVE = 5
N_LEVELS = 5
NOISE_SIGMA = 0.5
TRAIN_FRACTION = 0.75
K_FEAT_CAP = 10

# Set SEARCH_SEEDS = True to re-run the seed search (slow) and print the
# seeds where the two rank methods disagree if they differ from the current
# SEEDS_OF_INTEREST below; paste the result back into SEEDS_OF_INTEREST.
SEARCH_SEEDS = False
SEARCH_SEEDS_LIMIT = 500

# Seeds in [0, SEARCH_SEEDS_LIMIT) where rank_method="average" and
# rank_method="max" select different feature sets, found with SEARCH_SEEDS=True.
# fmt: off
SEEDS_OF_INTEREST = [
    0, 8, 10, 16, 49, 71, 86, 112, 115, 119, 158, 173, 195, 198, 210, 214, 217,
    218, 219, 225, 232, 271, 288, 302, 304, 309, 331, 336, 350, 354, 380, 386,
    393, 471
]
# fmt: on


def make_data(seed, n=N_SAMPLES, p=N_FEATURES, n_levels=N_LEVELS, sigma=NOISE_SIGMA):
    rng = np.random.RandomState(seed)
    X = rng.normal(size=(n, p))
    y_lat = (
        2.0 * (X[:, 0] ** 2 - 1.0)
        + 1.5 * np.sin(2.0 * X[:, 1])
        + 2.0 * np.exp(-X[:, 2] ** 2)
        + 1.5 * X[:, 3] * X[:, 4]
        + 1.0 * (X[:, 3] >= 0)
    )
    y_lat += sigma * rng.normal(size=n)
    q = np.quantile(y_lat, np.linspace(0, 1, n_levels + 1))
    q[0] -= 1e-9
    q[-1] += 1e-9
    y = np.digitize(y_lat, q[1:-1]).astype(float)
    return X, y


def select(rank_method, X_train, y_train):
    """Fit FOCISelector with the given rank method and return selected indices."""
    sel = FOCISelector(
        method="r_foci",
        max_features=K_FEAT_CAP,
        min_delta=0,
        nn_tie_breaking="mean",
        nn_strategy="grouping",
        standardize="normalize",
        rank_method=rank_method,
        random_state=0,
    )
    return sel.fit(X_train, y_train).selected_indices_


def test_r2(idx, X_train, X_test, y_train, y_test):
    """Test R² of the downstream model fitted on the selected columns."""
    if len(idx) == 0:
        return 0.0
    pred = RandomForestRegressor(n_estimators=100, random_state=0, n_jobs=-1)
    pred.fit(X_train[:, idx], y_train)
    return r2_score(y_test, pred.predict(X_test[:, idx]))


def evaluate_seed(seed):
    """Return (r2_max, r2_avg) if the two methods disagree, else None."""
    X, y = make_data(seed)
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, train_size=TRAIN_FRACTION, random_state=seed
    )
    idx_m = select("max", X_train, y_train)
    idx_a = select("average", X_train, y_train)
    if set(idx_m) == set(idx_a):
        return None
    r2_m = test_r2(idx_m, X_train, X_test, y_train, y_test)
    r2_a = test_r2(idx_a, X_train, X_test, y_train, y_test)
    return r2_m, r2_a


# -- evaluate the seeds of interest -----------------------------------------
if SEARCH_SEEDS:
    seeds_to_scan = range(SEARCH_SEEDS_LIMIT)
else:
    seeds_to_scan = SEEDS_OF_INTEREST

results = []
for seed in seeds_to_scan:
    outcome = evaluate_seed(seed)
    if outcome is not None:
        r2_m, r2_a = outcome
        results.append((seed, r2_m, r2_a))

if SEARCH_SEEDS:
    found = [seed for seed, _, _ in results]
    if found != SEEDS_OF_INTEREST:
        print("Seeds where the two rank methods disagree:")
        print(f"SEEDS_OF_INTEREST = {found}")
        print()
    else:
        print("Found seeds match SEEDS_OF_INTEREST; nothing to update.")
        print()

if not results:
    raise SystemExit("No seeds with differing selections were found.")

seeds = np.array([seed for seed, _, _ in results])
r2_m = np.array([r2_m for _, r2_m, _ in results])
r2_a = np.array([r2_a for _, _, r2_a in results])
dr2 = r2_a - r2_m

# -- plot ------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(6, 5))

tags = np.where(dr2 > 1e-9, "avg", np.where(dr2 < -1e-9, "max", "same"))
colors = {"avg": "tab:green", "max": "tab:red", "same": "tab:gray"}
labels = {"avg": "average better", "max": "max better", "same": "equal"}

for tag in ("avg", "max", "same"):
    mask = tags == tag
    if mask.any():
        ax.scatter(
            r2_m[mask],
            r2_a[mask],
            c=colors[tag],
            s=50,
            alpha=0.9,
            label=labels[tag],
            edgecolor="k",
            linewidth=0.4,
        )

lo = min(r2_m.min(), r2_a.min()) - 0.03
hi = max(r2_m.max(), r2_a.max()) + 0.03
ax.plot([lo, hi], [lo, hi], "k--", lw=0.8, label="equal R²")
ax.set_xlabel("Test R², rank='max'")
ax.set_ylabel("Test R², rank='average'")
ax.set_title("Seeds where average and max ranking disagree")
ax.set_xlim(lo, hi)
ax.set_ylim(lo, hi)
ax.legend(loc="lower right", fontsize=8)
ax.set_aspect("equal")

fig.tight_layout()
plt.show()

# -- numeric summary -------------------------------------------------------
n_avg = int((dr2 > 1e-9).sum())
n_max = int((dr2 < -1e-9).sum())
n_tie = int((np.abs(dr2) <= 1e-9).sum())

print(f"Seeds compared: {len(results)}")
print(f"  average better test R² : {n_avg}")
print(f"  max better test R²     : {n_max}")
print(f"  equal test R²          : {n_tie}")
print(f"  mean R² difference (average - max): {dr2.mean():+.4f}")