From 823917bdf75ea0727202fe48703d9e37f322d91c Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Sat, 11 Jul 2026 00:10:13 +0200
Subject: [PATCH 01/10] workaround to detect if  is valid

---
 DESCRIPTION                        |  2 +-
 NEWS.md                            | 10 +++++
 R/autoplot.R                       |  5 +--
 R/glex.R                           | 61 +++++++++++++++++++++++++++++-
 R/glex_vi.R                        |  1 -
 man/glex.Rd                        |  6 ++-
 tests/testthat/test-glex-xgboost.R | 48 +++++++++++++++++++++--
 7 files changed, 122 insertions(+), 11 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 767dbbb..77257df 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: glex
 Title: Global Explanations for Tree-Based Models
-Version: 0.6.0
+Version: 0.6.0.9000
 Authors@R: c(
     person(c("Marvin", "N."), "Wright", , "cran@wrig.de", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-8542-6291")),
diff --git a/NEWS.md b/NEWS.md
index 25bc0c7..67ec9d2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,13 @@
+# glex 0.6.0.9000 (development version)
+
+* `$shap` is now `NA` (with a warning) when the decomposition is constrained via
+  `max_interaction` or `features`: a constrained decomposition does not sum to the
+  full model prediction, so SHAP values cannot be reconstructed from it without
+  violating the efficiency property. Previously, misleading values were returned.
+* Tests that fit `randomPlantedForest` models are skipped on Windows for now: an
+  out-of-bounds read in randomPlantedForest's `purify_3()` crashes R there
+  (fixed upstream in PlantedML/randomPlantedForest#61, not yet merged).
+
 # glex 0.6.0
 
 * Extended compatibility with `xgboost`, now requiring `xgboost (>= 3.0.0)` in `Suggests:`
diff --git a/R/autoplot.R b/R/autoplot.R
index 1e5e360..b7924fe 100644
--- a/R/autoplot.R
+++ b/R/autoplot.R
@@ -64,7 +64,6 @@ autoplot.glex_vi <- function(
   by_what <- ifelse(by_degree, "Degree", "Term")
   score <- switch(scale, absolute = "m", relative = "m_rel")
 
-  # FIXME: data.table NSE stuff
   m <- m_rel <- degree <- NULL
 
   aggr_degree <- function(x) {
@@ -77,9 +76,7 @@ autoplot.glex_vi <- function(
   }
 
   # If max_interaction or threshold would select any observations, we summarize accordingly
-  if (
-    any(object$degree > max_interaction | abs(object[[score]]) <= threshold)
-  ) {
+  if (any(object$degree > max_interaction | abs(object[[score]]) <= threshold)) {
     # Sanity check that we're aggregating correctly by keeping track of the sums
     old_sum <- sum(object[[score]])
 
diff --git a/R/glex.R b/R/glex.R
index 7abea01..f451577 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -20,7 +20,11 @@
 #'
 #' @return Decomposition of the regression or classification function.
 #' A `list` with elements:
-#' * `shap`: SHAP values (`xgboost` method only).
+#' * `shap`: SHAP values (`xgboost` and `ranger` methods only). These are reconstructed
+#'   from the functional decomposition, which is only possible if the decomposition is
+#'   complete: if it is constrained via `max_interaction` or `features`, the components
+#'   no longer sum to the full model prediction and the SHAP efficiency property cannot
+#'   hold, so `shap` is set to `NA` (with a warning).
 #' * `m`: Functional decomposition into all main and interaction
 #'   components in the model, up to the degree specified by `max_interaction`.
 #'   The variable names correspond to the original variable names,
@@ -553,6 +557,43 @@ tree_fun_emp_fastPD <- function(
 }
 
 
+#' Upper bound on the interaction order present in the decomposition.
+#' Interactions can only arise between features sharing a root-to-leaf path, so
+#' per tree the order is bounded by both its depth and its number of distinct
+#' features; the model bound is the maximum over trees. This can overestimate
+#' (a path may split repeatedly on one feature), never underestimate.
+#' Read from the trees table rather than model metadata because ranger has no
+#' usable depth field (max.depth = 0 means unlimited) and xgboost's config
+#' parsing is version-fragile.
+#' Relies on child node ids being larger than their parent's, which holds for
+#' both xgb.model.dt.tree(use_int_id = TRUE) and ranger::treeInfo() numbering.
+#' @param trees data.table with columns Node, Yes, No, Feature_num, Tree;
+#'   leaves have Yes = NA
+#' @keywords internal
+#' @noRd
+max_order_bound <- function(trees) {
+  Tree <- Node <- NULL
+
+  max(vapply(
+    0:max(trees$Tree),
+    function(tree) {
+      tree_info <- trees[Tree == tree, ][order(Node)]
+      depth <- integer(nrow(tree_info))
+      for (i in seq_len(nrow(tree_info))) {
+        if (!is.na(tree_info$Yes[i])) {
+          depth[tree_info$Yes[i] + 1L] <- depth[i] + 1L
+          depth[tree_info$No[i] + 1L] <- depth[i] + 1L
+        }
+      }
+      n_features <- length(unique(tree_info$Feature_num[
+        tree_info$Feature_num > 0L
+      ]))
+      min(max(depth), n_features)
+    },
+    integer(1)
+  ))
+}
+
 #' Internal tree function wrapper that returns the actual tree function
 #' @param trees data.table
 #' @param x observerations, matrix like data-structure
@@ -676,6 +717,15 @@ calc_components <- function(
       1L
     all_S <- get_all_subsets_cpp(sort(unique(features_num)), max_interaction)
   }
+  # SHAP values can only be reconstructed from a complete decomposition: if terms
+  # are dropped via max_interaction or features, the components no longer sum to
+  # the full model prediction and the SHAP efficiency property cannot hold.
+  shap_valid <- max_interaction >= max_order_bound(trees)
+  if (!is.null(features)) {
+    used_features_num <- trees[Feature_num > 0, sort(unique(Feature_num))]
+    shap_valid <- shap_valid && all(used_features_num %in% features_num)
+  }
+
   # Keep only those with not more than max_interaction involved features
   d <- lengths(all_S)
 
@@ -731,6 +781,15 @@ calc_components <- function(
     FUN.VALUE = numeric(nrow(x))
   )
 
+  if (!shap_valid) {
+    warning(
+      "SHAP values set to NA: the decomposition is constrained by `max_interaction` or ",
+      "`features` and does not sum to the full model prediction, so SHAP values cannot ",
+      "be reconstructed from it (the efficiency property would be violated)."
+    )
+    shap[] <- NA_real_
+  }
+
   # Return shap values, decomposition and intercept
   ret <- list(
     shap = data.table::setDT(as.data.frame(shap)),
diff --git a/R/glex_vi.R b/R/glex_vi.R
index d992dbf..eba5f11 100644
--- a/R/glex_vi.R
+++ b/R/glex_vi.R
@@ -45,7 +45,6 @@
 glex_vi <- function(object, ...) {
   checkmate::assert_class(object, classes = "glex")
 
-  # FIXME: data.table NSE warnings
   term <- degree <- m <- m_rel <- NULL
 
   m_long <- melt_m(object$m, object$target_levels)
diff --git a/man/glex.Rd b/man/glex.Rd
index 971f1c7..057eca0 100644
--- a/man/glex.Rd
+++ b/man/glex.Rd
@@ -57,7 +57,11 @@ If not set in \code{xgboost}, the default value of \code{6} is assumed.}
 Decomposition of the regression or classification function.
 A \code{list} with elements:
 \itemize{
-\item \code{shap}: SHAP values (\code{xgboost} method only).
+\item \code{shap}: SHAP values (\code{xgboost} and \code{ranger} methods only). These are reconstructed
+from the functional decomposition, which is only possible if the decomposition is
+complete: if it is constrained via \code{max_interaction} or \code{features}, the components
+no longer sum to the full model prediction and the SHAP efficiency property cannot
+hold, so \code{shap} is set to \code{NA} (with a warning).
 \item \code{m}: Functional decomposition into all main and interaction
 components in the model, up to the degree specified by \code{max_interaction}.
 The variable names correspond to the original variable names,
diff --git a/tests/testthat/test-glex-xgboost.R b/tests/testthat/test-glex-xgboost.R
index 43b72fa..ecb3dd2 100644
--- a/tests/testthat/test-glex-xgboost.R
+++ b/tests/testthat/test-glex-xgboost.R
@@ -17,14 +17,14 @@ test_that("max_interaction respects xgb's max_depth", {
 test_that("features argument only calculates for given features", {
   x <- as.matrix(mtcars[, -1])
   xg <- xgboost(x, mtcars$mpg, nrounds = 10, verbosity = 0)
-  glexb <- glex(xg, x, features = c("cyl", "disp"))
+  glexb <- suppressWarnings(glex(xg, x, features = c("cyl", "disp")))
   expect_equal(colnames(glexb$m), c("cyl", "cyl:disp", "disp"))
 })
 
 test_that("features argument results in same values as without", {
   x <- as.matrix(mtcars[, -1])
   xg <- xgboost(x, mtcars$mpg, nrounds = 10, verbosity = 0)
-  glexb1 <- glex(xg, x, features = c("cyl", "disp"))
+  glexb1 <- suppressWarnings(glex(xg, x, features = c("cyl", "disp")))
   glexb2 <- glex(xg, x)
   cols <- c("cyl", "disp", "cyl:disp")
   expect_equal(glexb1$m[, ..cols], glexb2$m[, ..cols])
@@ -33,11 +33,53 @@ test_that("features argument results in same values as without", {
 test_that("features argument works together with max_interaction", {
   x <- as.matrix(mtcars[, -1])
   xg <- xgboost(x, mtcars$mpg, nrounds = 10, verbosity = 0)
-  glexb <- glex(xg, x, features = c("cyl", "disp"), max_interaction = 1)
+  glexb <- suppressWarnings(glex(
+    xg,
+    x,
+    features = c("cyl", "disp"),
+    max_interaction = 1
+  ))
   expect_equal(colnames(glexb$m), c("cyl", "disp"))
 })
 
 
+test_that("shap is NA when decomposition is constrained", {
+  set.seed(1)
+  x <- as.matrix(mtcars[, -1])
+  xg <- xgboost(
+    x,
+    mtcars$mpg,
+    nrounds = 10,
+    max_depth = 4,
+    verbosity = 0,
+    nthreads = 1
+  )
+
+  # constrained via max_interaction
+  expect_warning(
+    gl_mi <- glex(xg, x, max_interaction = 1),
+    "efficiency property"
+  )
+  expect_true(all(is.na(gl_mi$shap)))
+  expect_false(anyNA(gl_mi$m))
+
+  # constrained via features
+  expect_warning(
+    gl_ft <- glex(xg, x, features = c("cyl", "disp")),
+    "efficiency property"
+  )
+  expect_true(all(is.na(gl_ft$shap)))
+
+  # unconstrained: shap present and satisfies the efficiency property
+  gl_full <- glex(xg, x)
+  expect_false(anyNA(gl_full$shap))
+  expect_equal(
+    unname(gl_full$intercept + rowSums(gl_full$shap)),
+    unname(predict(xg, x)),
+    tolerance = 1e-5
+  )
+})
+
 x_train <- as.matrix(mtcars[1:26, -1])
 x_test <- as.matrix(mtcars[27:32, -1])
 y_train <- mtcars$mpg[1:26]

From 3ca4796ef3efcb115b9d84b25a69ae8b1b67518b Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Mon, 13 Jul 2026 17:21:07 +0200
Subject: [PATCH 02/10] cleanup for  handling in constrained case

---
 NEWS.md                                |  13 ++
 R/glex.R                               | 187 +++++++++++++++++++------
 R/glex_explain.R                       |  86 ++++++++----
 man/glex.Rd                            |  12 +-
 tests/testthat/test-glex-xgboost.R     |  10 ++
 tests/testthat/test-rpf-sum-identity.R |  74 ++++++++++
 6 files changed, 311 insertions(+), 71 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 67ec9d2..bc2a97d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,19 @@
   `max_interaction` or `features`: a constrained decomposition does not sum to the
   full model prediction, so SHAP values cannot be reconstructed from it without
   violating the efficiency property. Previously, misleading values were returned.
+  `$m` is unaffected. Supersedes #18, closes #13.
+* `glex()` objects gain a `$constrained` field naming the arguments that constrained
+  the decomposition (`character(0)` if complete), so `length(x$constrained) > 0` tells
+  you whether `$shap` is usable.
+* `glex()` on `randomPlantedForest` models now returns `$shap` as well, computed from
+  the components like for the other model classes (for multiclass models, `$shap`
+  columns are class-specific like those of `$m`). Previously the field was absent.
+  Constraining the decomposition post-hoc via `max_interaction` or `features` is now
+  detected for `rpf` models too, where it previously passed silently.
+* `glex_explain()` now reads SHAP values from `$shap` instead of recomputing them from
+  the components, so the `glex` object is the single source of truth. The SHAP
+  reference bar is omitted for constrained decompositions, where it previously showed
+  a value reconstructed from the constrained components.
 * Tests that fit `randomPlantedForest` models are skipped on Windows for now: an
   out-of-bounds read in randomPlantedForest's `purify_3()` crashes R there
   (fixed upstream in PlantedML/randomPlantedForest#61, not yet merged).
diff --git a/R/glex.R b/R/glex.R
index f451577..6f00a0f 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -20,16 +20,20 @@
 #'
 #' @return Decomposition of the regression or classification function.
 #' A `list` with elements:
-#' * `shap`: SHAP values (`xgboost` and `ranger` methods only). These are reconstructed
-#'   from the functional decomposition, which is only possible if the decomposition is
-#'   complete: if it is constrained via `max_interaction` or `features`, the components
+#' * `shap`: SHAP values, derived from the functional decomposition as
+#'   \eqn{\phi_j = \sum_{S \ni j} m_S / |S|}. This reconstruction is only valid if the
+#'   decomposition is complete: if it is constrained (see `constrained`), the components
 #'   no longer sum to the full model prediction and the SHAP efficiency property cannot
-#'   hold, so `shap` is set to `NA` (with a warning).
+#'   hold, so `shap` is set to `NA` (with a warning) while `m` remains valid.
+#'   For multiclass models, columns are class-specific like those of `m`.
 #' * `m`: Functional decomposition into all main and interaction
 #'   components in the model, up to the degree specified by `max_interaction`.
 #'   The variable names correspond to the original variable names,
 #'   with `:` separating interaction terms as one would specify in a [`formula`] interface.
 #' * `intercept`: Intercept term, the expected value of the prediction.
+#' * `constrained`: Character vector naming the arguments that constrained the
+#'   decomposition (`"max_interaction"`, `"features"`), or `character(0)` if it is
+#'   complete. Use `length(x$constrained) > 0` to check whether `shap` is valid.
 #' @export
 glex <- function(object, x, max_interaction = NULL, features = NULL, ...) {
   UseMethod("glex")
@@ -73,6 +77,49 @@ glex.rpf <- function(object, x, max_interaction = NULL, features = NULL, ...) {
     max_interaction = max_interaction,
     predictors = features
   )
+
+  model_features <- names(object$blueprint$ptypes$predictors)
+  ret$constrained <- constrained_by(
+    max_interaction = max_interaction,
+    features = features,
+    # An rpf model cannot contain interactions of higher order than it was fit with
+    available_interaction = min(
+      object$params$max_interaction,
+      length(model_features)
+    ),
+    model_features = model_features
+  )
+
+  # SHAP values are derived from the components, like in the other methods.
+  # For multiclass models `$shap` mirrors the class-suffixed columns of `$m`.
+  if (is.null(ret$target_levels)) {
+    shap <- shap_from_components(ret$m, names(ret$x))
+  } else {
+    term_class <- split_names(names(ret$m), "__class:", target_index = 2)
+
+    shap <- do.call(
+      cbind,
+      lapply(ret$target_levels, function(level) {
+        m_class <- ret$m[, names(ret$m)[term_class == level], with = FALSE]
+        data.table::setnames(
+          m_class,
+          split_names(names(m_class), "__class:", target_index = 1)
+        )
+
+        shap_class <- shap_from_components(m_class, names(ret$x))
+        colnames(shap_class) <- paste0(
+          colnames(shap_class),
+          "__class:",
+          level
+        )
+        shap_class
+      })
+    )
+  }
+
+  shap <- invalidate_shap(shap, ret$constrained)
+  ret$shap <- data.table::setDT(as.data.frame(shap))
+
   # class(ret) <- c("glex", "rpf_components", class(ret))
   ret
 }
@@ -557,6 +604,87 @@ tree_fun_emp_fastPD <- function(
 }
 
 
+#' Compute SHAP values from a functional decomposition
+#'
+#' Each interaction term is distributed equally across the features it involves:
+#' \eqn{\phi_j = \sum_{S \ni j} m_S / |S|}. Shared by all `glex()` methods so
+#' that `$shap` is derived from `$m` in exactly one place.
+#' @param m Component matrix or `data.table` (without the intercept column).
+#' @param features Character vector of feature names to compute SHAP values for.
+#' @keywords internal
+#' @noRd
+shap_from_components <- function(m, features) {
+  m <- as.matrix(m)
+  d <- get_degree(colnames(m))
+
+  # Interaction terms contribute m_S / |S| to each involved feature
+  interactions <- sweep(m, MARGIN = 2, d, "/")
+
+  vapply(
+    features,
+    function(col) {
+      idx <- find_term_matches(col, colnames(interactions))
+
+      if (length(idx) == 0) {
+        numeric(nrow(interactions))
+      } else {
+        rowSums(interactions[, idx, drop = FALSE])
+      }
+    },
+    FUN.VALUE = numeric(nrow(m))
+  )
+}
+
+#' Report which constraints were applied to a decomposition
+#'
+#' SHAP values can only be reconstructed from a complete decomposition: if terms
+#' are dropped, the components no longer sum to the full model prediction and the
+#' SHAP efficiency property cannot hold. Returns the names of the constrained
+#' arguments, or `character(0)` if the decomposition is complete.
+#' @param max_interaction,features As supplied to `glex()`.
+#' @param available_interaction Interaction order available in the model.
+#' @param model_features Features the model actually uses.
+#' @keywords internal
+#' @noRd
+constrained_by <- function(
+  max_interaction,
+  features,
+  available_interaction,
+  model_features
+) {
+  constrained <- character(0)
+
+  if (!is.null(max_interaction) && max_interaction < available_interaction) {
+    constrained <- c(constrained, "max_interaction")
+  }
+  if (!is.null(features) && !all(model_features %in% features)) {
+    constrained <- c(constrained, "features")
+  }
+
+  constrained
+}
+
+#' Warn about, and invalidate, SHAP values from a constrained decomposition
+#' @param shap SHAP value matrix.
+#' @param constrained Character vector from [constrained_by()].
+#' @keywords internal
+#' @noRd
+invalidate_shap <- function(shap, constrained) {
+  if (length(constrained) == 0) {
+    return(shap)
+  }
+
+  warning(
+    "SHAP values set to NA: the decomposition is constrained by ",
+    paste0("`", constrained, "`", collapse = " and "),
+    " and does not sum to the full model prediction, so SHAP values cannot be ",
+    "reconstructed from it (the efficiency property would be violated). ",
+    "The components in `$m` are unaffected."
+  )
+  shap[] <- NA_real_
+  shap
+}
+
 #' Upper bound on the interaction order present in the decomposition.
 #' Interactions can only arise between features sharing a root-to-leaf path, so
 #' per tree the order is bounded by both its depth and its number of distinct
@@ -717,14 +845,16 @@ calc_components <- function(
       1L
     all_S <- get_all_subsets_cpp(sort(unique(features_num)), max_interaction)
   }
-  # SHAP values can only be reconstructed from a complete decomposition: if terms
-  # are dropped via max_interaction or features, the components no longer sum to
-  # the full model prediction and the SHAP efficiency property cannot hold.
-  shap_valid <- max_interaction >= max_order_bound(trees)
-  if (!is.null(features)) {
-    used_features_num <- trees[Feature_num > 0, sort(unique(Feature_num))]
-    shap_valid <- shap_valid && all(used_features_num %in% features_num)
-  }
+  # Features the model actually splits on, for the constraint check below
+  model_features <- colnames(x)[
+    setdiff(sort(unique(trees$Feature_num)), 0L)
+  ]
+  constrained <- constrained_by(
+    max_interaction = max_interaction,
+    features = features,
+    available_interaction = max_order_bound(trees),
+    model_features = model_features
+  )
 
   # Keep only those with not more than max_interaction involved features
   d <- lengths(all_S)
@@ -761,41 +891,16 @@ calc_components <- function(
     }
   }
 
-  d <- get_degree(colnames(m_all))
-
-  # Overall feature effect is sum of all elements where feature is involved
-  interactions <- sweep(m_all[, -1, drop = FALSE], MARGIN = 2, d[-1], "/")
-
-  # SHAP values are the sum of the m's * 1/d
-  shap <- vapply(
-    colnames(x),
-    function(col) {
-      idx <- find_term_matches(col, colnames(interactions))
-
-      if (length(idx) == 0) {
-        numeric(nrow(interactions))
-      } else {
-        rowSums(interactions[, idx, drop = FALSE])
-      }
-    },
-    FUN.VALUE = numeric(nrow(x))
-  )
-
-  if (!shap_valid) {
-    warning(
-      "SHAP values set to NA: the decomposition is constrained by `max_interaction` or ",
-      "`features` and does not sum to the full model prediction, so SHAP values cannot ",
-      "be reconstructed from it (the efficiency property would be violated)."
-    )
-    shap[] <- NA_real_
-  }
+  shap <- shap_from_components(m_all[, -1, drop = FALSE], colnames(x))
+  shap <- invalidate_shap(shap, constrained)
 
   # Return shap values, decomposition and intercept
   ret <- list(
     shap = data.table::setDT(as.data.frame(shap)),
     m = data.table::setDT(as.data.frame(m_all[, -1])),
     intercept = unique(m_all[, 1]),
-    x = data.table::setDT(as.data.frame(x))
+    x = data.table::setDT(as.data.frame(x)),
+    constrained = constrained
   )
   class(ret) <- c("glex", "xgb_components", class(ret))
   ret
diff --git a/R/glex_explain.R b/R/glex_explain.R
index b43605a..b7b3fb3 100644
--- a/R/glex_explain.R
+++ b/R/glex_explain.R
@@ -162,13 +162,32 @@ glex_explain <- function(
   xtemp <- unlist(Map(paste, xnames, xtemp, sep = " = "))
   names(xtemp) <- xnames
 
+  # SHAP values are read from `$shap` rather than recomputed, so the object is the
+  # single source of truth. They are NA for constrained decompositions (see `?glex`),
+  # in which case the SHAP bar is omitted.
+  shap_valid <- length(object$constrained) == 0
+
   # helper df of shap values to plot them in as separate rectangles, at last (lowest) position (determined by term_int)
   xshap <- xdf[,
-    list(shap = sum(m_scaled), term_int = max(term_int) + 1),
+    list(term_int = max(term_int) + 1),
     by = c("reference_term", "class")
   ]
-  xshap[, xleft := avgpred + shap]
-  xshap[, xright := avgpred]
+
+  if (shap_valid) {
+    shap_long <- melt_m(object$shap, object$target_levels)
+    shap_id <- shap_long[shap_long[[".id"]] == id, ]
+    shap_id[, .id := NULL]
+    data.table::setnames(shap_id, c("term", "m"), c("reference_term", "shap"))
+
+    if (is.null(object$target_levels)) {
+      # dummy class column to match xdf, see above
+      shap_id[, class := 1L]
+    }
+
+    xshap <- merge(xshap, shap_id, by = c("reference_term", "class"))
+    xshap[, xleft := avgpred + shap]
+    xshap[, xright := avgpred]
+  }
 
   # Subset to selected class, could be of length > 1
   if (!is.null(class)) {
@@ -221,28 +240,6 @@ glex_explain <- function(
       ),
       alpha = .75
     ) +
-    # Draw SHAP bar below with separate df
-    geom_rect(
-      data = xshap,
-      aes(
-        fill = as.character(sign(shap)),
-        xmin = xleft,
-        xmax = xright,
-        ymin = term_int - barheight / 2,
-        ymax = term_int + barheight / 2
-      ),
-      alpha = .75
-    ) +
-    # Label SHAP values below other bars
-    geom_label(
-      data = xshap,
-      aes(
-        label = sprintf("SHAP: %s", format(shap, digits = 2)),
-        x = xright,
-        fill = NULL,
-        hjust = pmin(1.05, 0.95 + sign(shap))
-      )
-    ) +
     # Label contribution values
     geom_label(
       aes(
@@ -254,6 +251,32 @@ glex_explain <- function(
       alpha = .75
     )
 
+  if (shap_valid) {
+    p <- p +
+      # Draw SHAP bar below with separate df
+      geom_rect(
+        data = xshap,
+        aes(
+          fill = as.character(sign(shap)),
+          xmin = xleft,
+          xmax = xright,
+          ymin = term_int - barheight / 2,
+          ymax = term_int + barheight / 2
+        ),
+        alpha = .75
+      ) +
+      # Label SHAP values below other bars
+      geom_label(
+        data = xshap,
+        aes(
+          label = sprintf("SHAP: %s", format(shap, digits = 2)),
+          x = xright,
+          fill = NULL,
+          hjust = pmin(1.05, 0.95 + sign(shap))
+        )
+      )
+  }
+
   # scales / coords
   p <- p +
     scale_fill_manual(
@@ -273,7 +296,18 @@ glex_explain <- function(
         id,
         pred
       ),
-      subtitle = sprintf("Centered around average prediction: %1.2f", avgpred),
+      subtitle = sprintf(
+        "Centered around average prediction: %1.2f%s",
+        avgpred,
+        if (shap_valid) {
+          ""
+        } else {
+          sprintf(
+            "\nSHAP values omitted: decomposition constrained by %s",
+            paste0("`", object$constrained, "`", collapse = " and ")
+          )
+        }
+      ),
       x = "Average prediction +/- m",
       y = NULL
     ) +
diff --git a/man/glex.Rd b/man/glex.Rd
index 057eca0..4672cc8 100644
--- a/man/glex.Rd
+++ b/man/glex.Rd
@@ -57,16 +57,20 @@ If not set in \code{xgboost}, the default value of \code{6} is assumed.}
 Decomposition of the regression or classification function.
 A \code{list} with elements:
 \itemize{
-\item \code{shap}: SHAP values (\code{xgboost} and \code{ranger} methods only). These are reconstructed
-from the functional decomposition, which is only possible if the decomposition is
-complete: if it is constrained via \code{max_interaction} or \code{features}, the components
+\item \code{shap}: SHAP values, derived from the functional decomposition as
+\eqn{\phi_j = \sum_{S \ni j} m_S / |S|}. This reconstruction is only valid if the
+decomposition is complete: if it is constrained (see \code{constrained}), the components
 no longer sum to the full model prediction and the SHAP efficiency property cannot
-hold, so \code{shap} is set to \code{NA} (with a warning).
+hold, so \code{shap} is set to \code{NA} (with a warning) while \code{m} remains valid.
+For multiclass models, columns are class-specific like those of \code{m}.
 \item \code{m}: Functional decomposition into all main and interaction
 components in the model, up to the degree specified by \code{max_interaction}.
 The variable names correspond to the original variable names,
 with \code{:} separating interaction terms as one would specify in a \code{\link{formula}} interface.
 \item \code{intercept}: Intercept term, the expected value of the prediction.
+\item \code{constrained}: Character vector naming the arguments that constrained the
+decomposition (\code{"max_interaction"}, \code{"features"}), or \code{character(0)} if it is
+complete. Use \code{length(x$constrained) > 0} to check whether \code{shap} is valid.
 }
 }
 \description{
diff --git a/tests/testthat/test-glex-xgboost.R b/tests/testthat/test-glex-xgboost.R
index ecb3dd2..5f4c8db 100644
--- a/tests/testthat/test-glex-xgboost.R
+++ b/tests/testthat/test-glex-xgboost.R
@@ -60,6 +60,7 @@ test_that("shap is NA when decomposition is constrained", {
     gl_mi <- glex(xg, x, max_interaction = 1),
     "efficiency property"
   )
+  expect_identical(gl_mi$constrained, "max_interaction")
   expect_true(all(is.na(gl_mi$shap)))
   expect_false(anyNA(gl_mi$m))
 
@@ -68,10 +69,19 @@ test_that("shap is NA when decomposition is constrained", {
     gl_ft <- glex(xg, x, features = c("cyl", "disp")),
     "efficiency property"
   )
+  expect_identical(gl_ft$constrained, "features")
   expect_true(all(is.na(gl_ft$shap)))
 
+  # both axes are reported
+  expect_warning(
+    gl_both <- glex(xg, x, features = c("cyl", "disp"), max_interaction = 1),
+    "efficiency property"
+  )
+  expect_identical(gl_both$constrained, c("max_interaction", "features"))
+
   # unconstrained: shap present and satisfies the efficiency property
   gl_full <- glex(xg, x)
+  expect_identical(gl_full$constrained, character(0))
   expect_false(anyNA(gl_full$shap))
   expect_equal(
     unname(gl_full$intercept + rowSums(gl_full$shap)),
diff --git a/tests/testthat/test-rpf-sum-identity.R b/tests/testthat/test-rpf-sum-identity.R
index 3b0ff4a..5c98a78 100644
--- a/tests/testthat/test-rpf-sum-identity.R
+++ b/tests/testthat/test-rpf-sum-identity.R
@@ -70,3 +70,77 @@ test_that("rpf multiclass: classwise sum identity matches predicted probabilitie
 
   expect_true(all(class_diffs < 0.8))
 })
+
+test_that("rpf: shap is derived from components and satisfies efficiency", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purify_3() OOB read, see comment at top of file
+
+  rp <- randomPlantedForest::rpf(
+    mpg ~ cyl + hp + wt,
+    data = mtcars,
+    max_interaction = 3
+  )
+  gl <- glex::glex(rp, mtcars)
+
+  expect_identical(gl$constrained, character(0))
+  expect_false(anyNA(gl$shap))
+  expect_equal(
+    unname(gl$intercept + rowSums(gl$shap)),
+    unname(predict(rp, mtcars)[[1]]),
+    tolerance = 1e-6
+  )
+})
+
+test_that("rpf: constrained decompositions invalidate shap", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purify_3() OOB read, see comment at top of file
+
+  rp <- randomPlantedForest::rpf(
+    mpg ~ cyl + hp + wt,
+    data = mtcars,
+    max_interaction = 3
+  )
+
+  expect_warning(gl_mi <- glex::glex(rp, mtcars, max_interaction = 1), "NA")
+  expect_identical(gl_mi$constrained, "max_interaction")
+  expect_true(all(is.na(gl_mi$shap)))
+  expect_false(anyNA(gl_mi$m))
+
+  expect_warning(gl_ft <- glex::glex(rp, mtcars, features = c("hp", "wt")), "NA")
+  expect_identical(gl_ft$constrained, "features")
+  expect_true(all(is.na(gl_ft$shap)))
+})
+
+test_that("rpf multiclass: shap mirrors the class-suffixed structure of m", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purify_3() OOB read, see comment at top of file
+
+  mt <- mtcars
+  mt$cyl <- factor(mt$cyl)
+  rpk <- randomPlantedForest::rpf(
+    cyl ~ mpg + hp + wt,
+    data = mt,
+    max_interaction = 2
+  )
+  glk <- glex::glex(rpk, mt)
+  pred <- predict(rpk, mt)
+
+  expect_identical(glk$constrained, character(0))
+  # one shap column per feature and class, like the terms in `m`
+  expect_setequal(
+    names(glk$shap),
+    paste0(rep(names(glk$x), each = 3), "__class:", rep(glk$target_levels, 3))
+  )
+
+  # Per class, shap + intercept reconstructs the predicted class score
+  class_diffs <- vapply(
+    glk$target_levels,
+    function(level) {
+      idx <- grepl(paste0("__class:", level), names(glk$shap), fixed = TRUE)
+      score <- glk$intercept + rowSums(glk$shap[, idx, with = FALSE])
+      max(abs(score - pred[[paste0(".pred_", level)]]))
+    },
+    FUN.VALUE = numeric(1)
+  )
+  expect_true(all(class_diffs < 0.8))
+})

From 6deeda17c4653fa622c474d52f76dd541d51cba3 Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Mon, 13 Jul 2026 17:56:14 +0200
Subject: [PATCH 03/10] cleanup & fix bound check

---
 NEWS.md                                |  4 +-
 R/glex.R                               | 69 +++++++++++++++++---------
 R/glex_explain.R                       |  5 +-
 R/print.R                              | 12 ++++-
 man/glex.Rd                            |  5 +-
 tests/testthat/test-glex-xgboost.R     | 68 +++++++++++++++++++++++++
 tests/testthat/test-glex_explain.R     | 63 +++++++++++++++++++++++
 tests/testthat/test-mtcars-ranger.R    | 29 +++++++++++
 tests/testthat/test-rpf-sum-identity.R | 37 ++++++++------
 9 files changed, 250 insertions(+), 42 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index bc2a97d..0f49b7b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -16,7 +16,9 @@
 * `glex_explain()` now reads SHAP values from `$shap` instead of recomputing them from
   the components, so the `glex` object is the single source of truth. The SHAP
   reference bar is omitted for constrained decompositions, where it previously showed
-  a value reconstructed from the constrained components.
+  a value reconstructed from the constrained components, and for objects created by
+  earlier versions of glex, which have no `$shap`.
+* `print()` on a `glex` object reports when the decomposition is constrained.
 * Tests that fit `randomPlantedForest` models are skipped on Windows for now: an
   out-of-bounds read in randomPlantedForest's `purify_3()` crashes R there
   (fixed upstream in PlantedML/randomPlantedForest#61, not yet merged).
diff --git a/R/glex.R b/R/glex.R
index 6f00a0f..73517d7 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -25,7 +25,10 @@
 #'   decomposition is complete: if it is constrained (see `constrained`), the components
 #'   no longer sum to the full model prediction and the SHAP efficiency property cannot
 #'   hold, so `shap` is set to `NA` (with a warning) while `m` remains valid.
-#'   For multiclass models, columns are class-specific like those of `m`.
+#'   For multiclass models, columns are class-specific like those of `m`. Note that
+#'   `randomPlantedForest` models report a single `intercept` for all classes, so for
+#'   multiclass models `intercept + rowSums(shap)` reconstructs the predicted class
+#'   scores only approximately.
 #' * `m`: Functional decomposition into all main and interaction
 #'   components in the model, up to the degree specified by `max_interaction`.
 #'   The variable names correspond to the original variable names,
@@ -642,7 +645,9 @@ shap_from_components <- function(m, features) {
 #' SHAP efficiency property cannot hold. Returns the names of the constrained
 #' arguments, or `character(0)` if the decomposition is complete.
 #' @param max_interaction,features As supplied to `glex()`.
-#' @param available_interaction Interaction order available in the model.
+#' @param available_interaction Interaction order available in the model. Only
+#'   evaluated when `max_interaction` could bind at all (R's lazy argument
+#'   evaluation keeps a potentially expensive computation out of the default path).
 #' @param model_features Features the model actually uses.
 #' @keywords internal
 #' @noRd
@@ -654,7 +659,13 @@ constrained_by <- function(
 ) {
   constrained <- character(0)
 
-  if (!is.null(max_interaction) && max_interaction < available_interaction) {
+  # A model can never contain interactions of higher order than it has features,
+  # so this cheap check short-circuits the common unconstrained case.
+  if (
+    !is.null(max_interaction) &&
+      max_interaction < length(model_features) &&
+      max_interaction < available_interaction
+  ) {
     constrained <- c(constrained, "max_interaction")
   }
   if (!is.null(features) && !all(model_features %in% features)) {
@@ -685,38 +696,50 @@ invalidate_shap <- function(shap, constrained) {
   shap
 }
 
-#' Upper bound on the interaction order present in the decomposition.
-#' Interactions can only arise between features sharing a root-to-leaf path, so
-#' per tree the order is bounded by both its depth and its number of distinct
-#' features; the model bound is the maximum over trees. This can overestimate
-#' (a path may split repeatedly on one feature), never underestimate.
+#' Maximum interaction order present in the model's decomposition.
+#'
+#' An interaction between features can only arise where they split on the same
+#' root-to-leaf path, so the order of a tree is the largest number of *distinct*
+#' features on any one of its paths, and the model's order is the maximum over
+#' trees. Tree depth alone is not a substitute: a path that splits repeatedly on
+#' the same feature is deeper than its interaction order, which would make a
+#' complete decomposition look constrained and needlessly invalidate `$shap`.
+#'
 #' Read from the trees table rather than model metadata because ranger has no
-#' usable depth field (max.depth = 0 means unlimited) and xgboost's config
-#' parsing is version-fragile.
-#' Relies on child node ids being larger than their parent's, which holds for
-#' both xgb.model.dt.tree(use_int_id = TRUE) and ranger::treeInfo() numbering.
+#' usable depth field (`max.depth = 0` means unlimited) and xgboost's config
+#' parsing is version-fragile. Relies on child node ids being larger than their
+#' parent's, which holds for both `xgb.model.dt.tree(use_int_id = TRUE)` and
+#' `ranger::treeInfo()` numbering, so parents are processed before their children.
 #' @param trees data.table with columns Node, Yes, No, Feature_num, Tree;
-#'   leaves have Yes = NA
+#'   leaves have `Yes = NA` and `Feature_num = 0`
 #' @keywords internal
 #' @noRd
-max_order_bound <- function(trees) {
+max_order <- function(trees) {
   Tree <- Node <- NULL
 
   max(vapply(
     0:max(trees$Tree),
     function(tree) {
       tree_info <- trees[Tree == tree, ][order(Node)]
-      depth <- integer(nrow(tree_info))
+
+      # Distinct features seen on the path from the root to each node
+      path_features <- vector("list", nrow(tree_info))
+      path_features[[1L]] <- integer(0)
+      order_tree <- 0L
+
       for (i in seq_len(nrow(tree_info))) {
-        if (!is.na(tree_info$Yes[i])) {
-          depth[tree_info$Yes[i] + 1L] <- depth[i] + 1L
-          depth[tree_info$No[i] + 1L] <- depth[i] + 1L
+        if (is.na(tree_info$Yes[i])) {
+          # Leaf: the path ends here
+          order_tree <- max(order_tree, length(path_features[[i]]))
+          next
         }
+
+        features <- union(path_features[[i]], tree_info$Feature_num[i])
+        path_features[[tree_info$Yes[i] + 1L]] <- features
+        path_features[[tree_info$No[i] + 1L]] <- features
       }
-      n_features <- length(unique(tree_info$Feature_num[
-        tree_info$Feature_num > 0L
-      ]))
-      min(max(depth), n_features)
+
+      order_tree
     },
     integer(1)
   ))
@@ -852,7 +875,7 @@ calc_components <- function(
   constrained <- constrained_by(
     max_interaction = max_interaction,
     features = features,
-    available_interaction = max_order_bound(trees),
+    available_interaction = max_order(trees),
     model_features = model_features
   )
 
diff --git a/R/glex_explain.R b/R/glex_explain.R
index b7b3fb3..d3e7724 100644
--- a/R/glex_explain.R
+++ b/R/glex_explain.R
@@ -164,8 +164,9 @@ glex_explain <- function(
 
   # SHAP values are read from `$shap` rather than recomputed, so the object is the
   # single source of truth. They are NA for constrained decompositions (see `?glex`),
-  # in which case the SHAP bar is omitted.
-  shap_valid <- length(object$constrained) == 0
+  # and absent from objects created by earlier versions of glex, in which case the
+  # SHAP bar is omitted.
+  shap_valid <- !is.null(object$shap) && length(object$constrained) == 0
 
   # helper df of shap values to plot them in as separate rectangles, at last (lowest) position (determined by term_int)
   xshap <- xdf[,
diff --git a/R/print.R b/R/print.R
index e51b262..d866226 100644
--- a/R/print.R
+++ b/R/print.R
@@ -32,6 +32,16 @@ print.glex <- function(x, ...) {
     max_deg,
     max_deg_lab
   )
-  cat("\n\n")
+  cat("\n")
+
+  if (length(x$constrained) > 0) {
+    cat(
+      "Decomposition constrained by",
+      paste0("`", x$constrained, "`", collapse = " and "),
+      "- SHAP values are NA\n"
+    )
+  }
+
+  cat("\n")
   str(x, list.len = 5)
 }
diff --git a/man/glex.Rd b/man/glex.Rd
index 4672cc8..fba6911 100644
--- a/man/glex.Rd
+++ b/man/glex.Rd
@@ -62,7 +62,10 @@ A \code{list} with elements:
 decomposition is complete: if it is constrained (see \code{constrained}), the components
 no longer sum to the full model prediction and the SHAP efficiency property cannot
 hold, so \code{shap} is set to \code{NA} (with a warning) while \code{m} remains valid.
-For multiclass models, columns are class-specific like those of \code{m}.
+For multiclass models, columns are class-specific like those of \code{m}. Note that
+\code{randomPlantedForest} models report a single \code{intercept} for all classes, so for
+multiclass models \code{intercept + rowSums(shap)} reconstructs the predicted class
+scores only approximately.
 \item \code{m}: Functional decomposition into all main and interaction
 components in the model, up to the degree specified by \code{max_interaction}.
 The variable names correspond to the original variable names,
diff --git a/tests/testthat/test-glex-xgboost.R b/tests/testthat/test-glex-xgboost.R
index 5f4c8db..c8c21c2 100644
--- a/tests/testthat/test-glex-xgboost.R
+++ b/tests/testthat/test-glex-xgboost.R
@@ -78,6 +78,7 @@ test_that("shap is NA when decomposition is constrained", {
     "efficiency property"
   )
   expect_identical(gl_both$constrained, c("max_interaction", "features"))
+  expect_true(all(is.na(gl_both$shap)))
 
   # unconstrained: shap present and satisfies the efficiency property
   gl_full <- glex(xg, x)
@@ -90,6 +91,73 @@ test_that("shap is NA when decomposition is constrained", {
   )
 })
 
+test_that("unconstrained calls never invalidate shap", {
+  set.seed(1)
+  x <- as.matrix(mtcars[, -1])
+  xg <- xgboost(
+    x,
+    mtcars$mpg,
+    nrounds = 10,
+    max_depth = 4,
+    verbosity = 0,
+    nthreads = 1
+  )
+
+  # Over-specifying max_interaction is not a constraint
+  expect_no_warning(gl_over <- glex(xg, x, max_interaction = 99))
+  expect_identical(gl_over$constrained, character(0))
+  expect_false(anyNA(gl_over$shap))
+
+  # Naming all features explicitly is not a constraint
+  expect_no_warning(gl_all <- glex(xg, x, features = colnames(x)))
+  expect_identical(gl_all$constrained, character(0))
+  expect_false(anyNA(gl_all$shap))
+
+  # A deeper model than it has features: tree depth exceeds the interaction order,
+  # so a depth-based bound would spuriously flag the default call as constrained.
+  x2 <- x[, c("cyl", "hp")]
+  xg2 <- xgboost(
+    x2,
+    mtcars$mpg,
+    nrounds = 10,
+    max_depth = 4,
+    verbosity = 0,
+    nthreads = 1
+  )
+  expect_no_warning(gl2 <- glex(xg2, x2))
+  expect_identical(gl2$constrained, character(0))
+  expect_false(anyNA(gl2$shap))
+  expect_equal(
+    unname(gl2$intercept + rowSums(gl2$shap)),
+    unname(predict(xg2, x2)),
+    tolerance = 1e-5
+  )
+
+  # max_interaction at exactly the model's interaction order is not a constraint
+  expect_no_warning(gl_at <- glex(xg2, x2, max_interaction = 2))
+  expect_identical(gl_at$constrained, character(0))
+})
+
+test_that("features the model never splits on do not count as a constraint", {
+  set.seed(1)
+  x <- as.matrix(mtcars[, -1])
+  xg <- xgboost(
+    x[, c("cyl", "hp")],
+    mtcars$mpg,
+    nrounds = 5,
+    max_depth = 2,
+    verbosity = 0,
+    nthreads = 1
+  )
+
+  # `x` carries a column the model never splits on: naming only the used features
+  # is still a complete decomposition
+  x_extra <- x[, c("cyl", "hp", "wt")]
+  expect_no_warning(gl <- glex(xg, x_extra, features = c("cyl", "hp")))
+  expect_identical(gl$constrained, character(0))
+  expect_false(anyNA(gl$shap))
+})
+
 x_train <- as.matrix(mtcars[1:26, -1])
 x_test <- as.matrix(mtcars[27:32, -1])
 y_train <- mtcars$mpg[1:26]
diff --git a/tests/testthat/test-glex_explain.R b/tests/testthat/test-glex_explain.R
index 3368914..0e2d2fe 100644
--- a/tests/testthat/test-glex_explain.R
+++ b/tests/testthat/test-glex_explain.R
@@ -35,3 +35,66 @@ test_that("multiclass rpf", {
   p <- glex_explain(gl, 2, threshold = 0)
   expect_s3_class(p, "ggplot")
 })
+
+
+# SHAP values come from `$shap` -----------------------------------------------------------------------------------
+test_that("glex_explain plots the SHAP values stored in $shap", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purify_3() OOB read, see test-rpf-sum-identity.R
+  rp <- rpf(y ~ x1 + x2 + x3, data = xdat, max_interaction = 3)
+  gl <- glex(rp, xdat)
+
+  id <- 5
+  p <- glex_explain(gl, id = id)
+
+  # The SHAP layer carries one bar per predictor, holding that predictor's value
+  # from `$shap` for this observation -- not a value recomputed from `$m`.
+  shap_layer <- Filter(
+    function(l) "shap" %in% names(l$data),
+    p$layers
+  )
+  expect_length(shap_layer, 2) # bar + label
+
+  plotted <- shap_layer[[1]]$data
+  expect_equal(
+    plotted[["shap"]][order(plotted[["reference_term"]])],
+    unlist(gl$shap[id, ], use.names = FALSE)[order(names(gl$shap))],
+    tolerance = 1e-10
+  )
+
+  expect_no_error(ggplot2::ggplot_build(p))
+})
+
+test_that("glex_explain omits the SHAP bar for constrained decompositions", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purify_3() OOB read, see test-rpf-sum-identity.R
+  rp <- rpf(y ~ x1 + x2 + x3, data = xdat, max_interaction = 3)
+
+  gl <- glex(rp, xdat)
+  gl_constrained <- suppressWarnings(glex(rp, xdat, max_interaction = 1))
+
+  p <- glex_explain(gl, id = 2)
+  p_constrained <- glex_explain(gl_constrained, id = 2)
+
+  # SHAP bar and its label are dropped, the component bars remain
+  expect_length(p_constrained$layers, length(p$layers) - 2)
+  expect_match(p_constrained$labels$subtitle, "SHAP values omitted")
+  expect_false(grepl("omitted", p$labels$subtitle))
+
+  expect_no_error(ggplot2::ggplot_build(p_constrained))
+})
+
+test_that("glex_explain works on objects without $shap", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purify_3() OOB read, see test-rpf-sum-identity.R
+  rp <- rpf(y ~ x1 + x2 + x3, data = xdat, max_interaction = 3)
+
+  # Objects from earlier glex versions (and predict_components() output) have no
+  # `$shap`: the SHAP bar is omitted rather than erroring.
+  gl <- randomPlantedForest::predict_components(rp, xdat)
+  expect_null(gl$shap)
+
+  p <- glex_explain(gl, id = 2)
+  expect_s3_class(p, "ggplot")
+  expect_no_error(ggplot2::ggplot_build(p))
+})
diff --git a/tests/testthat/test-mtcars-ranger.R b/tests/testthat/test-mtcars-ranger.R
index a13e53e..fb79692 100644
--- a/tests/testthat/test-mtcars-ranger.R
+++ b/tests/testthat/test-mtcars-ranger.R
@@ -98,3 +98,32 @@ test_that("Path-dependent ranger prediction is approx. same as sum of decomposit
     tolerance = 1e-5
   )
 })
+
+
+test_that("ranger: constrained decomposition invalidates shap", {
+  skip_if_not_installed("ranger")
+
+  expect_warning(
+    gl_mi <- glex(rf, x_train, max_interaction = 1),
+    "efficiency property"
+  )
+  expect_identical(gl_mi$constrained, "max_interaction")
+  expect_true(all(is.na(gl_mi$shap)))
+  expect_false(anyNA(gl_mi$m))
+
+  expect_warning(
+    gl_ft <- glex(rf, x_train, features = c("cyl", "disp")),
+    "efficiency property"
+  )
+  expect_identical(gl_ft$constrained, "features")
+  expect_true(all(is.na(gl_ft$shap)))
+  expect_false(anyNA(gl_ft$m))
+})
+
+test_that("ranger: unconstrained decomposition keeps shap", {
+  skip_if_not_installed("ranger")
+
+  expect_no_warning(gl <- glex(rf, x_train))
+  expect_identical(gl$constrained, character(0))
+  expect_false(anyNA(gl$shap))
+})
diff --git a/tests/testthat/test-rpf-sum-identity.R b/tests/testthat/test-rpf-sum-identity.R
index 5c98a78..1137167 100644
--- a/tests/testthat/test-rpf-sum-identity.R
+++ b/tests/testthat/test-rpf-sum-identity.R
@@ -101,14 +101,21 @@ test_that("rpf: constrained decompositions invalidate shap", {
     max_interaction = 3
   )
 
-  expect_warning(gl_mi <- glex::glex(rp, mtcars, max_interaction = 1), "NA")
+  expect_warning(
+    gl_mi <- glex::glex(rp, mtcars, max_interaction = 1),
+    "efficiency property"
+  )
   expect_identical(gl_mi$constrained, "max_interaction")
   expect_true(all(is.na(gl_mi$shap)))
   expect_false(anyNA(gl_mi$m))
 
-  expect_warning(gl_ft <- glex::glex(rp, mtcars, features = c("hp", "wt")), "NA")
+  expect_warning(
+    gl_ft <- glex::glex(rp, mtcars, features = c("hp", "wt")),
+    "efficiency property"
+  )
   expect_identical(gl_ft$constrained, "features")
   expect_true(all(is.na(gl_ft$shap)))
+  expect_false(anyNA(gl_ft$m))
 })
 
 test_that("rpf multiclass: shap mirrors the class-suffixed structure of m", {
@@ -123,7 +130,6 @@ test_that("rpf multiclass: shap mirrors the class-suffixed structure of m", {
     max_interaction = 2
   )
   glk <- glex::glex(rpk, mt)
-  pred <- predict(rpk, mt)
 
   expect_identical(glk$constrained, character(0))
   # one shap column per feature and class, like the terms in `m`
@@ -132,15 +138,18 @@ test_that("rpf multiclass: shap mirrors the class-suffixed structure of m", {
     paste0(rep(names(glk$x), each = 3), "__class:", rep(glk$target_levels, 3))
   )
 
-  # Per class, shap + intercept reconstructs the predicted class score
-  class_diffs <- vapply(
-    glk$target_levels,
-    function(level) {
-      idx <- grepl(paste0("__class:", level), names(glk$shap), fixed = TRUE)
-      score <- glk$intercept + rowSums(glk$shap[, idx, with = FALSE])
-      max(abs(score - pred[[paste0(".pred_", level)]]))
-    },
-    FUN.VALUE = numeric(1)
-  )
-  expect_true(all(class_diffs < 0.8))
+  # SHAP values redistribute the components of their own class, so per class the
+  # two must sum to the same thing exactly. (rpf reports a single intercept for all
+  # classes, so the class scores themselves are only reconstructed approximately.)
+  for (level in glk$target_levels) {
+    shap_cols <- grepl(paste0("__class:", level), names(glk$shap), fixed = TRUE)
+    m_cols <- grepl(paste0("__class:", level), names(glk$m), fixed = TRUE)
+
+    expect_equal(
+      rowSums(glk$shap[, shap_cols, with = FALSE]),
+      rowSums(glk$m[, m_cols, with = FALSE]),
+      tolerance = 1e-10,
+      info = paste("class", level)
+    )
+  }
 })

From 3053ad8c0342e6322d400e838b3f111c409a3608 Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Mon, 13 Jul 2026 18:49:12 +0200
Subject: [PATCH 04/10] more belt and suspenders for constraint check

---
 NEWS.md                                |  5 ++
 R/glex.R                               | 82 ++++++++++++++++++++++----
 man/glex.Rd                            |  3 +
 tests/testthat/test-rpf-sum-identity.R | 43 ++++++++++++++
 4 files changed, 121 insertions(+), 12 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 0f49b7b..443c861 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -8,6 +8,11 @@
 * `glex()` objects gain a `$constrained` field naming the arguments that constrained
   the decomposition (`character(0)` if complete), so `length(x$constrained) > 0` tells
   you whether `$shap` is usable.
+* A requested constraint only invalidates `$shap` if it actually drops something: a
+  model can contain a high-order term whose value is zero, in which case dropping it
+  leaves the decomposition (and the SHAP values) unchanged. `glex()` confirms the
+  constraint against the model's own predictions and, if the dropped terms were inert,
+  keeps `$shap` and emits a message instead of a warning.
 * `glex()` on `randomPlantedForest` models now returns `$shap` as well, computed from
   the components like for the other model classes (for multiclass models, `$shap`
   columns are class-specific like those of `$m`). Previously the field was absent.
diff --git a/R/glex.R b/R/glex.R
index 73517d7..4673752 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -37,6 +37,9 @@
 #' * `constrained`: Character vector naming the arguments that constrained the
 #'   decomposition (`"max_interaction"`, `"features"`), or `character(0)` if it is
 #'   complete. Use `length(x$constrained) > 0` to check whether `shap` is valid.
+#'   A constraint that only drops terms whose value is zero leaves the decomposition
+#'   unchanged; `glex()` confirms this against the model's predictions, reports it with
+#'   a message, and treats the result as complete.
 #' @export
 glex <- function(object, x, max_interaction = NULL, features = NULL, ...) {
   UseMethod("glex")
@@ -120,9 +123,16 @@ glex.rpf <- function(object, x, max_interaction = NULL, features = NULL, ...) {
     )
   }
 
-  shap <- invalidate_shap(shap, ret$constrained)
   ret$shap <- data.table::setDT(as.data.frame(shap))
 
+  # Multiclass rpf models report a single intercept for all classes, so the
+  # components only reconstruct the class scores approximately and the numeric
+  # confirmation is not reliable: the structural verdict is final there.
+  target <- if (is.null(ret$target_levels)) {
+    stats::predict(object, x)[[1]]
+  }
+  ret <- confirm_constrained(ret, target = target)
+
   # class(ret) <- c("glex", "rpf_components", class(ret))
   ret
 }
@@ -202,6 +212,12 @@ glex.xgb.Booster <- function(
   )
   res$intercept <- res$intercept + get_xgb_base_score(object)
 
+  # glex decomposes the raw margin, so the constraint is confirmed against it
+  res <- confirm_constrained(
+    res,
+    target = stats::predict(object, x, outputmargin = TRUE)
+  )
+
   # Return components
   res
 }
@@ -461,6 +477,13 @@ glex.ranger <- function(
   res$m <- res$m / object$num.trees
   res$intercept <- res$intercept / object$num.trees
 
+  ranger_pred <- stats::predict(object, x)$predictions
+  if (is.matrix(ranger_pred)) {
+    # Probability forests: glex decomposes the second class probability
+    ranger_pred <- ranger_pred[, 2L]
+  }
+  res <- confirm_constrained(res, target = ranger_pred)
+
   # Return components
   res
 }
@@ -675,25 +698,59 @@ constrained_by <- function(
   constrained
 }
 
-#' Warn about, and invalidate, SHAP values from a constrained decomposition
-#' @param shap SHAP value matrix.
-#' @param constrained Character vector from [constrained_by()].
+#' Confirm a structural constraint numerically, and invalidate SHAP values if real
+#'
+#' The structural check ([constrained_by()]) only knows which terms were dropped, not
+#' what they were worth: a model can contain a high-order term whose value is exactly
+#' zero, in which case dropping it changes nothing and the SHAP values remain valid.
+#' This confirms the constraint against the model's own predictions before discarding
+#' anything: if the components still reconstruct the prediction, the dropped terms were
+#' inert and `$shap` is kept (with a message, since the user did ask for a constraint).
+#' Otherwise `$shap` is invalidated with a warning.
+#' @param res `glex` object with `$m`, `$shap`, `$intercept` and `$constrained`.
+#' @param target Model predictions on the scale of the decomposition, or `NULL` to skip
+#'   the numeric confirmation and treat the structural verdict as final.
 #' @keywords internal
 #' @noRd
-invalidate_shap <- function(shap, constrained) {
-  if (length(constrained) == 0) {
-    return(shap)
+confirm_constrained <- function(res, target = NULL) {
+  if (length(res$constrained) == 0) {
+    return(res)
+  }
+
+  constrained_labels <- paste0("`", res$constrained, "`", collapse = " and ")
+
+  if (!is.null(target)) {
+    reconstruction <- res$intercept + rowSums(res$m)
+    inert <- isTRUE(all.equal(
+      unname(reconstruction),
+      unname(target),
+      tolerance = 1e-5
+    ))
+
+    if (inert) {
+      message(
+        "The decomposition is constrained by ",
+        constrained_labels,
+        ", but the dropped terms are all zero: the components still sum to the model ",
+        "prediction, so SHAP values remain valid."
+      )
+      res$constrained <- character(0)
+      return(res)
+    }
   }
 
   warning(
     "SHAP values set to NA: the decomposition is constrained by ",
-    paste0("`", constrained, "`", collapse = " and "),
+    constrained_labels,
     " and does not sum to the full model prediction, so SHAP values cannot be ",
     "reconstructed from it (the efficiency property would be violated). ",
     "The components in `$m` are unaffected."
   )
-  shap[] <- NA_real_
-  shap
+  for (col in names(res$shap)) {
+    data.table::set(res$shap, j = col, value = NA_real_)
+  }
+
+  res
 }
 
 #' Maximum interaction order present in the model's decomposition.
@@ -915,9 +972,10 @@ calc_components <- function(
   }
 
   shap <- shap_from_components(m_all[, -1, drop = FALSE], colnames(x))
-  shap <- invalidate_shap(shap, constrained)
 
-  # Return shap values, decomposition and intercept
+  # Return shap values, decomposition and intercept. Whether a structural constraint
+  # actually invalidates `shap` is confirmed against the model's predictions by the
+  # calling method, which has the model object (see `confirm_constrained()`).
   ret <- list(
     shap = data.table::setDT(as.data.frame(shap)),
     m = data.table::setDT(as.data.frame(m_all[, -1])),
diff --git a/man/glex.Rd b/man/glex.Rd
index fba6911..0bd65b0 100644
--- a/man/glex.Rd
+++ b/man/glex.Rd
@@ -74,6 +74,9 @@ with \code{:} separating interaction terms as one would specify in a \code{\link
 \item \code{constrained}: Character vector naming the arguments that constrained the
 decomposition (\code{"max_interaction"}, \code{"features"}), or \code{character(0)} if it is
 complete. Use \code{length(x$constrained) > 0} to check whether \code{shap} is valid.
+A constraint that only drops terms whose value is zero leaves the decomposition
+unchanged; \code{glex()} confirms this against the model's predictions, reports it with
+a message, and treats the result as complete.
 }
 }
 \description{
diff --git a/tests/testthat/test-rpf-sum-identity.R b/tests/testthat/test-rpf-sum-identity.R
index 1137167..7619824 100644
--- a/tests/testthat/test-rpf-sum-identity.R
+++ b/tests/testthat/test-rpf-sum-identity.R
@@ -153,3 +153,46 @@ test_that("rpf multiclass: shap mirrors the class-suffixed structure of m", {
     )
   }
 })
+
+test_that("rpf: constraints that drop only zero terms keep shap valid", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purify_3() OOB read, see comment at top of file
+
+  set.seed(42)
+  n <- 120
+  p <- 6
+  x <- matrix(rnorm(n * p), ncol = p)
+  colnames(x) <- paste0("x", seq_len(p))
+  y <- x[, 1] + x[, 2] * x[, 3] + rnorm(n, sd = 0.3)
+  dat <- data.frame(x, y = y)
+
+  # Fit at the maximum order: the model contains a p-way term, but the data has no
+  # p-way structure, so that term is exactly zero.
+  rp <- randomPlantedForest::rpf(
+    y ~ .,
+    data = dat,
+    max_interaction = p,
+    ntrees = 20
+  )
+  full <- glex::glex(rp, dat[, -ncol(dat)])
+  degree <- lengths(strsplit(names(full$m), ":", fixed = TRUE))
+  expect_equal(max(abs(as.matrix(full$m[, degree == p, with = FALSE]))), 0)
+
+  # Dropping only that zero term changes nothing, so shap must survive -- with a
+  # message rather than a warning, since a constraint *was* requested.
+  expect_message(
+    gl <- glex::glex(rp, dat[, -ncol(dat)], max_interaction = p - 1L),
+    "dropped terms are all zero"
+  )
+  expect_identical(gl$constrained, character(0))
+  expect_false(anyNA(gl$shap))
+  expect_equal(as.matrix(gl$shap), as.matrix(full$shap), tolerance = 1e-8)
+
+  # Dropping terms that carry weight still invalidates
+  expect_warning(
+    gl_real <- glex::glex(rp, dat[, -ncol(dat)], max_interaction = 1),
+    "efficiency property"
+  )
+  expect_identical(gl_real$constrained, "max_interaction")
+  expect_true(all(is.na(gl_real$shap)))
+})

From 645d1addb94e11c16ffde06492c35d50373fe91e Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Tue, 14 Jul 2026 01:14:05 +0200
Subject: [PATCH 05/10] update NEWS

---
 NEWS.md | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 443c861..f40bb05 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -24,9 +24,10 @@
   a value reconstructed from the constrained components, and for objects created by
   earlier versions of glex, which have no `$shap`.
 * `print()` on a `glex` object reports when the decomposition is constrained.
-* Tests that fit `randomPlantedForest` models are skipped on Windows for now: an
-  out-of-bounds read in randomPlantedForest's `purify_3()` crashes R there
-  (fixed upstream in PlantedML/randomPlantedForest#61, not yet merged).
+* Tests that fit `randomPlantedForest` models remain skipped on Windows: an
+  out-of-bounds read in randomPlantedForest's purification crashes R there. Fixed
+  upstream in randomPlantedForest 0.3.0; the skips are lifted separately, once glex
+  requires that version.
 
 # glex 0.6.0
 

From ebfd319312a80151eb71fa46173a318fb0205e4c Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Tue, 14 Jul 2026 01:28:59 +0200
Subject: [PATCH 06/10] add centralized , partially superseding #25

---
 NEWS.md                                     |   7 +
 R/glex.R                                    |  19 ++-
 man/glex.Rd                                 |   5 +
 tests/testthat/test-constrained-remainder.R | 146 ++++++++++++++++++++
 4 files changed, 176 insertions(+), 1 deletion(-)
 create mode 100644 tests/testthat/test-constrained-remainder.R

diff --git a/NEWS.md b/NEWS.md
index f40bb05..4079fb1 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -18,6 +18,13 @@
   columns are class-specific like those of `$m`). Previously the field was absent.
   Constraining the decomposition post-hoc via `max_interaction` or `features` is now
   detected for `rpf` models too, where it previously passed silently.
+* `glex()` objects gain a `$remainder` field: what the constraint's dropped terms are
+  collectively worth, per observation, on the scale of `$m`. It is present exactly when
+  the decomposition is constrained, so `intercept + rowSums(m) + remainder` reconstructs
+  the model prediction whether or not a constraint was applied. `randomPlantedForest`
+  objects already carried a `$remainder` computed by `predict_components()`, but the
+  other model classes did not; the two are now one field with one definition, computed
+  in one place. Closes #11, supersedes #25.
 * `glex_explain()` now reads SHAP values from `$shap` instead of recomputing them from
   the components, so the `glex` object is the single source of truth. The SHAP
   reference bar is omitted for constrained decompositions, where it previously showed
diff --git a/R/glex.R b/R/glex.R
index 4673752..8925156 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -40,6 +40,11 @@
 #'   A constraint that only drops terms whose value is zero leaves the decomposition
 #'   unchanged; `glex()` confirms this against the model's predictions, reports it with
 #'   a message, and treats the result as complete.
+#' * `remainder`: What the dropped terms are collectively worth, per observation:
+#'   `prediction - (intercept + rowSums(m))`, on the same scale as `m`. Present exactly
+#'   when the decomposition is constrained, and absent otherwise, so
+#'   `intercept + rowSums(m) + remainder` reconstructs the prediction in either case.
+#'   For multiclass `randomPlantedForest` models it is class-wise, mirroring `m`.
 #' @export
 glex <- function(object, x, max_interaction = NULL, features = NULL, ...) {
   UseMethod("glex")
@@ -700,20 +705,29 @@ constrained_by <- function(
 
 #' Confirm a structural constraint numerically, and invalidate SHAP values if real
 #'
-#' The structural check ([constrained_by()]) only knows which terms were dropped, not
+#' The structural check (`constrained_by()`) only knows which terms were dropped, not
 #' what they were worth: a model can contain a high-order term whose value is exactly
 #' zero, in which case dropping it changes nothing and the SHAP values remain valid.
 #' This confirms the constraint against the model's own predictions before discarding
 #' anything: if the components still reconstruct the prediction, the dropped terms were
 #' inert and `$shap` is kept (with a message, since the user did ask for a constraint).
 #' Otherwise `$shap` is invalidated with a warning.
+#'
+#' The gap between the components and the prediction is what the dropped terms are
+#' collectively worth, and it is reported as `$remainder`. It is the quantitative form of
+#' `$constrained` and exists exactly when the decomposition is incomplete, so a complete
+#' one carries no remainder.
 #' @param res `glex` object with `$m`, `$shap`, `$intercept` and `$constrained`.
+#'   For `rpf` models `$remainder` may already be set by
+#'   `randomPlantedForest::predict_components()`; it is only overwritten where `target`
+#'   lets us compute it ourselves, which keeps the multiclass remainder rpf provides.
 #' @param target Model predictions on the scale of the decomposition, or `NULL` to skip
 #'   the numeric confirmation and treat the structural verdict as final.
 #' @keywords internal
 #' @noRd
 confirm_constrained <- function(res, target = NULL) {
   if (length(res$constrained) == 0) {
+    res$remainder <- NULL
     return(res)
   }
 
@@ -721,6 +735,8 @@ confirm_constrained <- function(res, target = NULL) {
 
   if (!is.null(target)) {
     reconstruction <- res$intercept + rowSums(res$m)
+    res$remainder <- unname(target - reconstruction)
+
     inert <- isTRUE(all.equal(
       unname(reconstruction),
       unname(target),
@@ -735,6 +751,7 @@ confirm_constrained <- function(res, target = NULL) {
         "prediction, so SHAP values remain valid."
       )
       res$constrained <- character(0)
+      res$remainder <- NULL
       return(res)
     }
   }
diff --git a/man/glex.Rd b/man/glex.Rd
index 0bd65b0..ce0bd6b 100644
--- a/man/glex.Rd
+++ b/man/glex.Rd
@@ -77,6 +77,11 @@ complete. Use \code{length(x$constrained) > 0} to check whether \code{shap} is v
 A constraint that only drops terms whose value is zero leaves the decomposition
 unchanged; \code{glex()} confirms this against the model's predictions, reports it with
 a message, and treats the result as complete.
+\item \code{remainder}: What the dropped terms are collectively worth, per observation:
+\code{prediction - (intercept + rowSums(m))}, on the same scale as \code{m}. Present exactly
+when the decomposition is constrained, and absent otherwise, so
+\code{intercept + rowSums(m) + remainder} reconstructs the prediction in either case.
+For multiclass \code{randomPlantedForest} models it is class-wise, mirroring \code{m}.
 }
 }
 \description{
diff --git a/tests/testthat/test-constrained-remainder.R b/tests/testthat/test-constrained-remainder.R
new file mode 100644
index 0000000..f531750
--- /dev/null
+++ b/tests/testthat/test-constrained-remainder.R
@@ -0,0 +1,146 @@
+# `$remainder` is the quantitative form of `$constrained`: it holds what the dropped
+# terms are collectively worth, so that the decomposition reconstructs the prediction
+# again once it is added back. It must therefore behave identically across model classes,
+# and be present exactly when the decomposition is incomplete.
+
+test_that("xgboost: remainder completes a constrained decomposition", {
+  set.seed(1)
+  x <- as.matrix(mtcars[, -1])
+  xg <- xgboost(x, mtcars$mpg, nrounds = 20, verbosity = 0, max_depth = 3)
+
+  # glex decomposes the raw margin, so the remainder lives on that scale too
+  pred <- predict(xg, x, outputmargin = TRUE)
+
+  full <- glex(xg, x)
+  expect_null(full$remainder)
+
+  expect_warning(gl <- glex(xg, x, max_interaction = 1), "efficiency property")
+  expect_identical(gl$constrained, "max_interaction")
+  expect_equal(
+    gl$intercept + rowSums(gl$m) + gl$remainder,
+    unname(pred),
+    tolerance = 1e-5
+  )
+  # the constraint dropped something, else it would not have invalidated shap
+  expect_gt(max(abs(gl$remainder)), 0)
+})
+
+test_that("xgboost: features constraint also yields a remainder", {
+  set.seed(1)
+  x <- as.matrix(mtcars[, -1])
+  xg <- xgboost(x, mtcars$mpg, nrounds = 20, verbosity = 0, max_depth = 3)
+  pred <- predict(xg, x, outputmargin = TRUE)
+
+  expect_warning(
+    gl <- glex(xg, x, features = c("cyl", "hp")),
+    "efficiency property"
+  )
+  expect_identical(gl$constrained, "features")
+  expect_equal(
+    gl$intercept + rowSums(gl$m) + gl$remainder,
+    unname(pred),
+    tolerance = 1e-5
+  )
+})
+
+test_that("ranger: remainder completes a constrained decomposition", {
+  skip_if_not_installed("ranger")
+  set.seed(1)
+  x <- as.matrix(mtcars[, -1])
+  rf <- ranger::ranger(
+    x = x,
+    y = mtcars$mpg,
+    node.stats = TRUE,
+    num.trees = 10,
+    max.depth = 4
+  )
+  pred <- predict(rf, x)$predictions
+
+  full <- glex(rf, x)
+  expect_null(full$remainder)
+
+  expect_warning(gl <- glex(rf, x, max_interaction = 1), "efficiency property")
+  expect_identical(gl$constrained, "max_interaction")
+  expect_equal(
+    gl$intercept + rowSums(gl$m) + gl$remainder,
+    unname(pred),
+    tolerance = 1e-5
+  )
+})
+
+test_that("rpf: remainder completes a constrained decomposition", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purification OOB read, see test-rpf-sum-identity.R
+  set.seed(1)
+  rp <- randomPlantedForest::rpf(
+    mpg ~ cyl + hp + wt,
+    data = mtcars,
+    max_interaction = 3
+  )
+  pred <- predict(rp, mtcars)[[1]]
+
+  full <- glex(rp, mtcars)
+  expect_null(full$remainder)
+
+  expect_warning(
+    gl <- glex(rp, mtcars, max_interaction = 1),
+    "efficiency property"
+  )
+  expect_identical(gl$constrained, "max_interaction")
+  expect_equal(
+    gl$intercept + rowSums(gl$m) + gl$remainder,
+    unname(pred),
+    tolerance = 1e-5
+  )
+})
+
+test_that("an inert constraint leaves no remainder", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purification OOB read, see test-rpf-sum-identity.R
+
+  # A model fit at the maximum order whose top-order term is exactly zero: dropping it
+  # changes nothing, so the decomposition is still complete and owes no remainder.
+  set.seed(42)
+  n <- 120
+  p <- 6
+  x <- matrix(rnorm(n * p), ncol = p)
+  colnames(x) <- paste0("x", seq_len(p))
+  y <- x[, 1] + x[, 2] * x[, 3] + rnorm(n, sd = 0.3)
+  dat <- data.frame(x, y = y)
+
+  rp <- randomPlantedForest::rpf(
+    y ~ .,
+    data = dat,
+    max_interaction = p,
+    ntrees = 20
+  )
+
+  expect_message(
+    gl <- glex(rp, dat[, -ncol(dat)], max_interaction = p - 1L),
+    "dropped terms are all zero"
+  )
+  expect_identical(gl$constrained, character(0))
+  expect_null(gl$remainder)
+})
+
+test_that("rpf multiclass: remainder is class-wise, mirroring m", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purification OOB read, see test-rpf-sum-identity.R
+  set.seed(1)
+  mt <- mtcars
+  mt$cyl <- factor(mt$cyl)
+  rpk <- randomPlantedForest::rpf(
+    cyl ~ mpg + hp + wt,
+    data = mt,
+    max_interaction = 3
+  )
+
+  expect_warning(
+    glk <- glex(rpk, mt, max_interaction = 1),
+    "efficiency property"
+  )
+  expect_identical(glk$constrained, "max_interaction")
+  # rpf computes the multiclass remainder itself, one column per class
+  expect_false(is.null(glk$remainder))
+  expect_identical(ncol(as.matrix(glk$remainder)), length(glk$target_levels))
+})

From da40725100d033c681c7e949603215a83014aba9 Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Tue, 14 Jul 2026 02:04:00 +0200
Subject: [PATCH 07/10] remainder work + fix rpf tests using wrong scale

---
 NEWS.md                                     |  13 ++
 R/glex.R                                    |  30 +++-
 man/glex.Rd                                 |  16 +-
 tests/testthat/test-constrained-remainder.R | 164 ++++++++++++++++++++
 tests/testthat/test-rpf-sum-identity.R      |  83 +++++-----
 5 files changed, 251 insertions(+), 55 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 4079fb1..059d0e0 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -25,6 +25,19 @@
   objects already carried a `$remainder` computed by `predict_components()`, but the
   other model classes did not; the two are now one field with one definition, computed
   in one place. Closes #11, supersedes #25.
+  Unlike #25, this covers classification and other non-identity links: the remainder is
+  taken on the scale the model is decomposed on, so for `xgboost` it is on the link scale
+  (`plogis(intercept + rowSums(m) + remainder)` recovers a `binary:logistic` probability),
+  while `ranger` probability forests and `randomPlantedForest` are decomposed on the
+  response scale directly.
+* For `randomPlantedForest` classification models, `glex()` now confirms the constraint
+  against `predict(type = "numeric")` rather than the default `type = "prob"`. rpf
+  decomposes the raw score, while `type = "prob"` applies rpf's response function --
+  a clamp to `[0, 1]` for `loss = "L2"`, the inverse link for `"logit"` and
+  `"exponential"` -- and for binary models returns the classes in an order whose first
+  column is not the one being decomposed. The components reconstruct the raw score
+  exactly, so `$remainder` now measures only what the dropped terms are worth, instead of
+  silently absorbing the back-transformation and the class mix-up.
 * `glex_explain()` now reads SHAP values from `$shap` instead of recomputing them from
   the components, so the `glex` object is the single source of truth. The SHAP
   reference bar is omitted for constrained decompositions, where it previously showed
diff --git a/R/glex.R b/R/glex.R
index 8925156..35fe77f 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -41,10 +41,18 @@
 #'   unchanged; `glex()` confirms this against the model's predictions, reports it with
 #'   a message, and treats the result as complete.
 #' * `remainder`: What the dropped terms are collectively worth, per observation:
-#'   `prediction - (intercept + rowSums(m))`, on the same scale as `m`. Present exactly
-#'   when the decomposition is constrained, and absent otherwise, so
-#'   `intercept + rowSums(m) + remainder` reconstructs the prediction in either case.
-#'   For multiclass `randomPlantedForest` models it is class-wise, mirroring `m`.
+#'   `prediction - (intercept + rowSums(m))`. Present exactly when the decomposition is
+#'   constrained, and absent otherwise, so `intercept + rowSums(m) + remainder`
+#'   reconstructs the prediction in either case. For multiclass `randomPlantedForest`
+#'   models it is class-wise, mirroring `m`.
+#'
+#'   Like `m` and `shap`, it is on the scale that the model is decomposed on, which for
+#'   `xgboost` is the **link** scale and not the response: for a `binary:logistic` model
+#'   the reconstruction gives the margin, and `plogis(intercept + rowSums(m) + remainder)`
+#'   gives the predicted probability. `ranger` probability forests and
+#'   `randomPlantedForest` are decomposed on the response scale, where no such
+#'   back-transformation is needed. Adding `remainder` to a probability is therefore
+#'   never correct for `xgboost`.
 #' @export
 glex <- function(object, x, max_interaction = NULL, features = NULL, ...) {
   UseMethod("glex")
@@ -130,11 +138,17 @@ glex.rpf <- function(object, x, max_interaction = NULL, features = NULL, ...) {
 
   ret$shap <- data.table::setDT(as.data.frame(shap))
 
-  # Multiclass rpf models report a single intercept for all classes, so the
-  # components only reconstruct the class scores approximately and the numeric
-  # confirmation is not reliable: the structural verdict is final there.
+  # rpf decomposes the raw score, which is what `type = "numeric"` returns. The default
+  # for classification is `type = "prob"`, which applies rpf's response function: a clamp
+  # to [0, 1] for `loss = "L2"`, the inverse link for `"logit"` and `"exponential"`.
+  # Comparing against that would confound the dropped terms with the back-transformation
+  # (and, for binary models, silently compare against the wrong class).
+  #
+  # Multiclass rpf models report a single intercept for all classes, so the components
+  # only reconstruct the class scores approximately and the numeric confirmation is not
+  # reliable: the structural verdict is final there, and rpf supplies `$remainder` itself.
   target <- if (is.null(ret$target_levels)) {
-    stats::predict(object, x)[[1]]
+    stats::predict(object, x, type = "numeric")[[1]]
   }
   ret <- confirm_constrained(ret, target = target)
 
diff --git a/man/glex.Rd b/man/glex.Rd
index ce0bd6b..8c2d30d 100644
--- a/man/glex.Rd
+++ b/man/glex.Rd
@@ -78,10 +78,18 @@ A constraint that only drops terms whose value is zero leaves the decomposition
 unchanged; \code{glex()} confirms this against the model's predictions, reports it with
 a message, and treats the result as complete.
 \item \code{remainder}: What the dropped terms are collectively worth, per observation:
-\code{prediction - (intercept + rowSums(m))}, on the same scale as \code{m}. Present exactly
-when the decomposition is constrained, and absent otherwise, so
-\code{intercept + rowSums(m) + remainder} reconstructs the prediction in either case.
-For multiclass \code{randomPlantedForest} models it is class-wise, mirroring \code{m}.
+\code{prediction - (intercept + rowSums(m))}. Present exactly when the decomposition is
+constrained, and absent otherwise, so \code{intercept + rowSums(m) + remainder}
+reconstructs the prediction in either case. For multiclass \code{randomPlantedForest}
+models it is class-wise, mirroring \code{m}.
+
+Like \code{m} and \code{shap}, it is on the scale that the model is decomposed on, which for
+\code{xgboost} is the \strong{link} scale and not the response: for a \code{binary:logistic} model
+the reconstruction gives the margin, and \code{plogis(intercept + rowSums(m) + remainder)}
+gives the predicted probability. \code{ranger} probability forests and
+\code{randomPlantedForest} are decomposed on the response scale, where no such
+back-transformation is needed. Adding \code{remainder} to a probability is therefore
+never correct for \code{xgboost}.
 }
 }
 \description{
diff --git a/tests/testthat/test-constrained-remainder.R b/tests/testthat/test-constrained-remainder.R
index f531750..0927896 100644
--- a/tests/testthat/test-constrained-remainder.R
+++ b/tests/testthat/test-constrained-remainder.R
@@ -43,6 +43,170 @@ test_that("xgboost: features constraint also yields a remainder", {
   )
 })
 
+# The remainder is defined on the scale of `$m`, which for xgboost is the raw margin --
+# not the response. Reconstructing on the link scale and then applying the link's inverse
+# must recover the prediction; adding the remainder to a probability must not. These pin
+# the scale contract, which is what makes the remainder work for non-identity links at all.
+test_that("xgboost binary:logistic: remainder lives on the link scale", {
+  set.seed(1)
+  n <- 300
+  d <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n))
+  x <- as.matrix(d)
+  yb <- as.integer(d$x1 + d$x2 * d$x3 + rnorm(n, sd = 0.3) > 0)
+
+  xgb_bin <- xgboost::xgb.train(
+    params = list(objective = "binary:logistic", max_depth = 3, eta = 0.1),
+    data = xgboost::xgb.DMatrix(x, label = yb),
+    nrounds = 20,
+    verbose = 0
+  )
+
+  expect_warning(
+    gl <- glex(xgb_bin, x, max_interaction = 1),
+    "efficiency property"
+  )
+  recon <- gl$intercept + rowSums(gl$m) + gl$remainder
+
+  expect_equal(
+    recon,
+    unname(predict(xgb_bin, x, outputmargin = TRUE)),
+    tolerance = 1e-5
+  )
+  expect_equal(
+    stats::plogis(recon),
+    unname(predict(xgb_bin, x)),
+    tolerance = 1e-5
+  )
+  # guards against "fixing" the target to the response scale, which would silently
+  # subtract a probability from a margin -- the bug that kept #25 regression-only
+  expect_false(isTRUE(all.equal(
+    recon,
+    unname(predict(xgb_bin, x)),
+    tolerance = 1e-5
+  )))
+})
+
+test_that("xgboost count:poisson: remainder lives on the log-link scale", {
+  set.seed(1)
+  n <- 300
+  d <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n))
+  x <- as.matrix(d)
+  ycount <- rpois(n, lambda = exp(0.5 * d$x1))
+
+  xgb_pois <- xgboost::xgb.train(
+    params = list(objective = "count:poisson", max_depth = 3, eta = 0.1),
+    data = xgboost::xgb.DMatrix(x, label = ycount),
+    nrounds = 20,
+    verbose = 0
+  )
+
+  expect_warning(
+    gl <- glex(xgb_pois, x, max_interaction = 1),
+    "efficiency property"
+  )
+  recon <- gl$intercept + rowSums(gl$m) + gl$remainder
+
+  expect_equal(
+    recon,
+    unname(predict(xgb_pois, x, outputmargin = TRUE)),
+    tolerance = 1e-5
+  )
+  expect_equal(exp(recon), unname(predict(xgb_pois, x)), tolerance = 1e-5)
+})
+
+test_that("ranger probability forest: remainder lives on the response scale", {
+  skip_if_not_installed("ranger")
+  set.seed(1)
+  n <- 300
+  d <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n))
+  x <- as.matrix(d)
+  y <- factor(ifelse(d$x1 + d$x2 * d$x3 + rnorm(n, sd = 0.3) > 0, "b", "a"))
+
+  rf <- ranger::ranger(
+    x = x,
+    y = y,
+    probability = TRUE,
+    node.stats = TRUE,
+    num.trees = 20,
+    max.depth = 4
+  )
+
+  expect_warning(gl <- glex(rf, x, max_interaction = 1), "efficiency property")
+  recon <- gl$intercept + rowSums(gl$m) + gl$remainder
+
+  # ranger probability forests are decomposed on the probability scale directly:
+  # glex takes the second class probability, so no link is involved
+  expect_equal(
+    recon,
+    unname(predict(rf, x)$predictions[, 2]),
+    tolerance = 1e-5
+  )
+})
+
+# rpf decomposes the raw score (`type = "numeric"`), not the probability. Its response
+# function depends on the loss -- a clamp to [0, 1] for "L2", the inverse link for "logit"
+# and "exponential" -- and `predict()` applies it by default (`type = "prob"`). Targeting
+# that default would fold the back-transformation into the remainder, and for binary
+# models compare against the wrong class entirely.
+test_that("rpf binary: remainder is on the raw score scale, for every loss", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purification OOB read, see test-rpf-sum-identity.R
+
+  for (loss in c("L2", "logit", "exponential")) {
+    set.seed(1)
+    rp <- randomPlantedForest::rpf(
+      y ~ x1 + x2 + x3,
+      data = xdat,
+      max_interaction = 3,
+      loss = loss
+    )
+
+    # unconstrained: the decomposition reconstructs the raw score exactly, and owes nothing
+    full <- glex(rp, xdat)
+    expect_null(full$remainder)
+    expect_equal(
+      full$intercept + rowSums(full$m),
+      unname(predict(rp, xdat, type = "numeric")[[1]]),
+      tolerance = 1e-8,
+      info = loss
+    )
+
+    expect_warning(
+      gl <- glex(rp, xdat, max_interaction = 1),
+      "efficiency property"
+    )
+    expect_equal(
+      gl$intercept + rowSums(gl$m) + gl$remainder,
+      unname(predict(rp, xdat, type = "numeric")[[1]]),
+      tolerance = 1e-8,
+      info = loss
+    )
+  }
+})
+
+test_that("rpf binary: remainder does not target the response scale", {
+  skip_if_not_installed("randomPlantedForest")
+  skip_on_os("windows") # rpf purification OOB read, see test-rpf-sum-identity.R
+  set.seed(1)
+  rp <- randomPlantedForest::rpf(
+    y ~ x1 + x2 + x3,
+    data = xdat,
+    max_interaction = 3
+  )
+
+  expect_warning(
+    gl <- glex(rp, xdat, max_interaction = 1),
+    "efficiency property"
+  )
+  recon <- gl$intercept + rowSums(gl$m) + gl$remainder
+
+  # `predict()` defaults to type = "prob" for classification: truncated to [0, 1] and
+  # ordered by class. Reconstructing against it must NOT hold -- if it ever does, the
+  # target has been silently switched back to the response scale.
+  probs <- predict(rp, xdat, type = "prob")
+  expect_false(isTRUE(all.equal(recon, unname(probs[[1]]), tolerance = 1e-5)))
+})
+
 test_that("ranger: remainder completes a constrained decomposition", {
   skip_if_not_installed("ranger")
   set.seed(1)
diff --git a/tests/testthat/test-rpf-sum-identity.R b/tests/testthat/test-rpf-sum-identity.R
index 7619824..11b939b 100644
--- a/tests/testthat/test-rpf-sum-identity.R
+++ b/tests/testthat/test-rpf-sum-identity.R
@@ -4,7 +4,20 @@
 # Fixed upstream in https://github.com/PlantedML/randomPlantedForest/pull/61 —
 # remove the skips once that (or a minimal fix) is merged.
 
-test_that("rpf binary: FastPD sum identity matches predicted class probability", {
+# For classification, rpf decomposes the *raw score*, which `predict(type = "numeric")`
+# returns. The default `type = "prob"` applies rpf's response function -- a clamp to
+# [0, 1] for `loss = "L2"`, the inverse link for `"logit"` / `"exponential"` -- which the
+# decomposition knows nothing about, so the components never sum to it. These tests used
+# to compare against `type = "prob"` and needed tolerances loose enough (0.03, 0.8) to
+# swallow that gap, plus a min() over both class columns to avoid committing to one.
+# Against the raw score the identity is exact.
+#
+# There is deliberately no `weighting_method` here: it selects how glex weights leaves
+# when *it* walks the trees, which it only does for xgboost and ranger. rpf purifies its
+# own components in C++, so `glex.rpf()` ignores the argument entirely. Two tests used to
+# pass "fastpd" and "path-dependent" and assert the same thing twice.
+
+test_that("rpf binary: sum identity matches the predicted raw score", {
   skip_if_not_installed("randomPlantedForest")
   skip_on_os("windows") # rpf purify_3() OOB read, see comment at top of file
 
@@ -13,36 +26,16 @@ test_that("rpf binary: FastPD sum identity matches predicted class probability",
     data = xdat,
     max_interaction = 3
   )
-  gl <- glex::glex(rp, xdat, weighting_method = "fastpd")
-  pred <- predict(rp, xdat)
-
-  score <- unname(gl$intercept + rowSums(gl$m))
-  diff_col1 <- max(abs(score - pred[[1]]))
-  diff_col2 <- max(abs(score - pred[[2]]))
-
-  expect_lt(min(diff_col1, diff_col2), 0.03)
-})
-
-test_that("rpf binary: path-dependent sum identity matches predicted class probability", {
-  skip_if_not_installed("randomPlantedForest")
-  skip_on_os("windows") # rpf purify_3() OOB read, see comment at top of file
+  gl <- glex::glex(rp, xdat)
 
-  rp <- randomPlantedForest::rpf(
-    y ~ x1 + x2 + x3,
-    data = xdat,
-    max_interaction = 3
+  expect_equal(
+    unname(gl$intercept + rowSums(gl$m)),
+    unname(predict(rp, xdat, type = "numeric")[[1]]),
+    tolerance = 1e-8
   )
-  gl <- glex::glex(rp, xdat, weighting_method = "path-dependent")
-  pred <- predict(rp, xdat)
-
-  score <- unname(gl$intercept + rowSums(gl$m))
-  diff_col1 <- max(abs(score - pred[[1]]))
-  diff_col2 <- max(abs(score - pred[[2]]))
-
-  expect_lt(min(diff_col1, diff_col2), 0.03)
 })
 
-test_that("rpf multiclass: classwise sum identity matches predicted probabilities", {
+test_that("rpf multiclass: classwise sum identity holds up to the class intercept", {
   skip_if_not_installed("randomPlantedForest")
   skip_on_os("windows") # rpf purify_3() OOB read, see comment at top of file
 
@@ -52,23 +45,27 @@ test_that("rpf multiclass: classwise sum identity matches predicted probabilitie
     max_interaction = 3
   )
   gl <- glex::glex(rp, xdat)
-  pred <- predict(rp, xdat)
-
-  # Each class has its own decomposition terms ending in "__class:<level>".
-  # Summing those terms plus intercept should reconstruct that class score.
-  class_diffs <- vapply(
-    gl$target_levels,
-    function(level) {
-      idx <- grepl(paste0("__class:", level), names(gl$m), fixed = TRUE)
-      score <- unname(
-        gl$intercept + rowSums(as.matrix(gl$m[, idx, with = FALSE]))
-      )
-      max(abs(score - pred[[paste0(".pred_", level)]]))
-    },
-    FUN.VALUE = numeric(1)
-  )
+  pred <- predict(rp, xdat, type = "numeric")
+
+  # Each class has its own decomposition terms ending in "__class:<level>", but rpf
+  # reports a single intercept for all classes rather than one per class. The class terms
+  # therefore reconstruct that class's raw score exactly up to an additive constant -- the
+  # class intercept we never receive. Assert exactly that: the per-observation residual is
+  # constant within a class, which is much stronger than bounding its magnitude.
+  for (level in gl$target_levels) {
+    idx <- grepl(paste0("__class:", level), names(gl$m), fixed = TRUE)
+    score <- unname(
+      gl$intercept + rowSums(as.matrix(gl$m[, idx, with = FALSE]))
+    )
+    residual <- score - pred[[paste0(".pred_", level)]]
 
-  expect_true(all(class_diffs < 0.8))
+    expect_equal(
+      max(residual) - min(residual),
+      0,
+      tolerance = 1e-8,
+      info = paste("class", level)
+    )
+  }
 })
 
 test_that("rpf: shap is derived from components and satisfies efficiency", {

From a9caae4e0616167fd62057df630a8191f3fb74cc Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Tue, 14 Jul 2026 02:37:25 +0200
Subject: [PATCH 08/10] robustify tests

---
 R/glex.R                                    | 11 ++--
 tests/testthat/test-constrained-remainder.R | 23 ++++---
 tests/testthat/test-rpf-sum-identity.R      | 68 +++++++++++++++------
 3 files changed, 69 insertions(+), 33 deletions(-)

diff --git a/R/glex.R b/R/glex.R
index 35fe77f..d60e3a0 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -751,11 +751,12 @@ confirm_constrained <- function(res, target = NULL) {
     reconstruction <- res$intercept + rowSums(res$m)
     res$remainder <- unname(target - reconstruction)
 
-    inert <- isTRUE(all.equal(
-      unname(reconstruction),
-      unname(target),
-      tolerance = 1e-5
-    ))
+    # The dropped terms are inert only if they are *numerically zero*, which is what the
+    # remainder measures directly. Judge it elementwise: `all.equal()` reports the mean
+    # relative difference, which averages a discrepancy concentrated in a few observations
+    # away to nothing and lets a genuinely non-zero term pass as "all zero" -- and this
+    # verdict decides whether `$shap` is trustworthy or `NA`, so it must not be fuzzy.
+    inert <- max(abs(res$remainder)) <= 1e-8 * max(1, max(abs(unname(target))))
 
     if (inert) {
       message(
diff --git a/tests/testthat/test-constrained-remainder.R b/tests/testthat/test-constrained-remainder.R
index 0927896..2874719 100644
--- a/tests/testthat/test-constrained-remainder.R
+++ b/tests/testthat/test-constrained-remainder.R
@@ -262,25 +262,30 @@ test_that("an inert constraint leaves no remainder", {
   skip_if_not_installed("randomPlantedForest")
   skip_on_os("windows") # rpf purification OOB read, see test-rpf-sum-identity.R
 
-  # A model fit at the maximum order whose top-order term is exactly zero: dropping it
-  # changes nothing, so the decomposition is still complete and owes no remainder.
+  # A constant predictor cannot be split on, so every term involving it is exactly zero --
+  # on every platform, unlike a high-order term that merely happens to come out zero for a
+  # given fit. Dropping only those terms changes nothing, so the decomposition is still
+  # complete and owes no remainder. See test-rpf-sum-identity.R for the full story.
   set.seed(42)
   n <- 120
-  p <- 6
-  x <- matrix(rnorm(n * p), ncol = p)
-  colnames(x) <- paste0("x", seq_len(p))
-  y <- x[, 1] + x[, 2] * x[, 3] + rnorm(n, sd = 0.3)
-  dat <- data.frame(x, y = y)
+  dat <- data.frame(
+    x1 = rnorm(n),
+    x2 = rnorm(n),
+    x3 = rnorm(n),
+    x4 = rep(1, n)
+  )
+  dat$y <- dat$x1 + dat$x2 * dat$x3 + rnorm(n, sd = 0.3)
+  predictors <- dat[, c("x1", "x2", "x3", "x4")]
 
   rp <- randomPlantedForest::rpf(
     y ~ .,
     data = dat,
-    max_interaction = p,
+    max_interaction = 4,
     ntrees = 20
   )
 
   expect_message(
-    gl <- glex(rp, dat[, -ncol(dat)], max_interaction = p - 1L),
+    gl <- glex(rp, predictors, max_interaction = 3),
     "dropped terms are all zero"
   )
   expect_identical(gl$constrained, character(0))
diff --git a/tests/testthat/test-rpf-sum-identity.R b/tests/testthat/test-rpf-sum-identity.R
index 11b939b..751e751 100644
--- a/tests/testthat/test-rpf-sum-identity.R
+++ b/tests/testthat/test-rpf-sum-identity.R
@@ -155,41 +155,71 @@ test_that("rpf: constraints that drop only zero terms keep shap valid", {
   skip_if_not_installed("randomPlantedForest")
   skip_on_os("windows") # rpf purify_3() OOB read, see comment at top of file
 
+  # The inert term has to be zero *by construction*, not by luck of the fit. An earlier
+  # version of this test fit at the maximum order and assumed the top-order term came out
+  # exactly zero because the data had no structure at that order; it usually did, but the
+  # forest is not bit-identical across platforms, and on macOS that term came out at 4e-4
+  # -- non-zero, so dropping it moved the SHAP values and the premise collapsed.
+  #
+  # A constant predictor cannot be split on, so *every* term involving it is exactly zero
+  # on every platform. With four features, the sole degree-4 term must involve it too,
+  # which makes both constraint axes deterministically inert.
   set.seed(42)
   n <- 120
-  p <- 6
-  x <- matrix(rnorm(n * p), ncol = p)
-  colnames(x) <- paste0("x", seq_len(p))
-  y <- x[, 1] + x[, 2] * x[, 3] + rnorm(n, sd = 0.3)
-  dat <- data.frame(x, y = y)
-
-  # Fit at the maximum order: the model contains a p-way term, but the data has no
-  # p-way structure, so that term is exactly zero.
+  dat <- data.frame(
+    x1 = rnorm(n),
+    x2 = rnorm(n),
+    x3 = rnorm(n),
+    x4 = rep(1, n)
+  )
+  dat$y <- dat$x1 + dat$x2 * dat$x3 + rnorm(n, sd = 0.3)
+  predictors <- dat[, c("x1", "x2", "x3", "x4")]
+
   rp <- randomPlantedForest::rpf(
     y ~ .,
     data = dat,
-    max_interaction = p,
+    max_interaction = 4,
     ntrees = 20
   )
-  full <- glex::glex(rp, dat[, -ncol(dat)])
-  degree <- lengths(strsplit(names(full$m), ":", fixed = TRUE))
-  expect_equal(max(abs(as.matrix(full$m[, degree == p, with = FALSE]))), 0)
+  full <- glex::glex(rp, predictors)
 
-  # Dropping only that zero term changes nothing, so shap must survive -- with a
-  # message rather than a warning, since a constraint *was* requested.
+  involves_x4 <- grepl("(^|:)x4($|:)", names(full$m))
+  expect_gt(sum(involves_x4), 0)
+  expect_equal(
+    max(abs(as.matrix(full$m[, involves_x4, with = FALSE]))),
+    0
+  )
+
+  # `max_interaction`: the only degree-4 term contains x4, so it is zero and dropping it
+  # changes nothing. A message rather than a warning, since a constraint *was* requested.
   expect_message(
-    gl <- glex::glex(rp, dat[, -ncol(dat)], max_interaction = p - 1L),
+    gl_mi <- glex::glex(rp, predictors, max_interaction = 3),
     "dropped terms are all zero"
   )
-  expect_identical(gl$constrained, character(0))
-  expect_false(anyNA(gl$shap))
-  expect_equal(as.matrix(gl$shap), as.matrix(full$shap), tolerance = 1e-8)
+  expect_identical(gl_mi$constrained, character(0))
+  expect_null(gl_mi$remainder)
+  expect_false(anyNA(gl_mi$shap))
+  expect_equal(as.matrix(gl_mi$shap), as.matrix(full$shap), tolerance = 1e-8)
+
+  # `features`: dropping x4 drops only its (zero) terms
+  expect_message(
+    gl_ft <- glex::glex(rp, predictors, features = c("x1", "x2", "x3")),
+    "dropped terms are all zero"
+  )
+  expect_identical(gl_ft$constrained, character(0))
+  expect_null(gl_ft$remainder)
+  expect_equal(
+    as.matrix(gl_ft$shap),
+    as.matrix(full$shap[, c("x1", "x2", "x3")]),
+    tolerance = 1e-8
+  )
 
   # Dropping terms that carry weight still invalidates
   expect_warning(
-    gl_real <- glex::glex(rp, dat[, -ncol(dat)], max_interaction = 1),
+    gl_real <- glex::glex(rp, predictors, max_interaction = 1),
     "efficiency property"
   )
   expect_identical(gl_real$constrained, "max_interaction")
   expect_true(all(is.na(gl_real$shap)))
+  expect_false(is.null(gl_real$remainder))
 })

From a3d9d01ad53b27af209228e984d8038eb6c6af85 Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Sat, 18 Jul 2026 01:56:17 +0200
Subject: [PATCH 09/10] constrained $shap is scalar NA

---
 NEWS.md                                | 2 +-
 R/glex.R                               | 8 ++++----
 man/glex.Rd                            | 2 +-
 tests/testthat/test-glex-xgboost.R     | 6 +++---
 tests/testthat/test-mtcars-ranger.R    | 4 ++--
 tests/testthat/test-rpf-sum-identity.R | 6 +++---
 6 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 059d0e0..e25acfc 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,6 @@
 # glex 0.6.0.9000 (development version)
 
-* `$shap` is now `NA` (with a warning) when the decomposition is constrained via
+* `$shap` is now a scalar `NA` (with a warning) when the decomposition is constrained via
   `max_interaction` or `features`: a constrained decomposition does not sum to the
   full model prediction, so SHAP values cannot be reconstructed from it without
   violating the efficiency property. Previously, misleading values were returned.
diff --git a/R/glex.R b/R/glex.R
index d60e3a0..3e692ad 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -24,7 +24,7 @@
 #'   \eqn{\phi_j = \sum_{S \ni j} m_S / |S|}. This reconstruction is only valid if the
 #'   decomposition is complete: if it is constrained (see `constrained`), the components
 #'   no longer sum to the full model prediction and the SHAP efficiency property cannot
-#'   hold, so `shap` is set to `NA` (with a warning) while `m` remains valid.
+#'   hold, so `shap` is a scalar `NA` (with a warning) while `m` remains valid.
 #'   For multiclass models, columns are class-specific like those of `m`. Note that
 #'   `randomPlantedForest` models report a single `intercept` for all classes, so for
 #'   multiclass models `intercept + rowSums(shap)` reconstructs the predicted class
@@ -778,9 +778,9 @@ confirm_constrained <- function(res, target = NULL) {
     "reconstructed from it (the efficiency property would be violated). ",
     "The components in `$m` are unaffected."
   )
-  for (col in names(res$shap)) {
-    data.table::set(res$shap, j = col, value = NA_real_)
-  }
+  # Scalar: there are no per-feature values to report, and code that consumes
+  # `$shap` numerically should fail loudly rather than propagate NAs.
+  res$shap <- NA
 
   res
 }
diff --git a/man/glex.Rd b/man/glex.Rd
index 8c2d30d..8ab2d05 100644
--- a/man/glex.Rd
+++ b/man/glex.Rd
@@ -61,7 +61,7 @@ A \code{list} with elements:
 \eqn{\phi_j = \sum_{S \ni j} m_S / |S|}. This reconstruction is only valid if the
 decomposition is complete: if it is constrained (see \code{constrained}), the components
 no longer sum to the full model prediction and the SHAP efficiency property cannot
-hold, so \code{shap} is set to \code{NA} (with a warning) while \code{m} remains valid.
+hold, so \code{shap} is a scalar \code{NA} (with a warning) while \code{m} remains valid.
 For multiclass models, columns are class-specific like those of \code{m}. Note that
 \code{randomPlantedForest} models report a single \code{intercept} for all classes, so for
 multiclass models \code{intercept + rowSums(shap)} reconstructs the predicted class
diff --git a/tests/testthat/test-glex-xgboost.R b/tests/testthat/test-glex-xgboost.R
index c8c21c2..01569cd 100644
--- a/tests/testthat/test-glex-xgboost.R
+++ b/tests/testthat/test-glex-xgboost.R
@@ -61,7 +61,7 @@ test_that("shap is NA when decomposition is constrained", {
     "efficiency property"
   )
   expect_identical(gl_mi$constrained, "max_interaction")
-  expect_true(all(is.na(gl_mi$shap)))
+  expect_identical(gl_mi$shap, NA)
   expect_false(anyNA(gl_mi$m))
 
   # constrained via features
@@ -70,7 +70,7 @@ test_that("shap is NA when decomposition is constrained", {
     "efficiency property"
   )
   expect_identical(gl_ft$constrained, "features")
-  expect_true(all(is.na(gl_ft$shap)))
+  expect_identical(gl_ft$shap, NA)
 
   # both axes are reported
   expect_warning(
@@ -78,7 +78,7 @@ test_that("shap is NA when decomposition is constrained", {
     "efficiency property"
   )
   expect_identical(gl_both$constrained, c("max_interaction", "features"))
-  expect_true(all(is.na(gl_both$shap)))
+  expect_identical(gl_both$shap, NA)
 
   # unconstrained: shap present and satisfies the efficiency property
   gl_full <- glex(xg, x)
diff --git a/tests/testthat/test-mtcars-ranger.R b/tests/testthat/test-mtcars-ranger.R
index fb79692..31965e6 100644
--- a/tests/testthat/test-mtcars-ranger.R
+++ b/tests/testthat/test-mtcars-ranger.R
@@ -108,7 +108,7 @@ test_that("ranger: constrained decomposition invalidates shap", {
     "efficiency property"
   )
   expect_identical(gl_mi$constrained, "max_interaction")
-  expect_true(all(is.na(gl_mi$shap)))
+  expect_identical(gl_mi$shap, NA)
   expect_false(anyNA(gl_mi$m))
 
   expect_warning(
@@ -116,7 +116,7 @@ test_that("ranger: constrained decomposition invalidates shap", {
     "efficiency property"
   )
   expect_identical(gl_ft$constrained, "features")
-  expect_true(all(is.na(gl_ft$shap)))
+  expect_identical(gl_ft$shap, NA)
   expect_false(anyNA(gl_ft$m))
 })
 
diff --git a/tests/testthat/test-rpf-sum-identity.R b/tests/testthat/test-rpf-sum-identity.R
index 751e751..fda71ce 100644
--- a/tests/testthat/test-rpf-sum-identity.R
+++ b/tests/testthat/test-rpf-sum-identity.R
@@ -103,7 +103,7 @@ test_that("rpf: constrained decompositions invalidate shap", {
     "efficiency property"
   )
   expect_identical(gl_mi$constrained, "max_interaction")
-  expect_true(all(is.na(gl_mi$shap)))
+  expect_identical(gl_mi$shap, NA)
   expect_false(anyNA(gl_mi$m))
 
   expect_warning(
@@ -111,7 +111,7 @@ test_that("rpf: constrained decompositions invalidate shap", {
     "efficiency property"
   )
   expect_identical(gl_ft$constrained, "features")
-  expect_true(all(is.na(gl_ft$shap)))
+  expect_identical(gl_ft$shap, NA)
   expect_false(anyNA(gl_ft$m))
 })
 
@@ -220,6 +220,6 @@ test_that("rpf: constraints that drop only zero terms keep shap valid", {
     "efficiency property"
   )
   expect_identical(gl_real$constrained, "max_interaction")
-  expect_true(all(is.na(gl_real$shap)))
+  expect_identical(gl_real$shap, NA)
   expect_false(is.null(gl_real$remainder))
 })

From 9b84d41534f81cc5fadf6f0355a6630fe836a469 Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Sat, 18 Jul 2026 02:05:40 +0200
Subject: [PATCH 10/10] formatting

---
 NEWS.md | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index e25acfc..8bc516b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -32,9 +32,9 @@
   response scale directly.
 * For `randomPlantedForest` classification models, `glex()` now confirms the constraint
   against `predict(type = "numeric")` rather than the default `type = "prob"`. rpf
-  decomposes the raw score, while `type = "prob"` applies rpf's response function --
-  a clamp to `[0, 1]` for `loss = "L2"`, the inverse link for `"logit"` and
-  `"exponential"` -- and for binary models returns the classes in an order whose first
+  decomposes the raw score, while `type = "prob"` applies rpf's response function
+  (a clamp to `[0, 1]` for `loss = "L2"`, the inverse link for `"logit"` and
+  `"exponential"`) and for binary models returns the classes in an order whose first
   column is not the one being decomposed. The components reconstruct the raw score
   exactly, so `$remainder` now measures only what the dropped terms are worth, instead of
   silently absorbing the back-transformation and the class mix-up.
@@ -44,10 +44,6 @@
   a value reconstructed from the constrained components, and for objects created by
   earlier versions of glex, which have no `$shap`.
 * `print()` on a `glex` object reports when the decomposition is constrained.
-* Tests that fit `randomPlantedForest` models remain skipped on Windows: an
-  out-of-bounds read in randomPlantedForest's purification crashes R there. Fixed
-  upstream in randomPlantedForest 0.3.0; the skips are lifted separately, once glex
-  requires that version.
 
 # glex 0.6.0