More consistent handling of constraints in glex(): $shap and $remainder

[jemus42]

Constraining a decomposition via max_interaction or features drops terms. The
components no longer sum to the prediction, so SHAP values reconstructed from them
violate efficiency. Previously they were returned anyway.

Supersedes #18 and #25, closes #13 and #11.

$shap is NA when the decomposition is constrained

gl <- glex(xg, x, max_interaction = 2)
#> Warning: SHAP values set to NA: the decomposition is constrained by
#> `max_interaction` and does not sum to the full model prediction [...]

gl$constrained
#> [1] "max_interaction"

gl$shap
#> [1] NA                   # scalar: numeric consumers fail loudly

gl$m                      # unaffected

Previously the same call warns about nothing, has no $constrained, and returns
plausible-looking SHAP values (-0.201 -0.189 1.083 ...) that are silently wrong.

Complete decomposition:

gl <- glex(xg, x)
gl$constrained
#> character(0)            # length(gl$constrained) > 0 tells you if $shap is usable

A constraint that only drops zero terms is not a constraint:

gl <- glex(rp, x, max_interaction = 3)
#> The decomposition is constrained by `max_interaction`, but the dropped terms
#> are all zero: the components still sum to the model prediction, so SHAP values
#> remain valid.

gl$constrained
#> character(0)

New $remainder

What the dropped terms are collectively worth. Present exactly when $constrained is
non-empty, so the prediction reconstructs either way:

gl <- glex(xg, x, max_interaction = 2)
all.equal(gl$intercept + rowSums(gl$m) + gl$remainder,
          predict(xg, x, outputmargin = TRUE))
#> TRUE

glex(xg, x)$remainder
#> NULL                    # complete decomposition owes nothing

rpf objects already carried a $remainder from predict_components(); the other
learners had none. Now one field, one definition, computed in one place.

On the scale the model is decomposed on, the link scale for xgboost:

gl <- glex(xgb_binary, x, max_interaction = 1)          # binary:logistic
plogis(gl$intercept + rowSums(gl$m) + gl$remainder) == predict(xgb_binary, x)
#> TRUE

Unlike #25, this covers classification and other non-identity links. #25 computed
predict(object, x) - rowSums(m) — a response-scale prediction minus a margin-scale
decomposition — and had to restrict itself to regression.

glex() on rpf now also returns $shap (class-suffixed for multiclass, like $m),
and detects post-hoc constraints, which previously passed silently.

Bug fixes

glex.rpf() compared against the wrong scale and the wrong class. rpf decomposes the
raw score; predict() defaults to type = "prob", which applies rpf's response function
(clamp to [0, 1] for loss = "L2", inverse link for "logit" / "exponential").

# glex decomposes .pred_b; predict()[[1]] is .pred_a, and truncated
max(abs(intercept + rowSums(m) - predict(rp, d)[[1]]))                    #> 1.08
max(abs(intercept + rowSums(m) - predict(rp, d, type = "numeric")[[1]]))  #> 1.3e-15

The inertness check averaged real terms away. all.equal() reports mean relative
difference, so a discrepancy concentrated in a few observations passed as "all zero",
glex() kept SHAP values that were wrong. Caught by a macOS CI failure. Now judged
elementwise against $remainder, at the noise floor.

Tests

• test-constrained-remainder.R: the $remainder invariant across xgboost, ranger and
  rpf; both xgboost link scales (binary:logistic, count:poisson); all three rpf losses.
  Includes negative assertions that reconstruction must not match the response scale, so
  the target cannot be quietly reverted.
• test-rpf-sum-identity.R: was comparing rowSums(m) (raw score) against predict()'s
  clamped probabilities, which is why it needed tolerance = 0.03 / 0.8 and a
  min() over both class columns. Against type = "numeric" the identity is exact
  (~1e-15). Multiclass now asserts the residual is constant within a class — the class
  intercept rpf does not report — instead of merely bounding its magnitude.
• The inert-constraint test relied on a top-order term happening to come out zero for a
  given fit. It did on Linux, not on macOS (4e-4). Now uses a constant predictor, which
  cannot be split on, so its terms are zero by construction on every platform.
• Dropped a duplicate rpf test: weighting_method only applies where glex walks the trees
  itself (xgboost, ranger). glex.rpf() ignores it, so the "fastpd" and "path-dependent"
  variants were asserting the same thing twice.

[jemus42]

Turns out it's a bit more complex, because

a) Multiclass, once again, makes everything annoying
b) glex_explain previously didn't assume $shap exists and reconstructed it from the components (since I mostly worked with rpf for simplicity, where fitting models is easier because xgboost needs more preprocessing, and rpf doesn't store $shap)

... anyway, the glex.rpf method now adds a $shap field for consistency, which is NA if some constraint applies (which is now also tracked in a $constraint field).

[jemus42]

This now also goes a bit further and adds the $remainder term in a more principled way to also apply to the xgboost (and any other, actually), thereby superseding the part of #25 I was in favor of.