From daa144e922dda6b86c76ef5c2cf478b3972f6a8b Mon Sep 17 00:00:00 2001
From: Emil Hvitfeldt <emilhhvitfeldt@gmail.com>
Date: Sun, 23 Aug 2026 22:53:59 -0700
Subject: [PATCH] Add tidypredict_combine_trees() and the last extractors

tidypredict_trees() returns per-tree expressions, but putting them back
together is not summing or averaging them, and the rule differs by
backend:

  blackboost   offset + nu * sum(trees)
  aorsf        mean(trees), inside a guard returning NA for a partial row
  catboost     scale * sum(trees) + bias, then an inverse link
  xgboost      base_score + sum(trees), then the objective's inverse link
  lightgbm     sum(trees), divided by n for boosting = "rf", then the link
  randomForest mean(trees), plus corr.bias rescaling when present
  cforest      mean(trees)

So exposing the trees alone was not enough. A caller doing the obvious
thing would drop blackboost's offset and shrinkage, and aorsf's missing-
row guard, and get plausible wrong numbers rather than an error.

Add tidypredict_combine_trees(x, trees), which holds that rule. `trees`
may be expressions or symbols naming columns the trees were written to,
since computing each tree into its own column and then combining
references to them is the case this exists for.

Every ensemble's tidypredict_fit() now routes through the two generics,
so the identity

  tidypredict_combine_trees(x, tidypredict_trees(x)) == tidypredict_fit(x)

holds by construction rather than by coincidence. The tests assert it at
the value level for all nine backends, and again via symbols.

Also adds tidypredict_trees() and tidypredict_n_trees() for cforest,
blackboost and aorsf.

Two backends deliberately have no trees exposed. Boosted C5.0 combines
by a confidence-weighted vote over class labels, so there are no per-tree
numbers to combine at all; tidypredict_combine_trees() refuses with an
explanation rather than a bare "no method". And xrf is not a tree
ensemble: it is a penalized linear model whose predictors are rule
indicators, already handled as a pm_regression.

Multiclass lightgbm and catboost also refuse: their tidypredict_trees()
output is a flat round-robin list whose adjacent trees belong to
different classes, so it is not summable as one ensemble.

