Note
Go to the end to download the full example code.
DirectRS Classification
This example demonstrates DirectRS post-processing on a pre-trained XGBoost classifier. DirectRS works on the raw score (logit) space, providing exact additive decomposition of the decision function while improving or maintaining classification accuracy.
We use the TaiwanCredit dataset with a depth-2 XGBoost base model.
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,
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)
[DirectRS] construction=C, trees=200, alpha=100.0, logistic=True
S' eigenvalues (C): [0.0178, 0.0083, 0.0075, 0.0064, 0.0058]
[DirectRS] Initial train LogLoss=0.5977, AUC=0.7972
Pass 1/1: train LogLoss=0.4148, AUC=0.8034
Accuracy Comparison
Compare AUC and accuracy between the base XGBoost model and DirectRS.
from sklearn.metrics import roc_auc_score, accuracy_score
y_test = ds.test_y.ravel()
base_proba = model.predict_proba(ds.test_x)[:, 1]
drs_proba = drs.predict_proba(ds.test_x)[:, 1]
base_acc = accuracy_score(y_test, model.predict(ds.test_x))
drs_acc = accuracy_score(y_test, drs.predict(ds.test_x))
print(f"{'Metric':<10s} {'Base XGB':>10s} {'DirectRS':>10s}")
print("-" * 32)
print(f"{'AUC':<10s} {roc_auc_score(y_test, base_proba):>10.4f} {roc_auc_score(y_test, drs_proba):>10.4f}")
print(f"{'Accuracy':<10s} {base_acc:>10.4f} {drs_acc:>10.4f}")
Metric Base XGB DirectRS
--------------------------------
AUC 0.7828 0.7834
Accuracy 0.8297 0.8277
S’ Stretch Analysis
Analyze the global stretch matrix S’ extracted from tree geometry.
result = drs.get_global_stretch_analysis(ds.feature_names)
Eigenvalue spectrum of S’.
result.plot("eigenvalue_spectrum")
Feature activity scores.
result.plot("feature_activity")
Local Explanation
For classification, the decomposition operates on raw scores (logits).
We verify using drs._core.predict which returns raw scores.
result = drs.explain_local(ds.test_x, feature_names=ds.feature_names)
local = result.value
raw_pred = drs._core.predict(ds.test_x)
recon = local['intercept'] + local['contributions'].sum(axis=1)
max_err = np.max(np.abs(raw_pred - recon))
print(f"Max |raw_score - (intercept + sum contributions)|: {max_err:.2e}")
print(f"Decomposition exact to machine precision: {max_err < 1e-10}")
Max |raw_score - (intercept + sum contributions)|: 9.73e-14
Decomposition exact to machine precision: True
Local explanation waterfall plot.
result.plot()
Global Feature Importance
Compute global feature importance using the default slope mode.
result = drs.importance_global(feature_names=ds.feature_names)
result.plot()
Main/Interaction Decomposition
Decompose model variance into main effects and interactions.
result = drs.importance_main_interaction(ds.test_x, feature_names=ds.feature_names)
mi = result.value
print(f"Orthogonalized variance split:")
print(f" eta2_main = {mi['eta2_main']:.4f} ({mi['eta2_main']*100:.1f}%)")
print(f" eta2_int = {mi['eta2_int']:.4f} ({mi['eta2_int']*100:.1f}%)")
print(f" rho(g, r) = {mi['rho']:.4f}")
Orthogonalized variance split:
eta2_main = 0.7213 (72.1%)
eta2_int = 0.2787 (27.9%)
rho(g, r) = -0.9878
Main vs interaction importance bar chart.
result.plot()
Geometric Interaction Traces
Trace feature interactions through the adjacency matrix A.
result = drs.geometric_interaction_traces(
feature_names=ds.feature_names, K=4, gamma=0.5
)
traces = result.value
print("Interaction spectrum:")
for k in range(len(traces['T'])):
print(f" k={k+1}: T_k = {traces['T'][k]:.6f}, E_k = {traces['E'][k]:.6f}")
Interaction spectrum:
k=1: T_k = 0.000000, E_k = 0.717975
k=2: T_k = 0.717975, E_k = 0.087608
k=3: T_k = 0.102000, E_k = 0.017712
k=4: T_k = 0.087608, E_k = 0.004206
Adjacency matrix heatmap.
result.plot("adjacency")
Interaction energy spectrum.
result.plot("spectrum")
Top Feature Interactions
Show the strongest off-diagonal entries in the stretch matrix.
result = drs.get_off_diagonal_analysis(ds.feature_names, top_k=15)
result.plot()
FANOVA Comparison
Compare DirectRS feature importance with FANOVA decomposition.
results = ts.interpret_fi()
results.plot()
results = ts.interpret_ei()
results.plot()
Side-by-side comparison of importance scores.
result_drs = drs.importance_global(feature_names=ds.feature_names, mode="slope")
result_fanova = ts.interpret_fi()
drs_imp = result_drs.value['importance']
fanova_table = result_fanova.table
fanova_imp = dict(zip(fanova_table["Name"], fanova_table["Score"]))
print(f"{'Feature':<16s} {'DirectRS':>10s} {'FANOVA':>10s}")
print("-" * 38)
for i, feat in enumerate(ds.feature_names):
print(f"{feat:<16s} {drs_imp[i]:>10.4f} {fanova_imp.get(feat, 0.0):>10.4f}")
Feature DirectRS FANOVA
--------------------------------------
LIMIT_BAL 0.0001 0.0847
SEX 0.0026 0.0046
EDUCATION 0.0090 0.0237
MARRIAGE 0.0059 0.0130
AGE 0.0386 0.0065
PAY_1 0.2744 0.5826
PAY_2 0.0905 0.0254
PAY_3 0.0531 0.0294
PAY_4 0.0351 0.0055
PAY_5 0.0361 0.0080
PAY_6 0.0526 0.0101
BILL_AMT1 0.1588 0.0395
BILL_AMT2 0.0111 0.0075
BILL_AMT3 0.0228 0.0194
BILL_AMT4 0.0199 0.0165
BILL_AMT5 0.0072 0.0029
BILL_AMT6 0.0072 0.0050
PAY_AMT1 0.0549 0.0315
PAY_AMT2 0.0629 0.0432
PAY_AMT3 0.0198 0.0208
PAY_AMT4 0.0115 0.0086
PAY_AMT5 0.0036 0.0015
PAY_AMT6 0.0222 0.0103
Total running time of the script: (0 minutes 8.655 seconds)