ICL-MoE Classification

This example demonstrates ICL-MoE post-processing on a fitted DirectRS classifier. ICL-MoE operates in logit space internally, adding kNN-based residual correction on top of Direct

# %%
# Setup
# -----
# Import libraries and suppress warnings.

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from modeva import DataSet, TestSuite
from modeva.models import MoXGBClassifier

# %%
# Load Dataset
# ------------
# Load the TaiwanCredit dataset and create a random train/test split.

ds = DataSet()
ds.load(name="TaiwanCredit")
ds.set_random_split()

# %%
# Train Base Model
# ----------------
# Train an XGBoost classifier with depth 2.

model = MoXGBClassifier(
    name="XGB-cls-depth2",
    n_estimators=200, max_depth=2, learning_rate=0.1,
    subsample=0.8, colsample_bytree=0.8,
    random_state=42, verbosity=0
)
model.fit(ds.train_x, ds.train_y.ravel())

ts = TestSuite(ds, model)
ts.diagnose_accuracy_table().table

# %%
# Fit DirectRS
# ------------
# Post-process the trained XGBoost classifier with DirectRS.

from modeva.models import MoDirectRSClassifier

drs = MoDirectRSClassifier(
    base_model=model, ridge_alpha=100.0, n_passes=1
)
drs.fit(ds.train_x, ds.train_y.ravel(), verbose=True)

# %%
# Fit ICL-MoE
# -----------
# Build ICL-MoE on top of the fitted DirectRS classifier. The hierarchical
# variant combines leaf expert predictions with a kNN residual correction.

from modeva.models import MoDirectRSICLClassifier

icl = MoDirectRSICLClassifier(
    directrs_model=drs, k=50, tau=1.0, ridge_lambda=1.0
)
icl.fit(ds.train_x, ds.train_y.ravel(), verbose=True)

# %%
# Accuracy Comparison
# -------------------
# Compare AUC and accuracy between XGBoost, DirectRS, and ICL-MoE.

from sklearn.metrics import roc_auc_score, accuracy_score, log_loss

y_test = ds.test_y.ravel()

base_proba = model.predict_proba(ds.test_x)[:, 1]
drs_proba = drs.predict_proba(ds.test_x)[:, 1]
icl_proba = icl.predict_proba(ds.test_x)[:, 1]

print(f"{'Metric':<10s} {'XGBoost':>10s} {'DirectRS':>10s} {'ICL-MoE':>10s}")
print("-" * 42)
print(f"{'AUC':<10s} {roc_auc_score(y_test, base_proba):>10.4f} "
      f"{roc_auc_score(y_test, drs_proba):>10.4f} "
      f"{roc_auc_score(y_test, icl_proba):>10.4f}")
print(f"{'Accuracy':<10s} {accuracy_score(y_test, model.predict(ds.test_x)):>10.4f} "
      f"{accuracy_score(y_test, drs.predict(ds.test_x)):>10.4f} "
      f"{accuracy_score(y_test, icl.predict(ds.test_x)):>10.4f}")
print(f"{'LogLoss':<10s} {log_loss(y_test, np.clip(base_proba, 1e-7, 1-1e-7)):>10.4f} "
      f"{log_loss(y_test, np.clip(drs_proba, 1e-7, 1-1e-7)):>10.4f} "
      f"{log_loss(y_test, np.clip(icl_proba, 1e-7, 1-1e-7)):>10.4f}")

# %%
# Probability Sanity Check
# ------------------------
# Verify ICL-MoE probabilities are valid (in [0, 1], rows sum to 1).

proba = icl.predict_proba(ds.test_x)

print(f"Shape: {proba.shape}")
print(f"Range: [{proba.min():.6f}, {proba.max():.6f}]")
print(f"Row sums: [{proba.sum(axis=1).min():.10f}, {proba.sum(axis=1).max():.10f}]")

# %%
# Local Explanation
# -----------------
# For classification, the decomposition operates on raw scores (logits).
# We verify using the internal core which returns raw logits.

result = icl.explain_local(ds.test_x, feature_names=ds.feature_names)
local = result.value

raw_logits = icl._icl_core.predict(ds.test_x)
recon = local['intercept'] + local['contributions'].sum(axis=1)
max_err = np.max(np.abs(raw_logits - recon))

print(f"Max |logit - (intercept + sum contributions)|: {max_err:.2e}")
print(f"Decomposition exact to machine precision: {max_err < 1e-10}")

# %%
# Local explanation waterfall plot for sample 0.
result.plot()

# %%
# Global Feature Importance
# -------------------------
# Compute global feature importance using mean absolute contributions.

result = icl.importance_global(
    ds.test_x, feature_names=ds.feature_names, mode="contrib_abs"
)
result.plot()

# %%
# Importance Comparison
# ---------------------
# Compare ICL-MoE, DirectRS, and FANOVA feature importance side-by-side.

icl_imp = icl.importance_global(
    ds.test_x, feature_names=ds.feature_names, mode="contrib_abs"
).value['importance']

drs_imp = drs.importance_global(
    ds.test_x, feature_names=ds.feature_names, mode="contrib_abs"
).value['importance']

result_fanova = ts.interpret_fi()
fanova_table = result_fanova.table
fanova_imp = dict(zip(fanova_table["Name"], fanova_table["Score"]))

print(f"{'Feature':<16s} {'ICL-MoE':>10s} {'DirectRS':>10s} {'FANOVA':>10s}")
print("-" * 48)
for i, feat in enumerate(ds.feature_names):
    print(f"{feat:<16s} {icl_imp[i]:>10.4f} {drs_imp[i]:>10.4f} {fanova_imp.get(feat, 0.0):>10.4f}")

# %%
# FANOVA feature importance plot.
result_fanova.plot()

# %%
# FANOVA interaction importance plot.
result_ei = ts.interpret_ei()
result_ei.plot()

# %%
# Neighbour Analysis
# ------------------
# Inspect the kNN neighbourhood used for residual correction.

result = icl.get_neighbor_analysis(
    ds.test_x, sample_index=0, feature_names=ds.feature_names
)

# %%
# Neighbour weights.
result.plot("weight_bar")

# %%
# PCA scatter of neighbours in z-space.
result.plot("neighbor_scatter")

# %%
# Leaf Gating Analysis
# --------------------
# Examine per-leaf routing for a specific tree.

result = icl.get_leaf_gating_analysis(
    ds.test_x, sample_index=0, tree_index=0,
    feature_names=ds.feature_names
)

active = result.value['active_leaf']
print(f"Active leaf: {active}")
print(f"Tree 0 contribution: {result.value['prediction']:.4f}")
print(f"Number of leaves: {len(result.value['leaf_weights'])}")

# %%
# Leaf gating weights for tree 0.
result.plot()