Model with conditional interaction
- Author: givasile
- Runtime: ~5 s
- Description: Global effects (PDP, ALE, RHALE) on a model with a
conditional interaction; each estimate is compared β and tested β against
its analytical formula from
effector.benchmarks. - π The whole notebook in one page: PDP report
In this example, we show global effects of a model with conditional interactions using PDP, ALE, and RHALE. In particular, we:
- show how to use
effectorto estimate the global effects using PDP, ALE, and RHALE - provide the analytical formulas for the global effects
- test that (1) and (2) match
We will use the following model:
where the features \(x_1, x_2, x_3\) are independent and uniformly distributed in the interval \([-1, 1]\).
The model has an interaction between \(x_1\) and \(x_2\) caused by the terms: \(f_{1,2}(x_1, x_2) = -x_1^2 \mathbb{1}_{x_2 <0} + x_1^2 \mathbb{1}_{x_2 \geq 0}\). This means that the effect of \(x_1\) on the output \(y\) depends on the value of \(x_2\) and vice versa. Therefore, there is no golden standard on how to split the effect of \(f_{1,2}\) to two parts, one that corresponds to \(x_1\) and one to \(x_2\). Each global effect method has a different strategy to handle this issue. Below we will see how PDP, ALE, and RHALE handle this interaction.
In contrast, \(x_3\) does not interact with any other feature, so its effect can be easily computed as \(e^{x_3}\).
import numpy as np
import matplotlib.pyplot as plt
import effector
np.random.seed(21)
bench = effector.benchmarks.ConditionalInteractionUniform()
model = bench.model
dataset = bench.dataset
x = bench.generate_data(1000)
PDP
Effector
Let's see below the PDP effects for each feature, using effector.
pdp = effector.PDP(x, model.predict, axis_limits=dataset.axis_limits)
pdp.fit(features="all", centering=True)
for feature in [0, 1, 2]:
pdp.plot(feature=feature, centering=True, y_limits=[-2, 2], heterogeneity=False)
PDP states that:
- \(x_1\) has a zero average effect on the model output
- \(x_2\) has a constant effect when \(x_2 <0\) or \(x_2 >0\), but when moving from \(x_2^-\) to \(x_2^+\) it adds (on average) \(+\frac{2}{3}\) units to \(y\)
- \(x_3\) has an effect of \(e^{x_3}\)
Derivations
How PDP leads to these explanations? Are they meaningfull? Let's have some analytical derivations. If you don't care about the derivations, skip the following three cells and go directly to the coclusions.
For \(x_1\):
For \(x_2\):
For \(x_3\):
Tests
# The closed form below lives in `effector.benchmarks` β the SAME function the
# test suite asserts against (tests/test_functional_conditional_interaction.py),
# so this notebook and the tests can never disagree about the right answer.
pdp_ground_truth = bench.pdp_gt
xx = np.linspace(-1, 1, 100)
y_pdp = []
for feature in [0, 1, 2]:
y_pdp.append(pdp_ground_truth(feature, xx))
plt.figure()
plt.title("PDP effects (ground truth)")
color_pallette = ["blue", "red", "green"]
feature_labels = ["Feature x1", "Feature x2", "Feature x3"]
for feature in [0, 1, 2]:
plt.plot(
xx,
y_pdp[feature],
color=color_pallette[feature],
linestyle="--",
label=feature_labels[feature],
)
plt.legend()
plt.xlim([-1.1, 1.1])
plt.ylim([-2, 2])
plt.xlabel("Feature Value")
plt.ylabel("PDP Effect")
plt.show()
# make a test
xx = np.linspace(-1, 1, 100)
for feature in [0, 1, 2]:
y_pdp = pdp.eval(feature=feature, xs=xx, centering=True)
y_gt = pdp_ground_truth(feature, xx)
np.testing.assert_allclose(y_pdp, y_gt, atol=1e-1)
Conclusions
Are the PDP effects intuitive?
- For \(x_1\) the effect is zero. The terms related to \(x_1\) are \(-x_1^2 \mathbb{1}_{x_2 <0}\) and \(x_1^2 \mathbb{1}_{x_2 \geq 0}\). Both terms involve an interaction with \(x_2\). Since \(x_2 \sim \mathcal{U}(-1,1)\), almost half of the instances have \(x_2^i < 0\) and the the other half \(x_2^i \geq 0\), so the the two terms cancel out.
- For \(x_2\), the effect is constant when \(x_2 < 0\) or \(x_2>0\) but has a positive jump of \(\frac{2}{3}\) when moving from \(x_2^-\) to \(x_2^+\). It makes sense; when \(x_2 < 0\) the active term is $-(x_1^i)^2 \mathbb{1}_{x_2 < 0} $ which adds a negative quantity to the output and when \(x_2 \geq 0\) the active term is \((x_1^i)^2 \mathbb{1}_{x_2 \geq 0}\) that adds something postive. Therefore in the transmission we observe a non-linearity.
- For \(x_3\), the effect is \(e^{x_3}\), as expected, since only the this term corresponds to \(x_3\) and has no interaction with other variables.
Feature importance and one-click explanation (new API)
Beyond plotting each effect, effector can summarize how much each feature
matters and produce a single auto-explanation:
fx.importances()returns a per-feature importance = the dispersion of the mean effect (the \(\mu\)-twin of heterogeneity). Here \(x_2\) (the jump) and \(x_3\) (the \(e^{x_3}\) curve) should rank above the flat-on-average \(x_1\).effector.explain(...)runs the whole pipeline and returns a serializableReportwith a self-contained HTML view.
# per-feature importance = dispersion of the mean effect (mu-twin of heterogeneity)
print("importances:", np.round(pdp.importances(), 3))
# one-click auto-explanation -> Report (serializable; self-contained HTML)
report = effector.explain(x, model.predict, method="pdp", nof_instances="all")
report.show()
# the whole notebook, in one page: the report published with this example
from pathlib import Path
_out = Path("reports") / "05_conditional_interaction_independent_uniform_global"
_out.mkdir(parents=True, exist_ok=True)
report.to_html(_out / "report_pdp.html")
importances: [0.007 0.326 0.668]
[effector] global effects (GAM) -> 86.7% of the model's variance
regional effects (CALM) -> 99.8%
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
PDP report Β· target: y
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
DATA & MODEL
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
instances 1,000
features 3 Β· 3 continuous
model output mean 1.17 Β· std 0.803 Β· range [-0.546, 3.38]
EXPLAINED VARIANCE
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
step split on solo ΞRΒ² RΒ² heter
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
GAM (all features global) β β 86.7% β
+ x_0 x_1 +13.1% +13.1% 99.8% 0.29 β 0.00
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
FINAL 99.8%
REJECTED SPLITS min gain 1.0%
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
feature split on solo ΞRΒ² reason
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
β x_1 x_0 +10.3% -10.2% redundant
β redundant: it would explain variance on its own (see solo),
but the accepted splits already account for it.
FEATURES ranked, in the selected snapshot
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
feature importance heter #regions
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
x_2 0.6679 ββββββββββββββββββ 0.0000 1
x_1 0.3261 βββββββββ 0.2919 1
x_0 0.2917 ββββββββ 0.0000 2
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
the features above carry 100% of the total importance mass
Feature 0 - Full partition tree:
π³ Full Tree Structure:
βββββββββββββββββββββββ
x_0 πΉ [id: 0 | heter: 0.29 | inst: 1000 | w: 1.00]
x_1 < 0.00 πΉ [id: 1 | heter: 0.00 | inst: 488 | w: 0.49]
x_1 β₯ 0.00 πΉ [id: 2 | heter: 0.00 | inst: 512 | w: 0.51]
--------------------------------------------------
Feature 0 - Statistics per tree level:
π³ Tree Summary:
βββββββββββββββββ
Level 0πΉheter: 0.29
Level 1πΉheter: 0.00 | π»0.29 (100.00%)
/home/givasile/github/packages/effector/effector/report.py:606: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
fig.tight_layout()
ALE
Effector
Let's see below the PDP effects for each feature, using effector.
ale = effector.ALE(x, model.predict, axis_limits=dataset.axis_limits)
ale.fit(features="all", centering=True, binning_method=effector.axis_partitioning.Fixed(nof_bins=31))
for feature in [0, 1, 2]:
ale.plot(feature=feature, centering=True, y_limits=[-2, 2], heterogeneity=False)
ALE states that: - \(x_1\) has a zero average effect on the model output (same as PDP) - \(x_2\) has a constant effect when \(x_2 <0\) or \(x_2 >0\), but when moving from \(x_2^-\) to \(x_2^+\) it adds (on average) \(+\frac{2}{3}\) units to \(y\) (same as PDP) - \(x_3\) has an effect of \(e^{x_3}\) (same as PDP)
Derivations
For all bins, except the central, it holds that bin limits are either both negative or both positive, so the effects cancel out. For central bin, i.e., the one from \(-\frac{2}{K}\) to \(\frac{2}{K}\), the effect is \(\frac{(2x^i_1)^2}{| \mathcal{S}_k |} = \frac{2}{3} \frac{K}{2} = \frac{K}{3}\). Therefore, the ALE effect is:
Tests
# The closed form below lives in `effector.benchmarks` β the SAME function the
# test suite asserts against (tests/test_functional_conditional_interaction.py),
# so this notebook and the tests can never disagree about the right answer.
ale_ground_truth = bench.ale_gt
xx = np.linspace(-1, 1, 100)
y_ale = []
for feature in [0, 1, 2]:
y_ale.append(ale_ground_truth(feature, xx))
plt.figure()
plt.title("ALE effects (ground truth)")
color_pallette = ["blue", "red", "green"]
feature_labels = ["Feature x1", "Feature x2", "Feature x3"]
for feature in [0, 1, 2]:
plt.plot(
xx,
y_ale[feature],
color=color_pallette[feature],
linestyle="--",
label=feature_labels[feature]
)
plt.legend()
plt.xlim([-1.1, 1.1])
plt.ylim([-2, 2])
plt.xlabel("Feature Value")
plt.ylabel("ALE Effect")
plt.show()
xx = np.linspace(-1, 1, 100)
for feature in [0, 1, 2]:
y_ale = ale.eval(feature=feature, xs=xx, centering=True)
y_gt = ale_ground_truth(feature, xx)
# hack to remove the effect at undefined region
if feature == 1:
K = 31
ind = np.logical_and(xx > -1/K, xx < 1/K)
y_ale[ind] = 0
y_gt[ind] = 0
np.testing.assert_allclose(y_ale, y_gt, atol=1e-1)
Conclusions
Are the ALE effects intuitive?
ALE effects are identical to PDP effects which, as discussed above, can be considered intutive.
RHALE
Effector
Let's see below the RHALE effects for each feature, using effector.
rhale = effector.RHALE(x, model.predict, model.jacobian, axis_limits=dataset.axis_limits)
rhale.fit(features="all", centering=True)
for feature in [0, 1, 2]:
rhale.plot(feature=feature, centering=True, y_limits=[-2, 2], heterogeneity=False)
RHALE states that: - \(x_1\) has a zero average effect on the model output (same as PDP) - \(x_2\) has a zero average effect on the model output (different than PDP and ALE) - \(x_3\) has an effect of \(e^{x_3}\) (same as PDP)
Derivations
Tests
# The closed form below lives in `effector.benchmarks` β the SAME function the
# test suite asserts against (tests/test_functional_conditional_interaction.py),
# so this notebook and the tests can never disagree about the right answer.
rhale_ground_truth = bench.rhale_gt
xx = np.linspace(-1, 1, 100)
y_rhale = []
for feature in [0, 1, 2]:
y_rhale.append(rhale_ground_truth(feature, xx))
plt.figure()
plt.title("RHALE effects (ground truth)")
color_pallette = ["blue", "red", "green"]
feature_labels = ["Feature x1", "Feature x2", "Feature x3"]
for feature in [0, 1, 2]:
plt.plot(
xx,
y_rhale[feature],
color=color_pallette[feature],
linestyle="-" if feature == 0 else "--",
label=feature_labels[feature]
)
plt.legend()
plt.xlim([-1.1, 1.1])
plt.ylim([-2, 2])
plt.xlabel("Feature Value")
plt.ylabel("RHALE Effect")
plt.show()
for feature in [0, 1, 2]:
y_ale = rhale.eval(feature=feature, xs=xx, centering=True)
y_gt = rhale_ground_truth(feature, xx)
np.testing.assert_allclose(y_ale, y_gt, atol=1e-1)
Conclusions
Are the RHALE effects intuitive?
RHALE does not add something new, compared to ALE and PDP, for features \(x_1\) and \(x_3\). For \(x_2\), however, it does not capture the abrupt increase by \(+\frac{2}{3}\) units when moving from \(x_2^-\) to \(x_2^+\), which can be considered as an error mode of RHALE. In fact, RHALE requires a differentiable black box model, and since \(f\) is not differentiable with respect to \(x_2\), that is why we get a slightly misleading result.
Cross-method sanity check
The one-liner effector.explain with every engine this notebook's model
supports. Everything must run end to end; the closing table puts the reads
side by side. Where methods disagree β ranking, accepted splits, RΒ² β that is
a property of the data/model worth a closer look, not an error.
from pathlib import Path
_out = Path("reports") / "05_conditional_interaction_independent_uniform_global"
_out.mkdir(parents=True, exist_ok=True)
# === cross-method sweep: effector.explain on every applicable engine ======
sweep_reports = {}
for _m in ["pdp", "derpdp", "ale", "rhale", "shapdp"]:
_kw = {"nof_instances": 300} if _m == "shapdp" else {}
print(f"--- {_m} " + "-" * 50)
sweep_reports[_m] = effector.explain(
x, model.predict, model.jacobian, method=_m, **_kw
)
if _m != "pdp": # the published report is the narrated one above
sweep_reports[_m].to_html(_out / f"report_{_m}.html")
print()
print(f"{'method':<8} {'ranking (plotted)':<44} {'GAM R2':>8} {'final R2':>9} splits")
for _m, _r in sweep_reports.items():
_rank = " > ".join(fr.name for fr in _r.features)
_ev = _r.explained_variance
if _ev:
_sp = "; ".join(f"{s['name']} on {s['on']}" for s in _ev["stages"]) or "none"
print(f"{_m:<8} {_rank:<44} {_ev['gam_r2']:>7.1%} {_ev['regional_r2']:>8.1%} {_sp}")
else:
print(f"{_m:<8} {_rank:<44} {'-':>7} {'-':>8} (derivative scale: no variance ledger)")
print(f"\nreports stored in {_out}/")
--- pdp --------------------------------------------------
[effector] global effects (GAM) -> 86.7% of the model's variance
regional effects (CALM) -> 99.8%
--- derpdp --------------------------------------------------
/home/givasile/github/packages/effector/effector/report.py:606: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
fig.tight_layout()
--- ale --------------------------------------------------
[effector] global effects (GAM) -> 85.3% of the model's variance
regional effects (CALM) -> 98.8%
/home/givasile/github/packages/effector/effector/report.py:606: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
fig.tight_layout()
--- rhale --------------------------------------------------
[effector] global effects (GAM) -> 70.0% of the model's variance
regional effects (CALM) -> 100.0%
/home/givasile/github/packages/effector/effector/report.py:606: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
fig.tight_layout()
--- shapdp --------------------------------------------------
[effector] global effects (GAM) -> 87.3% of the model's variance
regional effects (CALM) -> 96.6%
/home/givasile/github/packages/effector/effector/report.py:606: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
fig.tight_layout()
method ranking (plotted) GAM R2 final R2 splits
pdp x_2 > x_1 > x_0 86.7% 99.8% x_0 on x_1
derpdp x_2 - - (derivative scale: no variance ledger)
ale x_2 > x_1 > x_0 85.3% 98.8% x_0 on x_1
rhale x_2 > x_0 70.0% 100.0% x_0 on x_1
shapdp x_2 > x_1 > x_0 87.3% 96.6% x_0 on x_1
reports stored in reports/05_conditional_interaction_independent_uniform_global/











