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..8bc516b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,50 @@
+# glex 0.6.0.9000 (development version)
+
+* `$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.
+  `$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.
+* 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.
+  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.
+  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
+  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.
+
 # 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..3e692ad 100644
--- a/R/glex.R
+++ b/R/glex.R
@@ -20,12 +20,39 @@
 #'
 #' @return Decomposition of the regression or classification function.
 #' A `list` with elements:
-#' * `shap`: SHAP values (`xgboost` method only).
+#' * `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 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
+#'   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,
 #'   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.
+#'   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))`. 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")
@@ -69,6 +96,62 @@ 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
+      })
+    )
+  }
+
+  ret$shap <- data.table::setDT(as.data.frame(shap))
+
+  # 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, type = "numeric")[[1]]
+  }
+  ret <- confirm_constrained(ret, target = target)
+
   # class(ret) <- c("glex", "rpf_components", class(ret))
   ret
 }
@@ -148,6 +231,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
 }
@@ -407,6 +496,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
 }
@@ -553,6 +649,191 @@ 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. 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
+constrained_by <- function(
+  max_interaction,
+  features,
+  available_interaction,
+  model_features
+) {
+  constrained <- character(0)
+
+  # 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)) {
+    constrained <- c(constrained, "features")
+  }
+
+  constrained
+}
+
+#' 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.
+#'
+#' 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)
+  }
+
+  constrained_labels <- paste0("`", res$constrained, "`", collapse = " and ")
+
+  if (!is.null(target)) {
+    reconstruction <- res$intercept + rowSums(res$m)
+    res$remainder <- unname(target - reconstruction)
+
+    # 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(
+        "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)
+      res$remainder <- NULL
+      return(res)
+    }
+  }
+
+  warning(
+    "SHAP values set to NA: the decomposition is constrained by ",
+    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."
+  )
+  # 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
+}
+
+#' 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, so parents are processed before their children.
+#' @param trees data.table with columns Node, Yes, No, Feature_num, Tree;
+#'   leaves have `Yes = NA` and `Feature_num = 0`
+#' @keywords internal
+#' @noRd
+max_order <- function(trees) {
+  Tree <- Node <- NULL
+
+  max(vapply(
+    0:max(trees$Tree),
+    function(tree) {
+      tree_info <- trees[Tree == tree, ][order(Node)]
+
+      # 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])) {
+          # 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
+      }
+
+      order_tree
+    },
+    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 +957,17 @@ calc_components <- function(
       1L
     all_S <- get_all_subsets_cpp(sort(unique(features_num)), max_interaction)
   }
+  # 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(trees),
+    model_features = model_features
+  )
+
   # Keep only those with not more than max_interaction involved features
   d <- lengths(all_S)
 
@@ -711,32 +1003,17 @@ 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))
-  )
+  shap <- shap_from_components(m_all[, -1, drop = FALSE], colnames(x))
 
-  # 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])),
     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..d3e7724 100644
--- a/R/glex_explain.R
+++ b/R/glex_explain.R
@@ -162,13 +162,33 @@ 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`),
+  # 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[,
-    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 +241,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 +252,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 +297,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/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/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 971f1c7..8ab2d05 100644
--- a/man/glex.Rd
+++ b/man/glex.Rd
@@ -57,12 +57,39 @@ 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, 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 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
+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,
 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.
+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))}. 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
new file mode 100644
index 0000000..2874719
--- /dev/null
+++ b/tests/testthat/test-constrained-remainder.R
@@ -0,0 +1,315 @@
+# `$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
+  )
+})
+
+# 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)
+  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 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
+  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 = 4,
+    ntrees = 20
+  )
+
+  expect_message(
+    gl <- glex(rp, predictors, max_interaction = 3),
+    "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))
+})
diff --git a/tests/testthat/test-glex-xgboost.R b/tests/testthat/test-glex-xgboost.R
index 43b72fa..01569cd 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,131 @@ 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_identical(gl_mi$constrained, "max_interaction")
+  expect_identical(gl_mi$shap, NA)
+  expect_false(anyNA(gl_mi$m))
+
+  # constrained via features
+  expect_warning(
+    gl_ft <- glex(xg, x, features = c("cyl", "disp")),
+    "efficiency property"
+  )
+  expect_identical(gl_ft$constrained, "features")
+  expect_identical(gl_ft$shap, NA)
+
+  # 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"))
+  expect_identical(gl_both$shap, NA)
+
+  # 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)),
+    unname(predict(xg, x)),
+    tolerance = 1e-5
+  )
+})
+
+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..31965e6 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_identical(gl_mi$shap, NA)
+  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_identical(gl_ft$shap, NA)
+  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 3b0ff4a..fda71ce 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,60 +26,200 @@ 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]]))
+  gl <- glex::glex(rp, xdat)
 
-  expect_lt(min(diff_col1, diff_col2), 0.03)
+  expect_equal(
+    unname(gl$intercept + rowSums(gl$m)),
+    unname(predict(rp, xdat, type = "numeric")[[1]]),
+    tolerance = 1e-8
+  )
 })
 
-test_that("rpf binary: path-dependent sum identity matches predicted class probability", {
+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
 
   rp <- randomPlantedForest::rpf(
-    y ~ x1 + x2 + x3,
+    yk ~ x1 + x2 + x3,
     data = xdat,
     max_interaction = 3
   )
-  gl <- glex::glex(rp, xdat, weighting_method = "path-dependent")
-  pred <- predict(rp, xdat)
+  gl <- glex::glex(rp, xdat)
+  pred <- predict(rp, xdat, type = "numeric")
 
-  score <- unname(gl$intercept + rowSums(gl$m))
-  diff_col1 <- max(abs(score - pred[[1]]))
-  diff_col2 <- max(abs(score - pred[[2]]))
+  # 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_lt(min(diff_col1, diff_col2), 0.03)
+    expect_equal(
+      max(residual) - min(residual),
+      0,
+      tolerance = 1e-8,
+      info = paste("class", level)
+    )
+  }
 })
 
-test_that("rpf multiclass: classwise sum identity matches predicted probabilities", {
+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(
-    yk ~ x1 + x2 + x3,
-    data = xdat,
+    mpg ~ cyl + hp + wt,
+    data = mtcars,
     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)
-  )
-
-  expect_true(all(class_diffs < 0.8))
+  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),
+    "efficiency property"
+  )
+  expect_identical(gl_mi$constrained, "max_interaction")
+  expect_identical(gl_mi$shap, NA)
+  expect_false(anyNA(gl_mi$m))
+
+  expect_warning(
+    gl_ft <- glex::glex(rp, mtcars, features = c("hp", "wt")),
+    "efficiency property"
+  )
+  expect_identical(gl_ft$constrained, "features")
+  expect_identical(gl_ft$shap, NA)
+  expect_false(anyNA(gl_ft$m))
+})
+
+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)
+
+  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))
+  )
+
+  # 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)
+    )
+  }
+})
+
+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
+  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 = 4,
+    ntrees = 20
+  )
+  full <- glex::glex(rp, predictors)
+
+  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_mi <- glex::glex(rp, predictors, max_interaction = 3),
+    "dropped terms are all zero"
+  )
+  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, predictors, max_interaction = 1),
+    "efficiency property"
+  )
+  expect_identical(gl_real$constrained, "max_interaction")
+  expect_identical(gl_real$shap, NA)
+  expect_false(is.null(gl_real$remainder))
 })