ICL-MoE Regression

This example demonstrates how to use ICL-MoE to post-process a fitted DirectRS regression model. ICL-MoE adds a kNN-based residual correction on top of DirectRS leaf experts, provi

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

import warnings
warnings.filterwarnings("ignore")

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

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

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

# %%
# Train Base Model
# ----------------
# Train an XGBoost regressor with depth 2 (required for FANOVA comparison).

model = MoXGBRegressor(
    name="XGB-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 model with DirectRS.

from modeva.models import MoDirectRSRegressor

drs = MoDirectRSRegressor(
    base_model=model, ridge_alpha=100.0, n_passes=1
)
drs.fit(
    ds.train_x, ds.train_y.ravel(),
    X_val=ds.test_x, y_val=ds.test_y.ravel(),
    verbose=True
)

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

from modeva.models import MoDirectRSICLRegressor

icl = MoDirectRSICLRegressor(
    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 R-squared and MAE between XGBoost, DirectRS, and ICL-MoE.

from sklearn.metrics import r2_score, mean_absolute_error

y_train = ds.train_y.ravel()
y_test = ds.test_y.ravel()

print(f"{'Model':<12s} {'Train R2':>10s} {'Test R2':>10s} {'Train MAE':>10s} {'Test MAE':>10s}")
print("-" * 54)
for name, m in [("XGBoost", model), ("DirectRS", drs), ("ICL-MoE", icl)]:
    tr = m.predict(ds.train_x)
    te = m.predict(ds.test_x)
    print(f"{name:<12s} {r2_score(y_train, tr):>10.4f} {r2_score(y_test, te):>10.4f} "
          f"{mean_absolute_error(y_train, tr):>10.4f} {mean_absolute_error(y_test, te):>10.4f}")

# %%
# Local Explanation
# -----------------
# ICL-MoE provides exact additive decomposition: f(x) = c_0 + sum_j c_j.
# We verify the decomposition matches predictions to machine precision.

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

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

print(f"Max |predict - (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':<12s} {'ICL-MoE':>10s} {'DirectRS':>10s} {'FANOVA':>10s}")
print("-" * 44)
for i, feat in enumerate(ds.feature_names):
    print(f"{feat:<12s} {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()