Corrects a stale comment on tidypredict_n_trees.lgb.Booster() claiming
single-leaf trees are dropped by the extractor. add_lgb_stump_trees()
deliberately restores them, because multiclass class assignment is
positional and a gap shifts every later class (#419).
---
 NAMESPACE                               |  17 ++
 NEWS.md                                 |   8 +
 R/combine-trees.R                       |  71 ++++++
 R/model-C5.0.R                          |  28 +++
 R/model-aorsf.R                         |  42 +++-
 R/model-catboost.R                      |  57 ++++-
 R/model-cforest.R                       |  32 ++-
 R/model-lightgbm.R                      |  53 +++-
 R/model-mboost.R                        |  45 +++-
 R/model-ranger.R                        |  23 +-
 R/model-rf.R                            |  42 +++-
 R/model-xgboost.R                       |  37 ++-
 man/tidypredict_combine_trees.Rd        |  58 +++++
 tests/testthat/_snaps/combine-trees.md  |  33 +++
 tests/testthat/_snaps/model-catboost.md |   2 +-
 tests/testthat/_snaps/model-lightgbm.md |   6 +-
 tests/testthat/test-combine-trees.R     | 305 ++++++++++++++++++++++++
 17 files changed, 812 insertions(+), 47 deletions(-)
 create mode 100644 R/combine-trees.R
 create mode 100644 man/tidypredict_combine_trees.Rd
 create mode 100644 tests/testthat/_snaps/combine-trees.md
 create mode 100644 tests/testthat/test-combine-trees.R

diff --git a/NAMESPACE b/NAMESPACE
index b441dfcc..c7774909 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -72,6 +72,16 @@ S3method(tidypredict_class_exprs,rpart)
 S3method(tidypredict_class_trees,default)
 S3method(tidypredict_class_trees,randomForest)
 S3method(tidypredict_class_trees,ranger)
+S3method(tidypredict_combine_trees,C5.0)
+S3method(tidypredict_combine_trees,ObliqueForest)
+S3method(tidypredict_combine_trees,blackboost)
+S3method(tidypredict_combine_trees,catboost.Model)
+S3method(tidypredict_combine_trees,cforest)
+S3method(tidypredict_combine_trees,default)
+S3method(tidypredict_combine_trees,lgb.Booster)
+S3method(tidypredict_combine_trees,randomForest)
+S3method(tidypredict_combine_trees,ranger)
+S3method(tidypredict_combine_trees,xgb.Booster)
 S3method(tidypredict_fit,"_xgb.Booster")
 S3method(tidypredict_fit,C5.0)
 S3method(tidypredict_fit,H2OBinomialModel)
@@ -119,7 +129,10 @@ S3method(tidypredict_interval,default)
 S3method(tidypredict_interval,glm)
 S3method(tidypredict_interval,list)
 S3method(tidypredict_interval,lm)
+S3method(tidypredict_n_trees,ObliqueForest)
+S3method(tidypredict_n_trees,blackboost)
 S3method(tidypredict_n_trees,catboost.Model)
+S3method(tidypredict_n_trees,cforest)
 S3method(tidypredict_n_trees,default)
 S3method(tidypredict_n_trees,lgb.Booster)
 S3method(tidypredict_n_trees,randomForest)
@@ -226,7 +239,10 @@ S3method(tidypredict_test,rpart)
 S3method(tidypredict_test,sda)
 S3method(tidypredict_test,xgb.Booster)
 S3method(tidypredict_test,xrf)
+S3method(tidypredict_trees,ObliqueForest)
+S3method(tidypredict_trees,blackboost)
 S3method(tidypredict_trees,catboost.Model)
+S3method(tidypredict_trees,cforest)
 S3method(tidypredict_trees,default)
 S3method(tidypredict_trees,lgb.Booster)
 S3method(tidypredict_trees,randomForest)
@@ -260,6 +276,7 @@ export(set_catboost_categories)
 export(tidy)
 export(tidypredict_class_exprs)
 export(tidypredict_class_trees)
+export(tidypredict_combine_trees)
 export(tidypredict_fit)
 export(tidypredict_interval)
 export(tidypredict_load)
diff --git a/NEWS.md b/NEWS.md
index dc94b96f..50c29742 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,13 @@
 # tidypredict (development version)
 
+- `tidypredict_combine_trees()` is a new generic that turns per-tree expressions back into a model's prediction. `tidypredict_trees()` alone is not enough to do this: `mboost::blackboost()` combines as `offset + nu * sum(trees)`, `aorsf` averages inside a guard that returns `NA` for an incomplete row, CatBoost applies a scale and a bias, and boosters then apply their objective's inverse link. Summing or averaging the trees, as the shape of the list invites, is wrong for all of those. (#436)
+
+- `tidypredict_combine_trees()` has methods for `randomForest`, `ranger`, xgboost, LightGBM, CatBoost, `cforest`, `blackboost` and `aorsf`, so a caller that computes each tree into its own column can combine references to those columns without knowing which backend it is holding. Every one satisfies `tidypredict_combine_trees(x, tidypredict_trees(x))` computing the same values as `tidypredict_fit(x)`. (#436)
+
+- Boosted `C50::C5.0()` models deliberately have no `tidypredict_trees()` method, and `tidypredict_combine_trees()` refuses them with an explanation. Their trials vote with a class label and a confidence rather than contributing numbers, so there is nothing to sum or average and splitting the trees apart would only enable a wrong answer. (#436)
+
+- `tidypredict_trees()` and `tidypredict_n_trees()` gain methods for `partykit::cforest()`, `mboost::blackboost()` and `aorsf::orsf()`. (#436)
+
 - New generics describe what a model's fitted expressions compute, which the expressions themselves do not say: `tidypredict_output_type()` returns one of `"numeric"`, `"prob"`, `"decision"` or `"class"`, `tidypredict_outcome_levels()` returns the outcome levels in model order, and `tidypredict_normalized()` reports whether per-level probabilities already sum to one. See `?tidypredict_metadata`. (#435)
 
 - The distinctions these generics record cannot be recovered from the shape of a `tidypredict_fit()` result. A `LiblineaR` SVM classifier and a `LiblineaR` logistic regression both return a single expression, but the first is an uncalibrated decision value whose sign picks the class, so thresholding it at 0.5 as though it were a probability is wrong. A multiclass probability list and a `quantreg::rq()` fit with several `tau` are both named lists of expressions of the same length. (#435)
diff --git a/R/combine-trees.R b/R/combine-trees.R
new file mode 100644
index 00000000..c12639e0
--- /dev/null
+++ b/R/combine-trees.R
@@ -0,0 +1,71 @@
+#' Combine per-tree expressions into a model's prediction
+#'
+#' @description
+#' [tidypredict_trees()] returns one expression per tree. Turning those back
+#' into the model's prediction is not simply summing or averaging them, and the
+#' rule differs by backend: `mboost::blackboost()` needs an offset and a
+#' shrinkage factor, CatBoost needs a scale and a bias, `aorsf` needs a guard
+#' that returns `NA` for an incomplete row, and boosters then apply their
+#' objective's inverse link on top.
+#'
+#' This generic holds that rule, so a caller that has split the trees apart
+#' can put them back together without knowing which backend it is holding.
+#'
+#' @param x A fitted model object.
+#' @param trees A list of expressions, one per tree, in the order
+#'   [tidypredict_trees()] returns them. Typically either that return value
+#'   itself, or symbols naming the columns the individual trees were written
+#'   to.
+#' @param ... Additional arguments passed to methods.
+#'
+#' @returns A single language object.
+#'
+#' @details
+#' The point of separating `trees` from this function is that a caller can
+#' compute each tree into its own column, for a database to evaluate in
+#' parallel, and then pass symbols referring to those columns rather than the
+#' expressions themselves. The combination is the same either way.
+#'
+#' Every ensemble satisfies
+#' `tidypredict_combine_trees(x, tidypredict_trees(x))` computing the same
+#' values as `tidypredict_fit(x)`, and that identity is what the tests for
+#' these methods assert.
+#'
+#' Not every ensemble has a method. `C50::C5.0()` boosting combines its trees
+#' by a confidence-weighted vote that yields a class label, so there is no
+#' arithmetic to apply to per-tree numbers and no method is provided.
+#'
+#' @examplesIf rlang::is_installed("randomForest")
+#' model <- randomForest::randomForest(mpg ~ ., data = mtcars, ntree = 3)
+#'
+#' trees <- tidypredict_trees(model)
+#' tidypredict_combine_trees(model, trees)
+#'
+#' # Or referring to columns the trees were written to first.
+#' tidypredict_combine_trees(model, rlang::syms(c("t1", "t2", "t3")))
+#'
+#' @export
+tidypredict_combine_trees <- function(x, trees, ...) {
+  UseMethod("tidypredict_combine_trees")
+}
+
+#' @export
+tidypredict_combine_trees.default <- function(x, trees, ...) {
+  cli::cli_abort(
+    "{.fn tidypredict_combine_trees} is not available for models of class
+     {.cls {class(x)[[1]]}}.",
+    class = "tidypredict_no_combiner"
+  )
+}
+
+check_trees_arg <- function(trees, call = rlang::caller_env()) {
+  if (!is.list(trees) || length(trees) == 0) {
+    cli::cli_abort(
+      "{.arg trees} must be a non-empty list of expressions, not
+       {.obj_type_friendly {trees}}.",
+      call = call
+    )
+  }
+
+  invisible(trees)
+}
diff --git a/R/model-C5.0.R b/R/model-C5.0.R
index d0069754..f1ae1abf 100644
--- a/R/model-C5.0.R
+++ b/R/model-C5.0.R
@@ -1005,6 +1005,34 @@ c50_boosted_tree_info <- function(trees, classes) {
   lapply(trees, function(tree) c50_tree_info(tree, classes, tree$root_class))
 }
 
+# Extractors --------------------------------------
+
+# The per-trial trees are available, but they cannot be recombined by any
+# arithmetic on their values. Boosted C5.0 predicts by a confidence-weighted
+# vote across trials, where each tree contributes a *class label* plus a
+# confidence, and `SelectClassGen()` breaks ties using that trial's own root
+# class. Splitting the trees apart and adding or averaging them is not an
+# approximation of that, it is a different model.
+#
+# So `tidypredict_trees()` is deliberately not implemented here: exposing the
+# trees without a way to put them back together would only enable a wrong
+# answer. `tidypredict_combine_trees()` refuses for the same reason, with a
+# message saying why rather than the generic "no method" one.
+#' @export
+tidypredict_combine_trees.C5.0 <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+
+  cli::cli_abort(
+    c(
+      "Boosted {.pkg C5.0} trees cannot be recombined arithmetically.",
+      i = "Each trial votes with a class label and a confidence, so there are
+           no per-tree numbers to sum or average.",
+      i = "Use {.fn tidypredict_fit} for the whole model instead."
+    ),
+    class = "tidypredict_no_combiner"
+  )
+}
+
 # Test ---------------------------------------------
 
 # `C5.0()` keeps no copy of the training data, so `df` has no useful default.
diff --git a/R/model-aorsf.R b/R/model-aorsf.R
index 110418d1..c706e51c 100644
--- a/R/model-aorsf.R
+++ b/R/model-aorsf.R
@@ -108,17 +108,45 @@ parse_model.ObliqueForest <- function(model) {
 tidypredict_fit.ObliqueForest <- function(model, ...) {
   aorsf_check_supported(model)
 
-  n_trees <- length(model$forest$child_left)
-  tree_exprs <- map(
-    seq_len(n_trees),
+  tidypredict_combine_trees(model, tidypredict_trees(model))
+}
+
+# Extractors --------------------------------------
+
+#' @export
+tidypredict_trees.ObliqueForest <- function(x, ...) {
+  rlang::check_dots_empty()
+  aorsf_check_supported(x)
+
+  map(
+    seq_len(length(x$forest$child_left)),
     function(tree_no) {
-      generate_nested_case_when_tree(aorsf_tree_info_full(model, tree_no))
+      generate_nested_case_when_tree(aorsf_tree_info_full(x, tree_no))
     }
   )
+}
 
-  # `aorsf` refuses to predict from an incomplete row ("Please remove missing
-  # values from new data, or impute them."), so there is no value to match.
-  expr_na_if_incomplete(expr_mean(tree_exprs, n_trees), model$get_names_x())
+#' @export
+tidypredict_n_trees.ObliqueForest <- function(x, ...) {
+  rlang::check_dots_empty()
+
+  length(x$forest$child_left)
+}
+
+# The mean is only half of it. `aorsf` refuses to predict from an incomplete
+# row ("Please remove missing values from new data, or impute them."), so the
+# result is wrapped to return `NA` for such a row. A caller that averaged the
+# trees itself would silently produce a number where the model produces
+# nothing.
+#' @export
+tidypredict_combine_trees.ObliqueForest <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+  check_trees_arg(trees)
+
+  expr_na_if_incomplete(
+    expr_mean(trees, length(trees)),
+    x$get_names_x()
+  )
 }
 
 build_tree_formula.pm_tree_aorsf <- function(model) {
diff --git a/R/model-catboost.R b/R/model-catboost.R
index e99c3fc8..7cd1b279 100644
--- a/R/model-catboost.R
+++ b/R/model-catboost.R
@@ -717,7 +717,18 @@ build_fit_formula_catboost_nested <- function(parsedmodel) {
     cli::cli_abort("Model has no trees.")
   }
 
-  objective <- parsedmodel$general$params$objective %||% "RMSE"
+  objective <- catboost_parsed_objective(parsedmodel)
+  catboost_check_objective(parsedmodel)
+
+  if (objective %in% catboost_multiclass_objectives) {
+    return(build_fit_formula_catboost_multiclass_nested(parsedmodel, objective))
+  }
+
+  catboost_combine(extract_catboost_trees_nested(parsedmodel), parsedmodel)
+}
+
+catboost_check_objective <- function(parsedmodel) {
+  objective <- catboost_parsed_objective(parsedmodel)
 
   if (!objective %in% catboost_supported_objectives) {
     cli::cli_abort(
@@ -728,14 +739,15 @@ build_fit_formula_catboost_nested <- function(parsedmodel) {
     )
   }
 
-  if (objective %in% catboost_multiclass_objectives) {
-    return(build_fit_formula_catboost_multiclass_nested(parsedmodel, objective))
-  }
+  invisible(parsedmodel)
+}
 
-  # Extract nested trees
-  trees <- extract_catboost_trees_nested(parsedmodel)
-  f <- reduce_addition(trees)
-  f <- apply_catboost_scale_bias(f, parsedmodel)
+# Combine per-tree expressions into the model's prediction: an additive sum
+# rescaled by the model's `scale` and `bias`, then the objective's inverse link.
+catboost_combine <- function(trees, parsedmodel) {
+  objective <- catboost_parsed_objective(parsedmodel)
+
+  f <- apply_catboost_scale_bias(reduce_addition(trees), parsedmodel)
 
   if (objective %in% catboost_sigmoid_objectives) {
     f <- expr_logistic(f)
@@ -911,6 +923,35 @@ tidypredict_trees.catboost.Model <- function(x, ...) {
   extract_catboost_trees_nested(parse_model(x))
 }
 
+#' @export
+tidypredict_combine_trees.catboost.Model <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+  check_trees_arg(trees)
+
+  parsedmodel <- parse_model(x)
+  catboost_check_objective(parsedmodel)
+
+  # A multiclass fit is one expression per class rather than one for the model.
+  # `tidypredict_trees()` returns the trees of such a model as a flat list in
+  # which they belong to different classes round-robin, so summing them is not
+  # an approximation of the fit.
+  if (
+    catboost_parsed_objective(parsedmodel) %in% catboost_multiclass_objectives
+  ) {
+    cli::cli_abort(
+      c(
+        "Multiclass {.pkg catboost} trees cannot be combined into one
+         expression.",
+        i = "The fit is one expression per class.",
+        i = "Use {.fn tidypredict_fit} for the whole model instead."
+      ),
+      class = "tidypredict_no_combiner"
+    )
+  }
+
+  catboost_combine(trees, parsedmodel)
+}
+
 #' @export
 tidypredict_n_trees.catboost.Model <- function(x, ...) {
   rlang::check_dots_empty()
diff --git a/R/model-cforest.R b/R/model-cforest.R
index 1e232624..f6c6938f 100644
--- a/R/model-cforest.R
+++ b/R/model-cforest.R
@@ -55,16 +55,38 @@ parse_model.cforest <- function(model) {
 tidypredict_fit.cforest <- function(model, ...) {
   cforest_check_regression(model)
 
-  n_trees <- length(model$nodes)
-  tree_exprs <- map(
-    seq_len(n_trees),
+  tidypredict_combine_trees(model, tidypredict_trees(model))
+}
+
+# Extractors --------------------------------------
+
+#' @export
+tidypredict_trees.cforest <- function(x, ...) {
+  rlang::check_dots_empty()
+  cforest_check_regression(x)
+
+  map(
+    seq_len(length(x$nodes)),
     function(tree_no) {
-      tree_info <- partykit_tree_info_full(cforest_gettree(model, tree_no))
+      tree_info <- partykit_tree_info_full(cforest_gettree(x, tree_no))
       generate_nested_case_when_tree(tree_info, missing = "na")
     }
   )
+}
+
+#' @export
+tidypredict_n_trees.cforest <- function(x, ...) {
+  rlang::check_dots_empty()
+
+  length(x$nodes)
+}
+
+#' @export
+tidypredict_combine_trees.cforest <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+  check_trees_arg(trees)
 
-  expr_mean(tree_exprs, n_trees)
+  expr_mean(trees, length(trees))
 }
 
 build_tree_formula.pm_tree_cforest <- function(model) {
diff --git a/R/model-lightgbm.R b/R/model-lightgbm.R
index 0b0ac432..1c1564b5 100644
--- a/R/model-lightgbm.R
+++ b/R/model-lightgbm.R
@@ -683,7 +683,13 @@ assemble_lgb_formula <- function(parsedmodel, build_trees) {
     cli::cli_abort("Model has no trees.")
   }
 
-  objective <- parsedmodel$general$params$objective %||% "regression"
+  lgb_check_objective(parsedmodel)
+
+  lgb_combine(build_trees(), parsedmodel)
+}
+
+lgb_check_objective <- function(parsedmodel) {
+  objective <- lgb_parsed_objective(parsedmodel)
 
   if (!objective %in% lgb_supported_objectives) {
     cli::cli_abort(
@@ -694,7 +700,16 @@ assemble_lgb_formula <- function(parsedmodel, build_trees) {
     )
   }
 
-  trees <- build_trees()
+  invisible(parsedmodel)
+}
+
+# Combine per-tree expressions into the booster's prediction: an additive sum
+# (averaged instead when boosting is random forest), then the objective's
+# inverse link. A multiclass model instead groups the trees by class and
+# returns one expression per class, which is not a single language object.
+lgb_combine <- function(trees, parsedmodel) {
+  objective <- lgb_parsed_objective(parsedmodel)
+  n_trees <- length(trees)
 
   # A model of stumps mentions no column, so anchor it to one. The feature
   # names recorded at parse time are the columns `newdata` has to supply.
@@ -903,12 +918,42 @@ tidypredict_trees.lgb.Booster <- function(x, ...) {
   extract_lgb_trees_nested(x)
 }
 
+#' @export
+tidypredict_combine_trees.lgb.Booster <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+  check_trees_arg(trees)
+
+  parsedmodel <- parse_model(x)
+  lgb_check_objective(parsedmodel)
+
+  # A multiclass fit is one expression per class rather than one for the model,
+  # so there is no single language object to return. `tidypredict_trees()` does
+  # hand back the trees of such a model, and they are assigned to classes
+  # positionally, so a caller could otherwise sum trees belonging to different
+  # classes together.
+  if (lgb_parsed_objective(parsedmodel) %in% lgb_multiclass_objectives) {
+    cli::cli_abort(
+      c(
+        "Multiclass {.pkg lightgbm} trees cannot be combined into one
+         expression.",
+        i = "The fit is one expression per class.",
+        i = "Use {.fn tidypredict_fit} for the whole model instead."
+      ),
+      class = "tidypredict_no_combiner"
+    )
+  }
+
+  lgb_combine(trees, parsedmodel)
+}
+
 #' @export
 tidypredict_n_trees.lgb.Booster <- function(x, ...) {
   rlang::check_dots_empty()
 
-  # Trees with a single leaf are dropped by the extractor, so this counts the
-  # trees actually returned rather than the number LightGBM reports.
+  # Counts what the extractor returns, which includes single-leaf trees:
+  # `lgb.model.dt.tree()` omits them but `add_lgb_stump_trees()` puts them
+  # back, deliberately, because multiclass class assignment is positional and a
+  # gap shifts every later class (#419).
   length(tidypredict_trees(x))
 }
 
diff --git a/R/model-mboost.R b/R/model-mboost.R
index 7fb7b6ab..8f564c14 100644
--- a/R/model-mboost.R
+++ b/R/model-mboost.R
@@ -62,8 +62,11 @@ mboost_build_formula <- function(tree_info_list, nu, offset) {
     tree_info_list,
     \(tree_info) generate_nested_case_when_tree(tree_info, missing = "na")
   )
-  res <- reduce_addition(tree_exprs)
-  expr(!!offset + !!nu * !!res)
+  mboost_combine(tree_exprs, nu, offset)
+}
+
+mboost_combine <- function(tree_exprs, nu, offset) {
+  expr(!!offset + !!nu * !!reduce_addition(tree_exprs))
 }
 
 # Model parser -------------------------------------
@@ -90,8 +93,42 @@ parse_model.blackboost <- function(model) {
 tidypredict_fit.blackboost <- function(model, ...) {
   mboost_check_regression(model)
 
-  comps <- mboost_components(model)
-  mboost_build_formula(comps$tree_info_list, comps$nu, comps$offset)
+  tidypredict_combine_trees(model, tidypredict_trees(model))
+}
+
+# Extractors --------------------------------------
+
+#' @export
+tidypredict_trees.blackboost <- function(x, ...) {
+  rlang::check_dots_empty()
+  mboost_check_regression(x)
+
+  map(
+    mboost_components(x)$tree_info_list,
+    \(tree_info) generate_nested_case_when_tree(tree_info, missing = "na")
+  )
+}
+
+#' @export
+tidypredict_n_trees.blackboost <- function(x, ...) {
+  rlang::check_dots_empty()
+
+  # `mstop` rather than the raw ensemble length: subsetting a fitted model
+  # leaves `ens` at full length, and `mboost_components()` truncates to the
+  # iterations the model actually uses.
+  length(mboost_components(x)$tree_info_list)
+}
+
+# Boosting, so the trees are summed rather than averaged, then shrunk by `nu`
+# and offset. Summing them plainly, as a caller might assume, would be wrong on
+# both counts.
+#' @export
+tidypredict_combine_trees.blackboost <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+  check_trees_arg(trees)
+
+  comps <- mboost_components(x)
+  mboost_combine(trees, comps$nu, comps$offset)
 }
 
 build_tree_formula.pm_tree_blackboost <- function(model) {
diff --git a/R/model-ranger.R b/R/model-ranger.R
index ad0b0a6d..f6f59420 100644
--- a/R/model-ranger.R
+++ b/R/model-ranger.R
@@ -254,12 +254,14 @@ tidypredict_fit_ranger_nested <- function(model) {
     build_nested_ranger_tree(model, tree_no)
   })
 
-  # A forest of stumps mentions no column, so anchor it to one. These are the
-  # predictors `ranger:::predict.ranger()` itself requires in `newdata`.
-  expr_recycle_over_column(
-    expr_mean(tree_exprs, n_trees),
-    model$forest$independent.variable.names
-  )
+  ranger_combine(tree_exprs, n_trees, model$forest$independent.variable.names)
+}
+
+# A forest of stumps mentions no column, so the average is anchored to one.
+# These are the predictors `ranger:::predict.ranger()` itself requires in
+# `newdata`.
+ranger_combine <- function(tree_exprs, n_trees, predictors) {
+  expr_recycle_over_column(expr_mean(tree_exprs, n_trees), predictors)
 }
 
 # Build nested case_when for a single ranger tree
@@ -419,6 +421,15 @@ tidypredict_n_trees.ranger <- function(x, ...) {
   as.integer(x$num.trees)
 }
 
+#' @export
+tidypredict_combine_trees.ranger <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+  check_trees_arg(trees)
+  ranger_check_supported(x)
+
+  ranger_combine(trees, x$num.trees, x$forest$independent.variable.names)
+}
+
 build_tree_formula.pm_tree_ranger <- function(model) {
   expr_mean(map(
     model$tree_info_list,
diff --git a/R/model-rf.R b/R/model-rf.R
index 5edaa5bc..b49b046c 100644
--- a/R/model-rf.R
+++ b/R/model-rf.R
@@ -210,11 +210,15 @@ tidypredict_fit_rf_nested <- function(model) {
     build_nested_rf_tree(model, tree_no, term_labels)
   })
 
-  # `randomForest::predict()` returns `NA` for a row with any missing
-  # predictor, so the forest average is only defined on complete rows.
+  rf_combine(tree_exprs, n_trees, unname(model$coefs), term_labels)
+}
+
+# `randomForest::predict()` returns `NA` for a row with any missing predictor,
+# so the forest average is only defined on complete rows.
+rf_combine <- function(tree_exprs, n_trees, coefs, predictors) {
   expr_na_if_incomplete(
-    rf_correct_bias(expr_mean(tree_exprs, n_trees), unname(model$coefs)),
-    term_labels
+    rf_correct_bias(expr_mean(tree_exprs, n_trees), coefs),
+    predictors
   )
 }
 
@@ -384,6 +388,36 @@ tidypredict_trees.randomForest <- function(x, ...) {
   })
 }
 
+# Averaging the trees is not the whole prediction: `predict()` declines to score
+# a row with any missing predictor, so the average is wrapped in a guard that
+# returns `NA` for such a row.
+#' @export
+tidypredict_combine_trees.randomForest <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+  check_trees_arg(trees)
+
+  if (!is.null(x$classes)) {
+    cli::cli_abort(
+      c(
+        "Classification models are not supported.",
+        i = "Each tree votes for a class, so there are no per-tree numbers to
+             average."
+      ),
+      class = "tidypredict_no_combiner"
+    )
+  }
+  # As in `tidypredict_trees()`: the bias correction applies to the forest
+  # average, and the coefficients it needs cannot be recovered from the trees.
+  if (!is.null(x$coefs)) {
+    cli::cli_abort(
+      "Models fitted with {.code corr.bias = TRUE} are not supported.",
+      class = "tidypredict_no_combiner"
+    )
+  }
+
+  rf_combine(trees, x$ntree, NULL, names(x$forest$ncat))
+}
+
 build_tree_formula.pm_tree_randomForest <- function(model) {
   expr_na_if_incomplete(
     rf_correct_bias(build_tree_formula_forest(model), model$general$coefs),
diff --git a/R/model-xgboost.R b/R/model-xgboost.R
index 9db18d59..e83ef54f 100644
--- a/R/model-xgboost.R
+++ b/R/model-xgboost.R
@@ -299,16 +299,25 @@ assemble_xgb_formula <- function(
   objective
 ) {
   # Apply DART weight_drop if present
-  trees_nested <- apply_dart_weights(trees_nested, weight_drop)
-
-  # Additive model
-  f <- reduce_addition(trees_nested)
+  xgb_combine(
+    apply_dart_weights(trees_nested, weight_drop),
+    base_score,
+    objective
+  )
+}
 
+# Combine already-weighted per-tree expressions: an additive sum, then the
+# objective's inverse link and `base_score`.
+#
+# DART weighting is deliberately not applied here. `tidypredict_trees()` folds
+# it into the per-tree expressions it returns, so doing it again would square
+# the weights.
+xgb_combine <- function(trees, base_score, objective) {
   if (is.null(base_score)) {
     base_score <- 0.5 # nocov
   }
 
-  apply_xgb_objective(f, objective, base_score)
+  apply_xgb_objective(reduce_addition(trees), objective, base_score)
 }
 
 # Build nested formula from parsed xgboost model (version 3)
@@ -595,6 +604,24 @@ tidypredict_trees.xgb.Booster <- function(x, ...) {
   unname(apply_dart_weights(trees, json_params$weight_drop))
 }
 
+#' @export
+tidypredict_combine_trees.xgb.Booster <- function(x, trees, ...) {
+  rlang::check_dots_empty()
+  check_trees_arg(trees)
+
+  json_params <- get_xgb_json_params(x)
+  params <- attr(x, "param") %||% x$params
+
+  expr_recycle_over_column(
+    xgb_combine(
+      trees,
+      json_params$base_score,
+      params$objective %||% json_params$objective
+    ),
+    xgb_feature_names(x)
+  )
+}
+
 #' @export
 tidypredict_n_trees.xgb.Booster <- function(x, ...) {
   rlang::check_dots_empty()
diff --git a/man/tidypredict_combine_trees.Rd b/man/tidypredict_combine_trees.Rd
new file mode 100644
index 00000000..b358b1dc
--- /dev/null
+++ b/man/tidypredict_combine_trees.Rd
@@ -0,0 +1,58 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/combine-trees.R
+\name{tidypredict_combine_trees}
+\alias{tidypredict_combine_trees}
+\title{Combine per-tree expressions into a model's prediction}
+\usage{
+tidypredict_combine_trees(x, trees, ...)
+}
+\arguments{
+\item{x}{A fitted model object.}
+
+\item{trees}{A list of expressions, one per tree, in the order
+\code{\link[=tidypredict_trees]{tidypredict_trees()}} returns them. Typically either that return value
+itself, or symbols naming the columns the individual trees were written
+to.}
+
+\item{...}{Additional arguments passed to methods.}
+}
+\value{
+A single language object.
+}
+\description{
+\code{\link[=tidypredict_trees]{tidypredict_trees()}} returns one expression per tree. Turning those back
+into the model's prediction is not simply summing or averaging them, and the
+rule differs by backend: \code{mboost::blackboost()} needs an offset and a
+shrinkage factor, CatBoost needs a scale and a bias, \code{aorsf} needs a guard
+that returns \code{NA} for an incomplete row, and boosters then apply their
+objective's inverse link on top.
+
+This generic holds that rule, so a caller that has split the trees apart
+can put them back together without knowing which backend it is holding.
+}
+\details{
+The point of separating \code{trees} from this function is that a caller can
+compute each tree into its own column, for a database to evaluate in
+parallel, and then pass symbols referring to those columns rather than the
+expressions themselves. The combination is the same either way.
+
+Every ensemble satisfies
+\code{tidypredict_combine_trees(x, tidypredict_trees(x))} computing the same
+values as \code{tidypredict_fit(x)}, and that identity is what the tests for
+these methods assert.
+
+Not every ensemble has a method. \code{C50::C5.0()} boosting combines its trees
+by a confidence-weighted vote that yields a class label, so there is no
+arithmetic to apply to per-tree numbers and no method is provided.
+}
+\examples{
+\dontshow{if (rlang::is_installed("randomForest")) withAutoprint(\{ # examplesIf}
+model <- randomForest::randomForest(mpg ~ ., data = mtcars, ntree = 3)
+
+trees <- tidypredict_trees(model)
+tidypredict_combine_trees(model, trees)
+
+# Or referring to columns the trees were written to first.
+tidypredict_combine_trees(model, rlang::syms(c("t1", "t2", "t3")))
+\dontshow{\}) # examplesIf}
+}
diff --git a/tests/testthat/_snaps/combine-trees.md b/tests/testthat/_snaps/combine-trees.md
new file mode 100644
index 00000000..75cf69f8
--- /dev/null
+++ b/tests/testthat/_snaps/combine-trees.md
@@ -0,0 +1,33 @@
+# check_trees_arg() rejects a non-list and an empty list
+
+    Code
+      check_trees_arg(1:3)
+    Condition
+      Error:
+      ! `trees` must be a non-empty list of expressions, not an integer vector.
+    Code
+      check_trees_arg(list())
+    Condition
+      Error:
+      ! `trees` must be a non-empty list of expressions, not an empty list.
+
+# lightgbm refuses to combine a multiclass fit
+
+    Code
+      tidypredict_combine_trees(model, tidypredict_trees(model))
+    Condition
+      Error in `tidypredict_combine_trees()`:
+      ! Multiclass lightgbm trees cannot be combined into one expression.
+      i The fit is one expression per class.
+      i Use `tidypredict_fit()` for the whole model instead.
+
+# boosted C5.0 refuses to combine its trees
+
+    Code
+      tidypredict_combine_trees(model, list(1))
+    Condition
+      Error in `tidypredict_combine_trees()`:
+      ! Boosted C5.0 trees cannot be recombined arithmetically.
+      i Each trial votes with a class label and a confidence, so there are no per-tree numbers to sum or average.
+      i Use `tidypredict_fit()` for the whole model instead.
+
diff --git a/tests/testthat/_snaps/model-catboost.md b/tests/testthat/_snaps/model-catboost.md
index 64241c8a..d401510a 100644
--- a/tests/testthat/_snaps/model-catboost.md
+++ b/tests/testthat/_snaps/model-catboost.md
@@ -3,7 +3,7 @@
     Code
       tidypredict_fit(pm)
     Condition
-      Error in `build_fit_formula_catboost_nested()`:
+      Error in `catboost_check_objective()`:
       ! Unsupported objective: "UnsupportedObjective".
       i Supported objectives: "RMSE", "MAE", "Quantile", "MAPE", "Huber", "LogCosh", "Expectile", "Poisson", "Tweedie", "Logloss", "CrossEntropy", "MultiClass", and "MultiClassOneVsAll".
 
diff --git a/tests/testthat/_snaps/model-lightgbm.md b/tests/testthat/_snaps/model-lightgbm.md
index 6cc12e36..78e373fd 100644
--- a/tests/testthat/_snaps/model-lightgbm.md
+++ b/tests/testthat/_snaps/model-lightgbm.md
@@ -25,7 +25,7 @@
     Code
       tidypredict_fit(pm)
     Condition
-      Error in `assemble_lgb_formula()`:
+      Error in `lgb_check_objective()`:
       ! Unsupported objective: "unsupported_objective".
       i Supported objectives: "regression", "regression_l2", "regression_l1", "huber", "fair", "quantile", "mape", "poisson", "gamma", "tweedie", "binary", "cross_entropy", "multiclass", and "multiclassova".
 
@@ -42,7 +42,7 @@
     Code
       tidypredict_fit(pm)
     Condition
-      Error in `assemble_lgb_formula()`:
+      Error in `lgb_combine()`:
       ! Multiclass model must have num_class >= 2.
 
 # multiclass with NULL num_class throws error
@@ -50,7 +50,7 @@
     Code
       tidypredict_fit(pm)
     Condition
-      Error in `assemble_lgb_formula()`:
+      Error in `lgb_combine()`:
       ! Multiclass model must have num_class >= 2.
 
 # build_lgb_nested_condition errors on unknown type
diff --git a/tests/testthat/test-combine-trees.R b/tests/testthat/test-combine-trees.R
new file mode 100644
index 00000000..526cd7ba
--- /dev/null
+++ b/tests/testthat/test-combine-trees.R
@@ -0,0 +1,305 @@
+# Tree expressions call `case_when()` unqualified, which only resolves if dplyr
+# is on the search path. `tests/testthat.R` attaches testthat and tidypredict
+# only, so evaluate against dplyr's namespace explicitly rather than relying on
+# some earlier test file having attached it. Without this the tests pass under
+# `devtools::test()` and fail under `R CMD check`.
+eval_tree_expr <- function(expr, data) {
+  rlang::eval_tidy(expr, data, env = asNamespace("dplyr"))
+}
+
+# The identity every method has to satisfy: recombining the per-tree
+# expressions computes what the whole-model fit computes. Compared by value
+# rather than textually, since the two may legitimately arrange the arithmetic
+# differently.
+expect_combine_matches_fit <- function(model, data) {
+  trees <- tidypredict_trees(model)
+
+  expect_equal(
+    eval_tree_expr(tidypredict_combine_trees(model, trees), data),
+    eval_tree_expr(tidypredict_fit(model), data)
+  )
+}
+
+# The same identity, but built from symbols naming columns the trees were
+# written to first. This is the case the generic exists for.
+expect_combine_symbols_match_fit <- function(model, data) {
+  trees <- tidypredict_trees(model)
+  names <- paste0("tree_", seq_along(trees))
+  values <- lapply(trees, eval_tree_expr, data = data)
+  # Recycled so that a forest of stumps, whose trees evaluate to a scalar,
+  # still yields one row per observation.
+  cols <- as.data.frame(stats::setNames(values, names))
+
+  combined <- tidypredict_combine_trees(model, rlang::syms(names))
+
+  expect_equal(
+    eval_tree_expr(combined, cbind(cols, data)),
+    eval_tree_expr(tidypredict_fit(model), data)
+  )
+}
+
+test_that("check_trees_arg() rejects a non-list and an empty list", {
+  expect_snapshot(error = TRUE, {
+    check_trees_arg(1:3)
+    check_trees_arg(list())
+  })
+})
+
+test_that("tidypredict_combine_trees() has no default combination", {
+  expect_error(
+    tidypredict_combine_trees(lm(mpg ~ wt, data = mtcars), list(1)),
+    class = "tidypredict_no_combiner"
+  )
+})
+
+test_that("randomForest trees recombine into the fit", {
+  skip_if_not_installed("randomForest")
+
+  set.seed(1)
+  model <- randomForest::randomForest(mpg ~ wt + cyl + disp, mtcars, ntree = 3)
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("randomForest refuses to combine a corr.bias fit", {
+  skip_if_not_installed("randomForest")
+
+  set.seed(1)
+  model <- randomForest::randomForest(
+    mpg ~ wt + cyl + disp,
+    mtcars,
+    ntree = 3,
+    corr.bias = TRUE
+  )
+
+  expect_error(
+    tidypredict_combine_trees(model, list(1)),
+    class = "tidypredict_no_combiner"
+  )
+})
+
+test_that("randomForest refuses to combine a classification fit", {
+  skip_if_not_installed("randomForest")
+
+  set.seed(1)
+  df <- mtcars
+  df$vs <- factor(df$vs)
+  model <- randomForest::randomForest(vs ~ wt + cyl + disp, df, ntree = 3)
+
+  expect_error(
+    tidypredict_combine_trees(model, list(1)),
+    class = "tidypredict_no_combiner"
+  )
+})
+
+test_that("ranger trees recombine into the fit", {
+  skip_if_not_installed("ranger")
+
+  set.seed(1)
+  model <- ranger::ranger(mpg ~ wt + cyl + disp, mtcars, num.trees = 3)
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("xgboost trees recombine into the fit", {
+  skip_if_not_installed("xgboost")
+
+  data <- xgboost::xgb.DMatrix(
+    as.matrix(mtcars[, -9]),
+    label = mtcars$am
+  )
+  model <- xgboost::xgb.train(
+    params = list(
+      max_depth = 2L,
+      objective = "reg:squarederror",
+      base_score = 0.5
+    ),
+    data = data,
+    nrounds = 3L,
+    verbose = 0
+  )
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("xgboost trees recombine under a non-identity objective", {
+  skip_if_not_installed("xgboost")
+
+  data <- xgboost::xgb.DMatrix(
+    as.matrix(mtcars[, -9]),
+    label = mtcars$am
+  )
+  model <- xgboost::xgb.train(
+    params = list(
+      max_depth = 2L,
+      objective = "binary:logistic",
+      base_score = 0.5
+    ),
+    data = data,
+    nrounds = 3L,
+    verbose = 0
+  )
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("xgboost DART weights are not applied twice", {
+  skip_if_not_installed("xgboost")
+
+  # Shifted off the observed values so that a 32-bit split boundary cannot
+  # route a row differently than the booster does.
+  df <- mtcars
+  df[, -9] <- df[, -9] + 0.1
+
+  data <- xgboost::xgb.DMatrix(as.matrix(df[, -9]), label = df$am)
+  model <- xgboost::xgb.train(
+    params = list(
+      max_depth = 2L,
+      objective = "reg:squarederror",
+      base_score = 0.5,
+      booster = "dart",
+      rate_drop = 0.3,
+      seed = 123
+    ),
+    data = data,
+    nrounds = 4L,
+    verbose = 0
+  )
+
+  expect_combine_matches_fit(model, df)
+  expect_combine_symbols_match_fit(model, df)
+})
+
+test_that("lightgbm trees recombine into the fit", {
+  skip_if_not_installed("lightgbm")
+
+  set.seed(123)
+  cols <- c("mpg", "cyl", "disp")
+  dtrain <- lightgbm::lgb.Dataset(
+    data.matrix(mtcars[, cols]),
+    label = mtcars$hp,
+    colnames = cols
+  )
+  model <- lightgbm::lgb.train(
+    params = list(
+      num_leaves = 4L,
+      learning_rate = 1.0,
+      objective = "regression",
+      min_data_in_leaf = 1L
+    ),
+    data = dtrain,
+    nrounds = 3L,
+    verbose = -1L
+  )
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("lightgbm refuses to combine a multiclass fit", {
+  skip_if_not_installed("lightgbm")
+
+  set.seed(123)
+  cols <- c("Sepal.Length", "Sepal.Width", "Petal.Length")
+  dtrain <- lightgbm::lgb.Dataset(
+    data.matrix(iris[, cols]),
+    label = as.integer(iris$Species) - 1L,
+    colnames = cols
+  )
+  model <- lightgbm::lgb.train(
+    params = list(
+      num_leaves = 4L,
+      objective = "multiclass",
+      num_class = 3L,
+      min_data_in_leaf = 1L
+    ),
+    data = dtrain,
+    nrounds = 2L,
+    verbose = -1L
+  )
+
+  expect_snapshot(
+    error = TRUE,
+    tidypredict_combine_trees(model, tidypredict_trees(model))
+  )
+  expect_error(
+    tidypredict_combine_trees(model, tidypredict_trees(model)),
+    class = "tidypredict_no_combiner"
+  )
+})
+
+test_that("catboost trees recombine into the fit", {
+  skip_if_not_installed("catboost")
+
+  set.seed(123)
+  cols <- c("mpg", "cyl", "disp")
+  pool <- catboost_catboost.load_pool(
+    data.matrix(mtcars[, cols]),
+    label = mtcars$hp,
+    feature_names = as.list(cols)
+  )
+  model <- catboost_catboost.train(
+    pool,
+    params = list(
+      iterations = 3L,
+      depth = 3L,
+      learning_rate = 0.5,
+      loss_function = "RMSE",
+      logging_level = "Silent"
+    )
+  )
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("cforest trees recombine into the fit", {
+  skip_if_not_installed("partykit")
+
+  set.seed(1)
+  model <- partykit::cforest(mpg ~ wt + cyl, data = mtcars, ntree = 3)
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("blackboost trees recombine into the fit", {
+  skip_if_not_installed("mboost")
+
+  model <- mboost::blackboost(
+    mpg ~ wt + cyl,
+    data = mtcars,
+    control = mboost::boost_control(mstop = 3)
+  )
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("aorsf trees recombine into the fit", {
+  skip_if_not_installed("aorsf")
+
+  set.seed(1)
+  model <- aorsf::orsf(mtcars, mpg ~ wt + cyl + disp, n_tree = 3)
+
+  expect_combine_matches_fit(model, mtcars)
+  expect_combine_symbols_match_fit(model, mtcars)
+})
+
+test_that("boosted C5.0 refuses to combine its trees", {
+  skip_if_not_installed("C50")
+
+  df <- mtcars
+  df$vs <- factor(df$vs)
+  model <- C50::C5.0(df[, c("wt", "cyl", "mpg")], df$vs, trials = 3)
+
+  expect_snapshot(error = TRUE, tidypredict_combine_trees(model, list(1)))
+  expect_error(
+    tidypredict_combine_trees(model, list(1)),
+    class = "tidypredict_no_combiner"
+  )
+})