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.0181, 0.0084, 0.0073, 0.0066, 0.0061]
[DirectRS] Initial train LogLoss=0.5689, AUC=0.7978
Pass 1/1: train LogLoss=0.4143, AUC=0.8044
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.7834 0.7835
Accuracy 0.8292 0.8275
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)|: 7.19e-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.7206 (72.1%)
eta2_int = 0.2794 (27.9%)
rho(g, r) = -0.9877
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.612757
k=2: T_k = 0.612757, E_k = 0.075537
k=3: T_k = 0.098694, E_k = 0.015026
k=4: T_k = 0.075537, E_k = 0.003434
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.1120
SEX 0.0028 0.0038
EDUCATION 0.0095 0.0270
MARRIAGE 0.0078 0.0164
AGE 0.0405 0.0067
PAY_1 0.2842 0.5755
PAY_2 0.0874 0.0265
PAY_3 0.0498 0.0276
PAY_4 0.0332 0.0084
PAY_5 0.0246 0.0072
PAY_6 0.0752 0.0111
BILL_AMT1 0.1419 0.0345
BILL_AMT2 0.0144 0.0083
BILL_AMT3 0.0230 0.0142
BILL_AMT4 0.0195 0.0202
BILL_AMT5 0.0060 0.0028
BILL_AMT6 0.0051 0.0036
PAY_AMT1 0.0604 0.0227
PAY_AMT2 0.0566 0.0362
PAY_AMT3 0.0188 0.0181
PAY_AMT4 0.0126 0.0067
PAY_AMT5 0.0031 0.0015
PAY_AMT6 0.0235 0.0089
Total running time of the script: (0 minutes 12.619 seconds)