diff --git a/NEWS.md b/NEWS.md
index 1abb260..137bd30 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -69,6 +69,13 @@
 
 ## New features
 
+* Classification forests with `loss = "logit"` or `"exponential"` get a level: after the fit, one constant per score column, the maximum-likelihood intercept of the log-loss given the forest's out-of-bag scores, so that every class's mean out-of-bag probability is its frequency in the training rows.
+  It is stored in every tree family's null leaf, so it is part of the intercept of the forest, of its purification and of its pruned model; `$level` holds it and `print()` shows it (`?rpf`, "The level of the classification losses").
+  The families are averaged on the score scale, and where they disagree the probability of their mean score is below the mean of their probabilities for a class rarer than 1/2, so the forest predicted a rare class too rarely, the more so the more splits it made.
+  With the level the forest predicts each class about as often as it occurs, at any number of splits.
+  The out-of-bag choice of `n_splits` scores every number of splits at its own level, and the pruning starts from the forest with its level.
+  This changes every fit with these losses, binary and multiclass; `options(randomPlantedForest.level = FALSE)` fits without the level.
+
 * `rpf()` gains `keep_level_order`, to keep the level order of factor predictors instead of ordering their levels by the outcome.
   `TRUE` keeps it for every factor, ordered or not, and a character vector for the named factors only; splits then separate neighbouring levels in the order of `levels()`.
   The default, `FALSE`, orders every factor by the outcome as before, and fits are unchanged.
@@ -153,6 +160,12 @@
 
 ## Bug fixes
 
+* `prune_components()` (and so `rpf()` with pruning) left the value of a tree family that made no split out of the pruned model's intercept, so that rule `"D"` at `tau = 0` was not the forest where a family had found no split. The intercept now includes it, as `purify()`'s does.
+
+* The out-of-bag path (`oob_path`, and the out-of-bag choice of `n_splits`) counted a tree family that had made no split at the score 0 its rows start from, instead of at the fit of its rows that the family predicts (the mean for regression and the L1 / L2 losses, the log-odds of its sample for `"logit"` and `"exponential"`, so a probability of 1/2 instead).
+  It now counts such a family at its fit.
+  Families without a split are common only with a rare class and few rows, where they made the out-of-bag log-loss of few splits too high, so that the out-of-bag choice of `n_splits` could take far more splits than it should.
+
 * `split_structure = "tree_then_leaf"` samples leaves in proportion to their available thresholds, as documented.
   After a split, the sampling weights of new leaves counted rows instead of distinct values, and the cached prefix sums were built incorrectly.
   Leaves with tied values were then drawn too often, and draws were wasted on leaves that cannot be split.
diff --git a/R/criterion.R b/R/criterion.R
index c6346b3..6a6a672 100644
--- a/R/criterion.R
+++ b/R/criterion.R
@@ -118,6 +118,247 @@ crit_subset <- function(crit, idx) {
   crit
 }
 
+# The classification level. The tree families are averaged on the score (logit) scale, and the mean of their
+# probabilities lies above the probability of their mean score where the probability is convex in the score (a rare
+# class: Jensen), more so the more the families disagree, i.e. with more splits. So a forest predicts a rare class too
+# rarely. Its level is the constant c (one value per score column) that the forest's score gets after the fit (the C++
+# shift_level(), in every family's null leaf, so part of the intercept): the maximum-likelihood intercept of the
+# log-loss given the forest's out-of-bag scores, whose equations are the balance property
+#   sum_i p_k(c + g_i) - Delta_ik = n_k   (every free class k: the mean probability is the class's frequency)
+# over the rows `idx` with an out-of-bag family, g the OOB mean scores (G; the half means G0, G1, half counts m0, m1),
+# B the families. Delta is the second-order term of the OOB correction (criterion.R, above): an OOB mean of m families
+# scatters more than the forest's mean of B, which moves its expected probability by
+#   Delta_k = (1/m - 1/B) / 2 * sum_uv d2 p_k / dg_u dg_v V_uv = (1/m - 1/B) / 2 * w p_k ((d_k - p'd)^2 - (p'd^2 - (p'd)^2))
+# with V = w d d' the across-family covariance from the halves (d over all K classes, the reference class's 0); rows
+# without both halves get none. Solved by Newton's method from c = 0 with the Jacobian of the first term,
+# sum_i diag(p_i) - p_i p_i' (positive definite), halving the step while the residual's norm does not fall; it stops at
+# a largest residual below 1e-10 of a row (usually after 3 to 5 steps). Free classes: binary the second class;
+# multiclass logit the K - 1 classes after the first (the reference, score 0); multiclass exponential (K scores, the
+# softmax) the classes after the first, the first's score then fixed at 0, and c is centred at the end (the softmax
+# ignores what all classes share; SAMME's scores are centred). 0 where a free class or the reference class has no
+# out-of-bag row (no finite level). Returns V values, with the attribute `iterations`.
+#
+# Its arithmetic is plain double arithmetic in a fixed order (level_colsums(), level_rowsums(), level_solve()), so the
+# level's bits depend on its inputs and on exp() and log() only: no BLAS or LAPACK, whose crossprod() and solve() round
+# differently by library (OpenBLAS's triangular solve multiplies by the pivot's reciprocal where the reference BLAS
+# divides by it), and no sum(), colSums() or rowSums(), which add in long double (80 bits on x86_64, 64 on arm64 macOS).
+# Then the stopping test and the step halving take the same branches on every machine too.
+classif_level <- function(crit, idx, G, G0, G1, m0, m1, B) {
+  G <- as.matrix(G)
+  V <- crit$V
+  n <- nrow(G)
+  cls <- crit$cls[idx]
+  K <- crit$K
+  zero <- structure(numeric(V), iterations = 0L)
+  if (crit$type != "logloss" || n == 0) {
+    return(zero)
+  }
+  n_k <- tabulate(cls, K)
+  if (any(n_k == 0)) {
+    return(zero)
+  }
+  # the free classes and their score columns (exponential's softmax: the first score fixed)
+  free <- seq_len(K)[-1]
+  free_col <- if (crit$link == "softmax") free else free - 1L
+  D <- as.matrix(G0) - as.matrix(G1)
+  m <- m0 + m1
+  both <- m0 > 0 & m1 > 0
+  coef <- (1 / m - 1 / B) * 0.5 / (1 / pmax(m0, 1) + 1 / pmax(m1, 1))
+  coef[!both] <- 0
+  # the full probabilities (n x K) and the score differences of the halves over all K classes (reference 0)
+  full <- function(M) {
+    if (crit$link == "softmax") M else cbind(0, M)
+  }
+  Df <- full(D)
+  probs <- function(cv) {
+    Gc <- G + rep(cv, each = n)
+    if (crit$link == "sigmoid") {
+      p <- stats::plogis(Gc[, 1])
+      return(cbind(1 - p, p))
+    }
+    Gf <- full(Gc)
+    # log-sum-exp of the rows (crit_lse()), the sum in a fixed order
+    mx <- apply(Gf, 1, max)
+    exp(Gf - (mx + log(level_rowsums(exp(Gf - mx)))))
+  }
+  resid <- function(Pf) {
+    pd <- level_rowsums(Pf * Df)
+    vp <- level_rowsums(Pf * Df^2) - pd^2
+    Delta <- coef * Pf * ((Df - pd)^2 - vp)
+    n_k[free] - level_colsums(Pf[, free, drop = FALSE] - Delta[, free, drop = FALSE])
+  }
+  # the Jacobian sum_i diag(p_i) - p_i p_i' of the free classes' probabilities P: entry (u, v) sums p_u (1[u = v] - p_v)
+  jacobian <- function(P) {
+    nf <- ncol(P)
+    uv <- which(upper.tri(diag(nf), diag = TRUE), arr.ind = TRUE)
+    terms <- vapply(
+      seq_len(nrow(uv)),
+      function(j) {
+        u <- uv[j, 1]
+        v <- uv[j, 2]
+        P[, u] * (as.numeric(u == v) - P[, v])
+      },
+      numeric(nrow(P))
+    )
+    s <- level_colsums(matrix(terms, nrow(P)))
+    J <- matrix(0, nf, nf)
+    J[uv] <- s
+    J[uv[, 2:1, drop = FALSE]] <- s
+    J
+  }
+  cv <- numeric(V)
+  Pf <- probs(cv)
+  r <- resid(Pf)
+  it <- 0L
+  while (max(abs(r)) > 1e-10 * n && it < 50L) {
+    it <- it + 1L
+    step <- level_solve(jacobian(Pf[, free, drop = FALSE]), r)
+    if (is.null(step) || !all(is.finite(step))) {
+      break
+    }
+    t <- 1
+    improved <- FALSE
+    for (h in seq_len(30)) {
+      cv_new <- cv
+      cv_new[free_col] <- cv[free_col] + t * step
+      Pf_new <- probs(cv_new)
+      r_new <- resid(Pf_new)
+      if (all(is.finite(r_new)) && level_sum(r_new^2) < level_sum(r^2)) {
+        improved <- TRUE
+        break
+      }
+      t <- t / 2
+    }
+    if (!improved) {
+      break
+    }
+    cv <- cv_new
+    Pf <- Pf_new
+    r <- r_new
+  }
+  if (!all(is.finite(cv))) {
+    return(zero)
+  }
+  if (crit$link == "softmax") {
+    cv <- cv - level_sum(cv) / length(cv)
+  }
+  structure(cv, iterations = it)
+}
+
+# The level's sums (classif_level()), in double arithmetic in a fixed order, the same on every machine. Column sums
+# of an n x k matrix: pairwise, the first half of the rows plus the second half (the odd row out kept last) until one
+# row is left.
+level_colsums <- function(M) {
+  M <- as.matrix(M)
+  if (nrow(M) == 0L) {
+    return(numeric(ncol(M)))
+  }
+  while (nrow(M) > 1L) {
+    h <- nrow(M) %/% 2L
+    top <- M[seq_len(h), , drop = FALSE] + M[h + seq_len(h), , drop = FALSE]
+    M <- if (nrow(M) %% 2L == 1L) rbind(top, M[nrow(M), , drop = FALSE]) else top
+  }
+  M[1, ]
+}
+
+# The sum of a vector (level_colsums())
+level_sum <- function(x) {
+  level_colsums(matrix(as.numeric(x), ncol = 1L))
+}
+
+# Row sums of an n x k matrix: the columns added from the first
+level_rowsums <- function(M) {
+  s <- M[, 1]
+  for (j in seq_len(ncol(M))[-1]) {
+    s <- s + M[, j]
+  }
+  s
+}
+
+# The solution of J x = r for the level's Newton step (J the k x k Jacobian, positive definite): Gaussian elimination
+# with the diagonal as pivots in their order, in double arithmetic. NULL if a pivot is not positive and finite (J
+# numerically singular), as solve() fails.
+level_solve <- function(J, r) {
+  k <- length(r)
+  A <- J
+  b <- as.numeric(r)
+  for (p in seq_len(k)) {
+    if (!is.finite(A[p, p]) || A[p, p] <= 0) {
+      return(NULL)
+    }
+    for (i in seq_len(k - p) + p) {
+      l <- A[i, p] / A[p, p]
+      for (j in seq_len(k - p) + p) {
+        A[i, j] <- A[i, j] - l * A[p, j]
+      }
+      b[i] <- b[i] - l * b[p]
+    }
+  }
+  x <- numeric(k)
+  for (i in rev(seq_len(k))) {
+    s <- b[i]
+    for (j in seq_len(k - i) + i) {
+      s <- s - A[i, j] * x[j]
+    }
+    x[i] <- s / A[i, i]
+  }
+  x
+}
+
+# Whether a fit with these params gets the classification level (classif_level()): rpf()'s option
+# randomPlantedForest.level (TRUE by default) is stored in params$level; a forest fitted before the level existed has
+# no params$level and none. Classification with loss = "logit" or "exponential".
+level_on <- function(params, mode) {
+  mode == "classification" && isTRUE(params$level) && isTRUE(params$loss %in% c("logit", "exponential"))
+}
+
+# The level of the fitted forest `fit` (the C++ object, before shift_level()) with the classification criterion
+# `crit` of its training rows: from its OOB path's last checkpoint, at its last split (the fit tracks it,
+# rpf_impl()); for a forest without out-of-bag rows (deterministic) the balance on its own predictions of the
+# training rows X. B: the families.
+classif_fit_level <- function(fit, crit, X, B, deterministic) {
+  in_sample <- function() {
+    G <- matrix(fit$predict_matrix(X, 0, 0L), ncol = crit$V)
+    n <- nrow(G)
+    classif_level(crit, seq_len(n), G, G, G, rep(1L, n), rep(0L, n), 1)
+  }
+  if (deterministic) {
+    return(in_sample())
+  }
+  path <- fit$oob_path()
+  m0 <- path$count[, 1]
+  m1 <- path$count[, 2]
+  ok <- m0 + m1 > 0
+  if (!any(ok)) {
+    return(in_sample())
+  }
+  V <- crit$V
+  cols <- (length(path$checkpoints) - 1L) * V + seq_len(V)
+  s0 <- path$sum0[ok, cols, drop = FALSE]
+  s1 <- path$sum1[ok, cols, drop = FALSE]
+  m0 <- m0[ok]
+  m1 <- m1[ok]
+  classif_level(crit, which(ok), (s0 + s1) / (m0 + m1), s0 / pmax(m0, 1), s1 / pmax(m1, 1), m0, m1, B)
+}
+
+# The level of a forest with fixed splits, after its fit (rpf_impl() tracked its OOB path at its last split):
+# classif_fit_level(), then classif_level_finish(). Returns the level.
+classif_level_fixed <- function(fit, params, crit, X) {
+  level <- classif_fit_level(fit, crit, X, params$ntrees, isTRUE(params$deterministic))
+  classif_level_finish(fit, params, level)
+  level
+}
+
+# A forest after its fit with the level: shift_level(), then the purification that rpf(purify = TRUE) asked for (the
+# C++ fit does not purify such a forest itself, as the level comes after it; rpf_impl())
+classif_level_finish <- function(fit, params, level) {
+  fit$shift_level(as.numeric(level))
+  if (isTRUE(params$purify)) {
+    fit$purify_threads(0L, 0L, 2L)
+  }
+  invisible(fit)
+}
+
 # Log-sum-exp of the rows of the score columns G (n x V), with the reference score 0 joined if `ref`
 crit_lse <- function(G, ref) {
   mx <- apply(G, 1, max)
diff --git a/R/oob_path.R b/R/oob_path.R
index d9aebf7..5934e61 100644
--- a/R/oob_path.R
+++ b/R/oob_path.R
@@ -90,7 +90,7 @@ oob_path_point <- function(s0, s1, count, y, ntrees, rows = FALSE) {
 # oob_path_point() for a classification criterion `crit` (oob_score_y()): the mean corrected and plain criterion over
 # the rows with an OOB family and their standard errors, from s0 / s1 (n x V: the sums of the families of even / odd
 # index) and count. Rows whose OOB families are all even or all odd are not corrected (as for regression). rows: as
-# for oob_path_point().
+# for oob_path_point(). With the classification level the point has the attribute `level`, the forest's level there.
 oob_path_point_crit <- function(s0, s1, count, crit, ntrees, rows = FALSE) {
   m0 <- count[, 1]
   m1 <- count[, 2]
@@ -99,18 +99,23 @@ oob_path_point_crit <- function(s0, s1, count, crit, ntrees, rows = FALSE) {
   n_ok <- sum(ok)
   s0 <- s0[ok, , drop = FALSE]
   s1 <- s1[ok, , drop = FALSE]
-  r <- crit_rows(
-    crit,
-    which(ok),
-    (s0 + s1) / m[ok],
-    s0 / pmax(m0[ok], 1),
-    s1 / pmax(m1[ok], 1),
-    m0[ok],
-    m1[ok],
-    ntrees,
-    impute = FALSE
-  )
+  G <- (s0 + s1) / m[ok]
+  G0 <- s0 / pmax(m0[ok], 1)
+  G1 <- s1 / pmax(m1[ok], 1)
+  # with the classification level (crit$level, rpf_fit()) the forest is scored at its level here (classif_level()):
+  # the forest with these splits is returned with it
+  level <- NULL
+  if (isTRUE(crit$level)) {
+    level <- classif_level(crit, which(ok), G, G0, G1, m0[ok], m1[ok], ntrees)
+    shift <- rep(as.numeric(level), each = n_ok)
+    G <- G + shift
+    G0 <- G0 + shift
+    G1 <- G1 + shift
+  }
+  r <- crit_rows(crit, which(ok), G, G0, G1, m0[ok], m1[ok], ntrees, impute = FALSE)
   point <- c(mean(r$corr), mean(r$e), stats::sd(r$corr) / sqrt(n_ok), stats::sd(r$e) / sqrt(n_ok))
+  # (the level, for the adaptive fit's forest)
+  attr(point, "level") <- level
   if (rows) list(point = point, rows = r$corr) else point
 }
 
diff --git a/R/oob_stop.R b/R/oob_stop.R
index 1ea7575..954d2aa 100644
--- a/R/oob_stop.R
+++ b/R/oob_stop.R
@@ -65,6 +65,11 @@ tune_splits_oob <- function(
     )
   }
   res <- st$result()
+  # with the classification level the forest gets its level at the chosen checkpoint, as the fit with that many splits
+  # gets it (rpf_impl())
+  if (!is.null(res$level)) {
+    classif_level_finish(fit, params, res$level)
+  }
   if (res$capped) {
     cli::cli_warn(
       c(
@@ -85,7 +90,8 @@ tune_splits_oob <- function(
 # the lowest error) to stop: once the patience or the plateau and the floor allow it, or at the last checkpoint. The
 # error is the corrected OOB mean squared error, or the mean corrected criterion. result() gives the chosen number of
 # splits, the OOB path's table of the checkpoints grown (as oob_path_table() makes it), where growing stopped, and
-# whether it stopped at the cap with the error still falling (neither the patience nor the plateau reached).
+# whether it stopped at the cap with the error still falling (neither the patience nor the plateau reached), and with
+# the classification level (a criterion with `level`) the level at the chosen checkpoint (otherwise NULL).
 #
 # The plateau (tolerance > 0, from checkpoint patience + 1 on): with b the lowest error up to c and a the lowest up to
 # c - patience (the first of ties, as oob_select()), the gain e[a] - e[b] over the last `patience` checkpoints is at
@@ -107,6 +113,8 @@ oob_stopper <- function(y, ntrees, checkpoints, patience = 4L, tolerance = 0, fl
   st$no_oob <- FALSE
   st$points <- matrix(NA_real_, C, 4)
   st$rows <- vector("list", C)
+  # (the classification level: the forest's level at each checkpoint, classif_level())
+  st$levels <- vector("list", C)
   st$step <- function(c, path) {
     if (c != st$c + 1L) {
       cli::cli_abort("The adaptive fit's checkpoints came out of order.", .internal = TRUE)
@@ -130,6 +138,9 @@ oob_stopper <- function(y, ntrees, checkpoints, patience = 4L, tolerance = 0, fl
       st$rows[[c]] <- pt$rows
       pt <- pt$point
     }
+    if (!is.null(attr(pt, "level"))) {
+      st$levels[[c]] <- attr(pt, "level")
+    }
     st$points[c, ] <- pt
     e <- st$points[seq_len(c), 1]
     best <- oob_select(e)
@@ -161,7 +172,8 @@ oob_stopper <- function(y, ntrees, checkpoints, patience = 4L, tolerance = 0, fl
       splits = as.integer(checkpoints[st$sel]),
       path = oob_path_frame(checkpoints[g], st$points[g, , drop = FALSE], st$n_rows, crit),
       stopped_at = as.integer(checkpoints[st$c]),
-      capped = st$capped
+      capped = st$capped,
+      level = st$levels[[st$sel]]
     )
   }
   st
diff --git a/R/print.R b/R/print.R
index e9e5e68..56826fd 100644
--- a/R/print.R
+++ b/R/print.R
@@ -136,7 +136,21 @@ format.rpf <- function(x, ...) {
         if (identical(params$loss, "logit") && !is.null(params$exact_step)) "exact_step",
         if (params$loss %in% c("logit", "exponential") && !is.null(params$prior_rows)) "prior_rows"
       )
-      print_params(params[loss_params])
+      # the level (classif_level()), on the score scale; multiclass per class, centred as the scores are reported
+      level <- NULL
+      if (!is.null(x$level)) {
+        lev <- x$level
+        if (length(lev) > 1) {
+          if (identical(params$loss, "logit")) {
+            lev <- c(0, lev)
+          }
+          lev <- stats::setNames(lev - mean(lev), levels(x$blueprint$ptypes$outcomes[[1]]))
+          level <- list(level = paste0(names(lev), " ", sprintf("%+.3f", lev), collapse = ", "))
+        } else {
+          level <- list(level = sprintf("%+.3f", lev))
+        }
+      }
+      print_params(c(params[loss_params], level))
     }
 
     cli::cli_text("")
diff --git a/R/prune.R b/R/prune.R
index 7267bbd..592a8c7 100644
--- a/R/prune.R
+++ b/R/prune.R
@@ -885,6 +885,10 @@ prune_refit <- function(object, predictors, y, rows, seed, nthreads) {
     mode = object$mode,
     params = params
   )
+  # (a forest with the classification level gets its own, on these rows)
+  if (level_on(params, object$mode)) {
+    classif_level_fixed(fit, params, oob_criterion(object$mode, params$loss, y[rows]), pre$predictors_matrix)
+  }
   new_rpf(
     fit = fit,
     blueprint = object$blueprint,
diff --git a/R/rpf.R b/R/rpf.R
index 644e124..90853a5 100644
--- a/R/rpf.R
+++ b/R/rpf.R
@@ -287,6 +287,29 @@
 #' three or more classes. `"exponential"` can need many more splits with many
 #' classes.
 #'
+#' ## The level of the classification losses
+#'
+#' `"logit"` and `"exponential"` average their tree families on the score
+#' scale. Where the families disagree, the probability of their mean score is
+#' not the mean of their probabilities: for a rare class (probability below
+#' 1/2, where the probability is convex in the score) it is lower, the more so
+#' the more splits the families make. So the forest would predict a rare class
+#' too rarely. After the fit the forest's score therefore gets one constant per
+#' score column, its level, which is part of its intercept: the
+#' maximum-likelihood intercept of the log-loss given the forest's out-of-bag
+#' scores, so that the mean out-of-bag probability of every class is its
+#' frequency in the training rows (corrected, as the out-of-bag criterion is,
+#' for an out-of-bag mean averaging fewer families than the forest). With
+#' `deterministic = TRUE`, which has no out-of-bag rows, the forest's own
+#' predictions of the training rows are balanced instead. The out-of-bag choice
+#' of `n_splits` scores every number of splits at its own level, and the
+#' pruning starts from the forest with its level. `$level` holds it (one value
+#' per score column) and `print()` shows it.
+#' `options(randomPlantedForest.level = FALSE)` fits without it. The setting
+#' is stored in the fit's `$params$level`, and refits (the pruning's
+#' cross-validation) follow it; a fit without the setting is refitted without
+#' a level.
+#'
 #' ## Factor predictors
 #'
 #' Factor predictors, and character predictors, which are converted to
@@ -1082,6 +1105,10 @@ rpf_bridge <- function(
     # logit's steps (exponential's are exact anyway), and the losses' pseudo-rows, by an option for studies only
     params$exact_step <- exact_step
     params$prior_rows <- rpf_prior_rows(call = call)
+    # the level of loss = "logit" and "exponential" (classif_level()), by the option randomPlantedForest.level
+    if (loss %in% c("logit", "exponential")) {
+      params$level <- rpf_level(call = call)
+    }
   }
 
   # a pruned model without its forest holds no R list of it either (export_forest)
@@ -1107,9 +1134,14 @@ rpf_bridge <- function(
 # `checkpoints` those of rpf(oob_path = ) (fixed splits) or NULL. The OOB error is the regression one or the
 # classification criterion (oob_score_y()).
 rpf_fit <- function(processed, outcomes, predictors, params, adaptive, checkpoints, export_forest, call) {
-  score_y <- if (adaptive || !is.null(checkpoints)) {
+  level <- level_on(params, outcomes$mode)
+  score_y <- if (adaptive || !is.null(checkpoints) || level) {
     oob_score_y(outcomes$mode, params$loss, processed$outcomes[[1]])
   }
+  # with the classification level every checkpoint of the OOB path is scored at its own level (oob_path_point_crit())
+  if (level) {
+    score_y$level <- TRUE
+  }
   if (adaptive) {
     chosen <- tune_splits_oob(
       outcomes$outcomes,
@@ -1129,6 +1161,9 @@ rpf_fit <- function(processed, outcomes, predictors, params, adaptive, checkpoin
       params = params,
       oob_path = checkpoints
     )
+    if (level) {
+      chosen <- list(level = classif_level_fixed(fit, params, score_y, predictors$predictors_matrix))
+    }
   }
 
   # Optionally export a compact R list representation of the forest.
@@ -1146,11 +1181,25 @@ rpf_fit <- function(processed, outcomes, predictors, params, adaptive, checkpoin
     params = params,
     forest = forest
   )
+  # the level the forest got (classif_level()), one value per score column
+  if (level) {
+    res$level <- as.numeric(chosen$level)
+  }
   if (adaptive) {
     res$oob_path <- chosen$path
     res$splits_oob <- list(n_splits = chosen$splits, stopped_at = chosen$stopped_at, capped = chosen$capped)
   } else if (!is.null(checkpoints)) {
-    res$oob_path <- oob_path_table(fit$oob_path(), score_y, params$ntrees, call = call)
+    path <- fit$oob_path()
+    # (with the level the fit also tracks its last split: only the checkpoints asked for are reported)
+    keep <- which(path$checkpoints %in% checkpoints)
+    if (length(keep) < length(path$checkpoints)) {
+      V <- ncol(path$sum0) %/% length(path$checkpoints)
+      cols <- as.integer(outer(seq_len(V), (keep - 1L) * V, "+"))
+      path$checkpoints <- path$checkpoints[keep]
+      path$sum0 <- path$sum0[, cols, drop = FALSE]
+      path$sum1 <- path$sum1[, cols, drop = FALSE]
+    }
+    res$oob_path <- oob_path_table(path, score_y, params$ntrees, call = call)
   }
   res
 }
@@ -1202,6 +1251,20 @@ rpf_exact_draw <- function(call = rlang::caller_env()) {
   exact_draw
 }
 
+# The level of loss = "logit" and "exponential" (classif_level()): the option randomPlantedForest.level (TRUE, the
+# default: on; FALSE: off). Stored in the fit's params, so that a refit (the pruning's cross-validation) and a saved
+# fit keep it; a fit without params$level gets no level (level_on()).
+rpf_level <- function(call = rlang::caller_env()) {
+  level <- getOption("randomPlantedForest.level", TRUE)
+  if (!rlang::is_bool(level)) {
+    cli::cli_abort(
+      "The option {.code randomPlantedForest.level} must be {.code TRUE} or {.code FALSE}.",
+      call = call
+    )
+  }
+  level
+}
+
 # The pseudo-rows a of the logit and exponential losses (1): the undocumented
 # option randomPlantedForest.prior_rows, for studies; NULL if it is not set.
 rpf_prior_rows <- function(call = rlang::caller_env()) {
@@ -1261,11 +1324,22 @@ rpf_param_vector <- function(params, mode) {
 # `oob_path`: NULL, or the checkpoints of the OOB path (see check_oob_path())
 # `adaptive`: NULL, or the adaptive fit's checkpoints, stopper, patience, tolerance, floor and pause_gb (see
 # tune_splits_oob())
+# With the classification level (level_on()) the forest gets its level after this (classif_level_fixed(), or for the
+# adaptive fit from the stopper, tune_splits_oob()), so the C++ fit does not purify it (classif_level_finish() does),
+# and a fit with fixed splits also tracks its OOB path at its last split, from which the level is computed.
 rpf_impl <- function(Y, X, mode = c("regression", "classification"), params, oob_path = NULL, adaptive = NULL) {
   # Final input validation, should be superfluous
   checkmate::assert_matrix(X, mode = "numeric", any.missing = FALSE)
   mode <- match.arg(mode)
-  pars <- rpf_param_vector(params, mode)
+  level <- level_on(params, mode)
+  if (level) {
+    pars <- rpf_param_vector(utils::modifyList(params, list(purify = FALSE)), mode)
+    if (is.null(adaptive) && !isTRUE(params$deterministic)) {
+      oob_path <- sort(unique(c(oob_path, as.integer(params$splits))))
+    }
+  } else {
+    pars <- rpf_param_vector(params, mode)
+  }
   # allowed feature sets as 1-based column indices (NULL: no constraint, the constructors without them)
   constraints <- if (!is.null(params$interaction_constraints)) {
     lapply(params$interaction_constraints, match, table = colnames(X))
diff --git a/man/rpf.Rd b/man/rpf.Rd
index 9fb2a7a..5aaf9e6 100644
--- a/man/rpf.Rd
+++ b/man/rpf.Rd
@@ -439,6 +439,30 @@ three or more classes. \code{"exponential"} can need many more splits with many
 classes.
 }
 
+\subsection{The level of the classification losses}{
+
+\code{"logit"} and \code{"exponential"} average their tree families on the score
+scale. Where the families disagree, the probability of their mean score is
+not the mean of their probabilities: for a rare class (probability below
+1/2, where the probability is convex in the score) it is lower, the more so
+the more splits the families make. So the forest would predict a rare class
+too rarely. After the fit the forest's score therefore gets one constant per
+score column, its level, which is part of its intercept: the
+maximum-likelihood intercept of the log-loss given the forest's out-of-bag
+scores, so that the mean out-of-bag probability of every class is its
+frequency in the training rows (corrected, as the out-of-bag criterion is,
+for an out-of-bag mean averaging fewer families than the forest). With
+\code{deterministic = TRUE}, which has no out-of-bag rows, the forest's own
+predictions of the training rows are balanced instead. The out-of-bag choice
+of \code{n_splits} scores every number of splits at its own level, and the
+pruning starts from the forest with its level. \verb{$level} holds it (one value
+per score column) and \code{print()} shows it.
+\code{options(randomPlantedForest.level = FALSE)} fits without it. The setting
+is stored in the fit's \verb{$params$level}, and refits (the pruning's
+cross-validation) follow it; a fit without the setting is refitted without
+a level.
+}
+
 \subsection{Factor predictors}{
 
 Factor predictors, and character predictors, which are converted to
diff --git a/src/include/path_tracker.hpp b/src/include/path_tracker.hpp
index a1884ec..d984ea2 100644
--- a/src/include/path_tracker.hpp
+++ b/src/include/path_tracker.hpp
@@ -37,6 +37,7 @@ struct RandomPlantedForest::PathTracker {
   std::map<std::set<int>, int32_t> dims_id; /**< the split trees' dims, numbered as they first occur */
   std::vector<std::set<int>> dims_list;
   std::vector<double> null_value; /**< the fit of the null leaf if the family makes no split (path_start()) */
+  bool has_split = false;         /**< it has made a split (path_split()) */
 };
 
 #endif // RPF_PATH_TRACKER_HPP
diff --git a/src/include/rpf.hpp b/src/include/rpf.hpp
index 08a8aa1..0e7e461 100644
--- a/src/include/rpf.hpp
+++ b/src/include/rpf.hpp
@@ -240,6 +240,10 @@ namespace rpf_cnt {
 class SubsetCounter;
 }
 
+// The value of a family's null tree, whose one leaf holds all rows: what the family adds to the intercept before its
+// other trees are purified. (lib/purify.cpp)
+std::vector<double> rpf_null_tree_value(const TreeFamily &family, size_t value_size);
+
 // Runs fn(family index) for n_families families, at most threads_to_use at a time, in batches of increasing index;
 // after_batch(start, batch) runs on the calling thread after each batch. An exception in fn is rethrown on the calling
 // thread. (lib/purify.cpp)
@@ -285,6 +289,10 @@ class RandomPlantedForest {
   );
   // Store training data without fitting and without recomputing bounds.
   void set_training_data(const Matrix2D &samples_Y, const Matrix2D &samples_X);
+  // Add c (value_size values, one per score column) to every family's null leaf: the forest's score moves by c
+  // everywhere, as part of its intercept (predict, purify, the pruning's components and get_model() see it). The
+  // classification level (R: classif_level()) is set so. An earlier purification is dropped.
+  void shift_level(const std::vector<double> &c);
   const Matrix2D &get_X() const { return X; }
   const Matrix2D &get_Y() const { return Y; }
   const std::vector<double> &get_lower_bounds() const { return lower_bounds; }
@@ -543,6 +551,8 @@ class RandomPlantedForest {
   // `k` splits are done: snapshot the predictions if k is the next checkpoint. Returns whether the family grows on
   // (always, except in an adaptive fit)
   bool path_step(PathTracker &trk, int k) const;
+  // The tracked predictions for a checkpoint: a family without a split so far at its null leaf's fit
+  std::vector<double> path_snapshot(const PathTracker &trk) const;
   // Family n is grown: its last snapshot, then add its snapshots to the sums and free its tracker
   void path_merge(size_t n);
   // The adaptive fit (oob_path.cpp; set_adaptive()): its state, while it runs and (the records) until the rebuild
diff --git a/src/lib/oob_path.cpp b/src/lib/oob_path.cpp
index e442c0b..c09ef64 100644
--- a/src/lib/oob_path.cpp
+++ b/src/lib/oob_path.cpp
@@ -54,9 +54,10 @@ RandomPlantedForest::path_start(size_t n, TreeFamily *fam, const std::vector<std
   if (path_checkpoints_.empty())
     return nullptr;
   auto t = std::make_shared<PathTracker>();
-  // the adaptive fit rebuilds the families from their records: one without a split at the chosen checkpoint gets the
-  // fit of its null leaf, as fit_null_leaf() gives it to a family that made no split (from the same rows and Y)
-  if (alone_ && Y && fam) {
+  // the fit of the null leaf, as fit_null_leaf() gives it to a family that made no split (from the same rows and Y):
+  // the path tracks a family without a split at it (path_snapshot()), and the adaptive fit's rebuild gives it to a
+  // family without a split at the chosen checkpoint
+  if (Y && fam) {
     auto null_tree = fam->find(std::set<int>{0});
     if (null_tree != fam->end() && !null_tree->second->leaves.empty() &&
         !null_tree->second->leaves[0].individuals.empty())
@@ -81,6 +82,7 @@ void RandomPlantedForest::path_split(
     PathTracker &trk, const Split &split, const std::vector<double> *inc_s, const std::vector<double> *inc_b
 ) const
 {
+  trk.has_split = true;
   const std::vector<double> &Ms = inc_s ? *inc_s : split.M_s, &Mb = inc_b ? *inc_b : split.M_b;
   const DecisionTree &T = *split.tree_index;
   auto &lists = trk.in_leaf[T.split_dims];
@@ -127,12 +129,28 @@ void RandomPlantedForest::path_split(
   }
 }
 
+// The family's predictions at its OOB rows now, for a checkpoint. A family that has made no split predicts its null
+// leaf's fit (fit_null_leaf(): the fit of its rows, e.g. their mean, or logit's log-odds of its sample), not the 0
+// its rows start from: the path holds what a forest with so many splits predicts. A family with a split: the tracked
+// predictions (its null leaf stays at 0).
+std::vector<double> RandomPlantedForest::path_snapshot(const PathTracker &trk) const
+{
+  if (trk.has_split || trk.null_value.size() != value_size)
+    return trk.pred;
+  std::vector<double> out(trk.pred.size());
+  for (size_t p = 0; p < trk.rows.size(); ++p)
+    for (size_t v = 0; v < value_size; ++v)
+      out[p * value_size + v] = trk.null_value[v];
+  return out;
+}
+
 bool RandomPlantedForest::path_step(PathTracker &trk, int k) const
 {
   if (alone_)
     return alone_step(trk, (size_t)k);
   if (trk.next < path_checkpoints_.size() && path_checkpoints_[trk.next] == k) {
-    trk.snaps.insert(trk.snaps.end(), trk.pred.begin(), trk.pred.end());
+    const std::vector<double> now = path_snapshot(trk);
+    trk.snaps.insert(trk.snaps.end(), now.begin(), now.end());
     ++trk.next;
   }
   return true;
@@ -570,7 +588,7 @@ bool RandomPlantedForest::alone_step(PathTracker &trk, size_t k) const
     if (trk.next >= C || (size_t)path_checkpoints_[trk.next] != k)
       return !A.cancel;
     const size_t c = trk.next++;
-    trk.csnaps.push_back(trk.pred);
+    trk.csnaps.push_back(path_snapshot(trk));
     trk.rec_at.push_back((uint32_t)trk.rec.size());
     if (A.cancel)
       return false;
diff --git a/src/lib/prune_columns.cpp b/src/lib/prune_columns.cpp
index 105c0db..e370899 100644
--- a/src/lib/prune_columns.cpp
+++ b/src/lib/prune_columns.cpp
@@ -73,7 +73,8 @@ void RandomPlantedForest::family_columns(
     bounds.erase(std::unique(bounds.begin(), bounds.end()), bounds.end());
     lim_list[d - 1] = bounds;
   }
-  intercept.assign(VS, 0.0);
+  // the null tree's value (a family without a split), as purify_fast_exact_family()
+  intercept = rpf_null_tree_value(curr_family, value_size);
   std::map<std::set<int>, double, utils::setComp> tot_cache;
   auto get_tot_for_U = [&](const std::set<int> &U) -> double {
     auto it = tot_cache.find(U);
diff --git a/src/lib/purify.cpp b/src/lib/purify.cpp
index 43f8b44..f10e9d3 100644
--- a/src/lib/purify.cpp
+++ b/src/lib/purify.cpp
@@ -13,7 +13,7 @@
 
 // The value of the family's null tree: 0, unless the family made no split, when its one leaf holds the family's
 // fit (fit_null_leaf()). The intercept starts from it.
-static std::vector<double> null_tree_value(const TreeFamily &family, size_t value_size)
+std::vector<double> rpf_null_tree_value(const TreeFamily &family, size_t value_size)
 {
   std::vector<double> value(value_size, 0);
   auto null_tree = family.find(std::set<int>{0});
@@ -68,7 +68,7 @@ void RandomPlantedForest::purify_fast_exact_family(
     lim_list[d - 1] = bounds;
   }
 
-  acc.intercept = null_tree_value(curr_family, value_size);
+  acc.intercept = rpf_null_tree_value(curr_family, value_size);
 
   // Precompute tot(U) with half-open domain [front, back)
   std::map<std::set<int>, double, utils::setComp> tot_cache;
@@ -874,7 +874,8 @@ void RandomPlantedForest::purify_3_family(TreeFamily &curr_family, int maxp_inte
   // ------------- create new trees -------------
   grids.insert(grids.begin(), grid::NDGrid());
   values.insert(
-      values.begin(), utils::Matrix<std::vector<double>>(std::vector<int>{1}, null_tree_value(curr_family, value_size))
+      values.begin(),
+      utils::Matrix<std::vector<double>>(std::vector<int>{1}, rpf_null_tree_value(curr_family, value_size))
   );
   values_old.insert(
       values_old.begin(), utils::Matrix<std::vector<double>>(std::vector<int>{1}, std::vector<double>(value_size, 0))
diff --git a/src/lib/rpf.cpp b/src/lib/rpf.cpp
index 5b2c0fd..2d887a5 100644
--- a/src/lib/rpf.cpp
+++ b/src/lib/rpf.cpp
@@ -173,6 +173,29 @@ void RandomPlantedForest::set_training_data(const Matrix2D &samples_Y, const Mat
   drop_purify_counts(); // (counts of the earlier rows)
 }
 
+void RandomPlantedForest::shift_level(const std::vector<double> &c)
+{
+  if (c.size() != value_size)
+    throw std::invalid_argument("The level needs one value per score column.");
+  for (double v : c)
+    if (!std::isfinite(v))
+      throw std::invalid_argument("The level must be finite.");
+  for (auto &family : tree_families) {
+    auto null_tree = family.find(std::set<int>{0});
+    if (null_tree == family.end() || null_tree->second->leaves.empty())
+      throw std::logic_error("A tree family has no null leaf.");
+    std::vector<double> &value = null_tree->second->leaves[0].value;
+    if (value.size() != value_size)
+      throw std::logic_error("A null leaf has the wrong number of values.");
+    for (size_t k = 0; k < value_size; ++k)
+      value[k] += c[k];
+  }
+  purified = false;
+  purify_capped = false;
+  rect_forest_.reset();
+  purify_cap = 0;
+}
+
 void RandomPlantedForest::set_model(const std::vector<FamilyExport> &model)
 {
   tree_families = std::vector<TreeFamily>(model.size());
diff --git a/src/randomPlantedForest.cpp b/src/randomPlantedForest.cpp
index 302dcea..9d68885 100644
--- a/src/randomPlantedForest.cpp
+++ b/src/randomPlantedForest.cpp
@@ -7,6 +7,7 @@ template <class T> void expose_methods(Rcpp::class_<T> &cls)
       .method("fit", &T::fit)
       .method("set_shape", &T::set_shape)
       .method("set_training_data", &T::set_training_data)
+      .method("shift_level", &T::shift_level)
       .method("get_data", &T::get_data)
       .method("get_bounds", &T::get_bounds)
       .method("get_shape", &T::get_shape)
diff --git a/src/rcpp_interface.h b/src/rcpp_interface.h
index 75dc4d9..3136a93 100644
--- a/src/rcpp_interface.h
+++ b/src/rcpp_interface.h
@@ -193,6 +193,7 @@ template <class Core> class RcppForest {
   {
     core.set_training_data(to_matrix2d(samples_Y), to_matrix2d(samples_X));
   }
+  void shift_level(const Rcpp::NumericVector &c) { core.shift_level(Rcpp::as<std::vector<double>>(c)); }
   Rcpp::List get_data()
   {
     return Rcpp::List::create(
diff --git a/tests/testthat/_snaps/print.md b/tests/testthat/_snaps/print.md
index 3545798..be2bc1e 100644
--- a/tests/testthat/_snaps/print.md
+++ b/tests/testthat/_snaps/print.md
@@ -60,6 +60,7 @@
       -- Loss 
               loss: logit
         exact_step: TRUE
+             level: setosa -0.279, versicolor +0.086, virginica +0.193
       
       i Fit using 1 thread, also the default for `predict()` and `purify()`.
 
diff --git a/tests/testthat/test-classif-level.R b/tests/testthat/test-classif-level.R
new file mode 100644
index 0000000..cfb40c3
--- /dev/null
+++ b/tests/testthat/test-classif-level.R
@@ -0,0 +1,313 @@
+# The classification level (classif_level()): one constant per score column, the maximum-likelihood intercept of the
+# log-loss given the forest's out-of-bag scores (the balance property), in every family's null leaf. The option
+# randomPlantedForest.level (TRUE by default) switches it; off, a fit is the fit of earlier versions.
+
+level_data <- function(n = 400, seed = 3, K = 2) {
+  set.seed(seed)
+  x <- data.frame(x1 = runif(n), x2 = runif(n), x3 = runif(n))
+  if (K == 2) {
+    p <- plogis(-2.5 + 2 * x$x1 + sin(2 * pi * x$x2))
+    x$y <- factor(ifelse(runif(n) < p, "yes", "no"), levels = c("no", "yes"))
+  } else {
+    eta <- cbind(0, 2 * x$x1 - 1, -2.5 + 2 * x$x2)
+    cls <- apply(eta, 1, function(e) sample.int(3, 1, prob = exp(e)))
+    x$y <- factor(c("a", "b", "c")[cls])
+  }
+  x
+}
+
+level_raw <- function(fit, d) {
+  X <- preprocess_predictors_predict(fit, hardhat::forge(d, fit$blueprint)$predictors)
+  matrix(fit$fit$predict_matrix(X, 0, 0L), nrow = nrow(d))
+}
+
+with_level <- function(on, code) {
+  old <- options(randomPlantedForest.level = on)
+  on.exit(options(old))
+  code
+}
+
+test_that("classif_level() solves the balance equations (binary, multiclass logit, SAMME)", {
+  set.seed(11)
+  n <- 500
+  for (case in list(c("logit", 2), c("logit", 4), c("exponential", 2), c("exponential", 3))) {
+    K <- as.integer(case[2])
+    y <- factor(sample(letters[1:K], n, TRUE, prob = c(0.6, rep(0.4 / (K - 1), K - 1))), levels = letters[1:K])
+    crit <- oob_criterion("classification", case[1], y)
+    G <- matrix(rnorm(n * crit$V, sd = 0.7), n)
+    # no correction (B = m): the mean probability of every class is its frequency
+    lev <- classif_level(crit, seq_len(n), G, G, G, rep(1L, n), rep(1L, n), 2)
+    expect_length(lev, crit$V)
+    expect_lte(attr(lev, "iterations"), 10L)
+    P <- crit_logloss(crit, G + rep(lev, each = n), crit$cls)$P
+    Pf <- if (crit$link == "sigmoid") {
+      cbind(1 - P[, 1], P[, 1])
+    } else if (crit$link == "softmax_ref") {
+      cbind(1 - rowSums(P), P)
+    } else {
+      P
+    }
+    expect_equal(colMeans(Pf), as.numeric(table(y)) / n, tolerance = 1e-10)
+    if (crit$link == "softmax") {
+      expect_equal(sum(lev), 0, tolerance = 1e-12)
+    }
+    # the level is the maximum-likelihood intercept: the log-loss is lowest there
+    ll <- function(c) sum(crit_logloss(crit, G + rep(c, each = n), crit$cls)$e)
+    for (v in seq_len(crit$V)) {
+      e <- replace(numeric(crit$V), v, 1e-3)
+      expect_lt(ll(lev), ll(lev + e))
+      expect_lt(ll(lev), ll(lev - e))
+    }
+  }
+  # a class without an out-of-bag row: no finite level, 0
+  y <- factor(c(rep("a", 10), "b"), levels = c("a", "b", "c"))
+  crit <- oob_criterion("classification", "logit", y)
+  G <- matrix(0, 10, 2)
+  expect_identical(as.numeric(classif_level(crit, 1:10, G, G, G, rep(1L, 10), rep(1L, 10), 2)), c(0, 0))
+})
+
+test_that("the level moves the forest's score by a constant; off, no level and the same splits", {
+  d <- level_data()
+  fits <- lapply(c(FALSE, TRUE), function(on) {
+    with_level(on, {
+      set.seed(4)
+      rpf(y ~ ., data = d, loss = "logit", max_interaction = 2, n_families = 20, n_splits = 30, prune = FALSE)
+    })
+  })
+  expect_null(fits[[1]]$level)
+  expect_false(fits[[1]]$params$level)
+  expect_true(fits[[2]]$params$level)
+  lev <- fits[[2]]$level
+  expect_length(lev, 1)
+  # a rare class: the forest's mean probability is below the frequency, the level raises it
+  expect_gt(lev, 0)
+  diff <- level_raw(fits[[2]], d) - level_raw(fits[[1]], d)
+  expect_equal(diff, matrix(lev, nrow(d), 1), tolerance = 1e-12)
+  # the same leaves but the null leaves
+  m0 <- fits[[1]]$fit$get_model()
+  m1 <- fits[[2]]$fit$get_model()
+  for (b in seq_along(m0)) {
+    null <- which(vapply(m0[[b]]$variables, function(v) identical(as.numeric(v), 0), TRUE))
+    expect_identical(m0[[b]]$values[-null], m1[[b]]$values[-null])
+    expect_equal(unlist(m1[[b]]$values[null]), unlist(m0[[b]]$values[null]) + lev, tolerance = 1e-15)
+  }
+  # the out-of-bag balance at the forest's last split: the OOB mean probability is the frequency up to the
+  # second-order term of the OOB correction, which is large with 20 families (1 / m - 1 / B about 0.08)
+  path <- fits[[2]]$fit$oob_path()
+  m0 <- path$count[, 1]
+  m1 <- path$count[, 2]
+  ok <- m0 + m1 > 0
+  last <- length(path$checkpoints)
+  s0 <- path$sum0[ok, last, drop = FALSE]
+  s1 <- path$sum1[ok, last, drop = FALSE]
+  g <- (s0 + s1) / (m0 + m1)[ok]
+  expect_equal(mean(plogis(g + lev)), mean(d$y[ok] == "yes"), tolerance = 0.05)
+  crit <- oob_criterion("classification", "logit", d$y)
+  expect_equal(
+    as.numeric(classif_level(crit, which(ok), g, s0 / pmax(m0[ok], 1), s1 / pmax(m1[ok], 1), m0[ok], m1[ok], 20)),
+    lev
+  )
+  # without the second-order term (B = m: as many families as the OOB means), the plain balance
+  lev0 <- classif_level(crit, which(ok), g, g, g, rep(1L, sum(ok)), rep(1L, sum(ok)), 2)
+  expect_equal(mean(plogis(g + as.numeric(lev0))), mean(d$y[ok] == "yes"), tolerance = 1e-10)
+  expect_match(paste(format(fits[[2]]), collapse = "\n"), "level: +", fixed = TRUE)
+  expect_no_match(paste(format(fits[[1]]), collapse = "\n"), "level:", fixed = TRUE)
+})
+
+test_that("multiclass logit and exponential get one level per score column", {
+  d <- level_data(K = 3)
+  for (loss in c("logit", "exponential")) {
+    set.seed(5)
+    fit <- rpf(y ~ ., data = d, loss = loss, max_interaction = 2, n_families = 20, n_splits = 25, prune = FALSE)
+    expect_length(fit$level, if (loss == "logit") 2L else 3L)
+    if (loss == "exponential") {
+      expect_equal(sum(fit$level), 0, tolerance = 1e-12)
+    }
+    with_level(FALSE, {
+      set.seed(5)
+      fit0 <- rpf(y ~ ., data = d, loss = loss, max_interaction = 2, n_families = 20, n_splits = 25, prune = FALSE)
+    })
+    expect_equal(
+      level_raw(fit, d) - level_raw(fit0, d),
+      matrix(fit$level, nrow(d), length(fit$level), byrow = TRUE),
+      tolerance = 1e-12
+    )
+  }
+})
+
+test_that("L1 / L2 classification and regression get no level", {
+  d <- level_data()
+  set.seed(6)
+  fit <- rpf(y ~ ., data = d, loss = "L2", n_families = 10, n_splits = 10, prune = FALSE)
+  expect_null(fit$level)
+  expect_null(fit$params$level)
+  d$z <- as.numeric(d$y == "yes") + d$x1
+  set.seed(6)
+  fit <- rpf(z ~ x1 + x2 + x3, data = d, n_families = 10, n_splits = 10, prune = FALSE)
+  expect_null(fit$level)
+})
+
+test_that("the adaptive fit scores every checkpoint at its level and returns the forest at the chosen one's", {
+  d <- level_data(n = 300)
+  set.seed(7)
+  fit <- rpf(y ~ ., data = d, loss = "logit", max_interaction = 2, n_families = 20, prune = FALSE)
+  # the chosen checkpoint's level, recomputed from the fit's path at that checkpoint (a fit with the chosen splits)
+  set.seed(7)
+  fixed <- rpf(
+    y ~ .,
+    data = d,
+    loss = "logit",
+    max_interaction = 2,
+    n_families = 20,
+    prune = FALSE,
+    n_splits = fit$params$splits
+  )
+  expect_equal(fit$level, fixed$level, tolerance = 1e-12)
+  expect_equal(level_raw(fit, d), level_raw(fixed, d), tolerance = 1e-12)
+  # the path of a fit with fixed splits scores each checkpoint at its own level; the last one is the adaptive path's
+  set.seed(7)
+  path_fit <- rpf(
+    y ~ .,
+    data = d,
+    loss = "logit",
+    max_interaction = 2,
+    n_families = 20,
+    prune = FALSE,
+    n_splits = fit$params$splits,
+    oob_path = TRUE
+  )
+  k <- match(fit$params$splits, fit$oob_path$n_splits)
+  expect_equal(tail(path_fit$oob_path$oob_loss, 1), fit$oob_path$oob_loss[k], tolerance = 1e-12)
+})
+
+test_that("a path of fewer checkpoints than the splits reports those only; the level is the last split's", {
+  d <- level_data(n = 300)
+  set.seed(8)
+  fit <- rpf(
+    y ~ .,
+    data = d,
+    loss = "logit",
+    max_interaction = 2,
+    n_families = 20,
+    n_splits = 30,
+    prune = FALSE,
+    oob_path = c(5, 10)
+  )
+  expect_identical(fit$oob_path$n_splits, c(5L, 10L))
+  set.seed(8)
+  fit30 <- rpf(y ~ ., data = d, loss = "logit", max_interaction = 2, n_families = 20, n_splits = 30, prune = FALSE)
+  expect_identical(fit$level, fit30$level)
+})
+
+test_that("purify = TRUE purifies the forest with its level; the pruned model's intercept holds it", {
+  d <- level_data()
+  set.seed(9)
+  fit <- rpf(
+    y ~ .,
+    data = d,
+    loss = "logit",
+    max_interaction = 2,
+    n_families = 20,
+    n_splits = 30,
+    prune = FALSE,
+    purify = TRUE
+  )
+  expect_true(fit$fit$is_purified())
+  set.seed(9)
+  raw <- rpf(y ~ ., data = d, loss = "logit", max_interaction = 2, n_families = 20, n_splits = 30, prune = FALSE)
+  expect_equal(predict(fit, d, type = "link"), predict(raw, d, type = "link"), tolerance = 1e-12)
+  # rule D at tau = 0 keeps every component and the intercept: the forest
+  pr <- prune_components(raw, d, rule = "D", tau = 0, select = "min")
+  expect_equal(predict(pr, d, type = "prob"), predict(raw, d, type = "prob"), tolerance = 1e-10)
+})
+
+test_that("deterministic: the level balances the forest's own predictions", {
+  d <- level_data()
+  fit <- rpf(y ~ ., data = d, loss = "logit", max_interaction = 2, n_families = 1, n_splits = 15, deterministic = TRUE)
+  expect_length(fit$level, 1)
+  expect_equal(mean(predict(fit, d, type = "prob")$.pred_yes), mean(d$y == "yes"), tolerance = 1e-10)
+})
+
+test_that("the pruning's cross-validation refits each fold with its level; a fit without params$level gets none", {
+  d <- level_data(n = 300)
+  set.seed(10)
+  fit <- rpf(y ~ ., data = d, loss = "logit", max_interaction = 2, n_families = 10, n_splits = 20, prune = FALSE)
+  refit <- prune_refit(fit, d[c("x1", "x2", "x3")], d$y, seq_len(200), seed = 1, nthreads = 1L)
+  # its level: the one a fit on those rows gets
+  set.seed(1)
+  direct <- rpf(
+    y ~ .,
+    data = d[seq_len(200), ],
+    loss = "logit",
+    max_interaction = 2,
+    n_families = 10,
+    n_splits = 20,
+    prune = FALSE
+  )
+  expect_equal(level_raw(refit, d), level_raw(direct, d), tolerance = 1e-12)
+  # a fit saved before the level existed has no params$level: its refits get none
+  old <- fit
+  old$params$level <- NULL
+  expect_false(level_on(old$params, "classification"))
+  refit_old <- prune_refit(old, d[c("x1", "x2", "x3")], d$y, seq_len(200), seed = 1, nthreads = 1L)
+  with_level(FALSE, {
+    set.seed(1)
+    direct0 <- rpf(
+      y ~ .,
+      data = d[seq_len(200), ],
+      loss = "logit",
+      max_interaction = 2,
+      n_families = 10,
+      n_splits = 20,
+      prune = FALSE
+    )
+  })
+  expect_identical(level_raw(refit_old, d), level_raw(direct0, d))
+})
+
+test_that("the option must be TRUE or FALSE", {
+  d <- level_data(n = 100)
+  with_level("yes", {
+    expect_error(
+      rpf(y ~ ., data = d, loss = "logit", n_families = 2, n_splits = 2, prune = FALSE),
+      "randomPlantedForest.level"
+    )
+  })
+})
+
+test_that("families without a split: the level is computed from their fits, and is small", {
+  d <- level_data(n = 100)
+  # 100 rows: no threshold keeps 51 rows on each side, so no family splits
+  fit_path <- function(level) {
+    with_level(level, {
+      set.seed(2)
+      rpf(
+        y ~ .,
+        data = d,
+        loss = "logit",
+        n_families = 10,
+        n_splits = 5,
+        min_leaf_rows = 51,
+        prune = FALSE,
+        oob_path = TRUE
+      )
+    })
+  }
+  on <- fit_path(TRUE)
+  expect_true(all(vapply(on$fit$get_model(), function(f) length(f$variables) == 1, TRUE)))
+  path <- on$fit$oob_path()
+  m <- path$count[, 1] + path$count[, 2]
+  ok <- m > 0
+  g <- (path$sum0[ok, 1] + path$sum1[ok, 1]) / m[ok]
+  # the OOB path holds the families' fits (test-no-split-family.R), without the level
+  nulls <- vapply(on$fit$get_model(), function(f) f$values[[1]][[1]][1], 0) - on$level
+  expect_true(all(g >= min(nulls) - 1e-12 & g <= max(nulls) + 1e-12))
+  expect_lt(max(abs(g - stats::qlogis(mean(d$y == "yes")))), 0.5)
+  # the OOB fits of the families are balanced already, up to the bootstrap: a small level
+  expect_lt(abs(on$level), 0.2)
+  # off: the same path and the same families, without the level
+  off <- fit_path(FALSE)
+  expect_identical(off$fit$oob_path(), path)
+  nulls_off <- vapply(off$fit$get_model(), function(f) f$values[[1]][[1]][1], 0)
+  expect_equal(nulls_off, nulls, tolerance = 1e-12)
+})
diff --git a/tests/testthat/test-classification-pipeline.R b/tests/testthat/test-classification-pipeline.R
index 0227b03..b3b65e9 100644
--- a/tests/testthat/test-classification-pipeline.R
+++ b/tests/testthat/test-classification-pipeline.R
@@ -390,6 +390,8 @@ test_that("the log-loss Hessian is the numerical one; logit's centred shares do
 })
 
 test_that("n_splits = NULL does not over-grow on a flat out-of-bag curve (three-class logit)", {
+  # the curve described above is the one without the classification level (a levelled forest has its own)
+  rlang::local_options(randomPlantedForest.level = FALSE)
   skip_on_cran()
   # A family pauses at most one doubling past the first family's estimate and grows on when the stop comes later, so
   # on a long, flat out-of-bag curve each grows to about the stop (before that fix, the estimate from part of the
@@ -400,7 +402,7 @@ test_that("n_splits = NULL does not over-grow on a flat out-of-bag curve (three-
   # curves: there the bounds are checked at their stops
   d <- cp_data(300)
   old <- options(randomPlantedForest.exact_draw = TRUE, randomPlantedForest.adaptive_tolerance = 0)
-  on.exit(options(old))
+  on.exit(options(old), add = TRUE)
   set.seed(8)
   f <- rpf(
     d$x,
@@ -454,6 +456,8 @@ test_that("n_splits = NULL does not over-grow on a flat out-of-bag curve (three-
 })
 
 test_that("the OOB path of a classification forest: prefix property, n_splits = NULL = the fixed fit", {
+  # the path's sums are the families' predictions before the level, which a fit adds after its last split
+  rlang::local_options(randomPlantedForest.level = FALSE)
   skip_on_cran()
   d <- cp_data(300)
   tot <- function(f) {
diff --git a/tests/testthat/test-driftfree-losses.R b/tests/testthat/test-driftfree-losses.R
index 56f1d18..d11b49f 100644
--- a/tests/testthat/test-driftfree-losses.R
+++ b/tests/testthat/test-driftfree-losses.R
@@ -180,6 +180,8 @@ test_that("logit, 3 classes: the Newton step is the full K x K one; exact_step r
 })
 
 test_that("one split from the scores 0: the leaves get the exact minimisers, or one Newton step", {
+  # the leaves' steps, before the level (test-classif-level.R)
+  rlang::local_options(randomPlantedForest.level = FALSE)
   set.seed(4)
   n <- 200
   x <- data.frame(x1 = runif(n))
@@ -222,6 +224,8 @@ test_that("one split from the scores 0: the leaves get the exact minimisers, or
 })
 
 test_that("separable classes: one split each, then no candidate gains and the fit stays as it is", {
+  # the leaves' steps, before the level (test-classif-level.R)
+  rlang::local_options(randomPlantedForest.level = FALSE)
   n <- 120
   x <- data.frame(x1 = seq(0.005, 0.995, length.out = n))
   y <- factor(ifelse(x$x1 < 0.3 | x$x1 >= 0.7, "a", "b"))
@@ -289,6 +293,8 @@ test_that("the logit forest does not depend on which class is the reference", {
 })
 
 test_that("the running sums choose the splits of the per-candidate pass, in every split structure", {
+  # a refit through the C++ object (set_parameters()) gets no level: compared without it
+  rlang::local_options(randomPlantedForest.level = FALSE)
   # The searches score a split point from the sums of the rows' terms along the leaf's sorted rows (or per bin, or by
   # one pass per try); set_parameters("scan", 0) (internal) scores every try by a pass over its leaf's rows instead.
   # The leaves' steps come from the winner's rows either way, so on data without ties (no bootstrap: deterministic;
@@ -326,6 +332,8 @@ test_that("the running sums choose the splits of the per-candidate pass, in ever
 })
 
 test_that("logit keeps its rows' terms or computes them: the same fit", {
+  # a refit through the C++ object (set_parameters()) gets no level: compared without it
+  rlang::local_options(randomPlantedForest.level = FALSE)
   set.seed(8)
   n <- 200
   x <- data.frame(x1 = runif(n), x2 = runif(n))
diff --git a/tests/testthat/test-exponential-samme.R b/tests/testthat/test-exponential-samme.R
index 94c9b0f..4f47d73 100644
--- a/tests/testthat/test-exponential-samme.R
+++ b/tests/testthat/test-exponential-samme.R
@@ -80,6 +80,8 @@ test_that("binary exponential: the weight is exp(-y g / 2) of the predicted scor
 })
 
 test_that("the exponential's first leaf step is log(n_k + a / K) - mean_j log(n_j + a / K), a = 1 pseudo-row", {
+  # the leaves' steps, before the level (test-classif-level.R)
+  rlang::local_options(randomPlantedForest.level = FALSE)
   d <- samme_data(3)
   set.seed(1)
   m <- rpf(
diff --git a/tests/testthat/test-fit-digests.R b/tests/testthat/test-fit-digests.R
index 97b00d3..0227de4 100644
--- a/tests/testthat/test-fit-digests.R
+++ b/tests/testthat/test-fit-digests.R
@@ -7,11 +7,14 @@
 # exponential fits (`loss3_exponential`, `pur_cls3_*`: SAMME, the multiclass exponential loss), the three-class
 # L2 fit (`loss3_L2`: L2 classification scores a split with regression's arithmetic, which breaks a tie otherwise), and
 # every logit and exponential fit (`cls_*`, `loss*_logit`, `loss*_exponential`, `pur_cls*` but their names: the
-# losses' penalised step without clamps). A change meant to leave every fit as it was keeps them all, and a change
-# that moves fits re-pins exactly those.
+# losses' penalised step without clamps), and again every logit and exponential fit (the classification level,
+# test-classif-level.R; with options(randomPlantedForest.level = FALSE) they are master's). A change meant to leave
+# every fit as it was keeps them all, and a change that moves fits re-pins exactly those.
 #
 # The digests depend on the platform's floating point and random integers (libm's exp and log draw the candidates),
-# so they are checked only where they were computed: macOS on arm64.
+# so they are checked only where they were computed: macOS on arm64. They do not depend on the BLAS / LAPACK that R
+# uses (these fits call none; the level is computed without them, classif_level()), so R's reference BLAS (CRAN's
+# macOS R), Accelerate and OpenBLAS give the same.
 
 test_that("small fits are bitwise as pinned", {
   skip_on_cran()
@@ -23,31 +26,31 @@ test_that("small fits are bitwise as pinned", {
     reg_leaves_det0_nt1 = "4d69f27123f4bf5b78bcb9228a33ae42",
     reg_leaves_det0_nt2 = "4d69f27123f4bf5b78bcb9228a33ae42",
     reg_leaves_det1_nt1 = "eb4a7b16933b5f1c88b2a6df1a8d5f63",
-    cls_leaves = "ea00222afc71241ba4cd8d68b1386b4e",
+    cls_leaves = "bc74cd1927cf1a839d4ce242455bb3a9",
     reg_cur_trees_1_det0_nt1 = "f6bb739f4bc8010a420cda3570ca266d",
     reg_cur_trees_1_det0_nt2 = "f6bb739f4bc8010a420cda3570ca266d",
     reg_cur_trees_1_det1_nt1 = "f61d46f5d76d988905d0992f11b0be5c",
-    cls_cur_trees_1 = "8342865d55b2727867ed04041ae31cc4",
+    cls_cur_trees_1 = "11719e5c548a12f4d1631048bac4a583",
     reg_cur_trees_2_det0_nt1 = "40d7845817eec7abf4547b51c2b25af7",
     reg_cur_trees_2_det0_nt2 = "40d7845817eec7abf4547b51c2b25af7",
     reg_cur_trees_2_det1_nt1 = "153661c73977f6b9a979e066f57fa3b0",
-    cls_cur_trees_2 = "41dd6f683e1cc359fb7120de1bb604f6",
+    cls_cur_trees_2 = "67c9e6f8d931b4da212e205f55272d35",
     reg_res_trees_det0_nt1 = "df541a66a06983f90c60bc21fe4cc317",
     reg_res_trees_det0_nt2 = "df541a66a06983f90c60bc21fe4cc317",
     reg_res_trees_det1_nt1 = "f3c531da7d664a9a8bff4ce0de0e1e9a",
-    cls_res_trees = "cd41fe2592fa698efb964c42e0585ccd",
+    cls_res_trees = "0236fb00aedd61d08c4b41b74cd253a7",
     reg_hist_det0_nt1 = "bda13fd7237b1bf2bcb8358e75ad659e",
     reg_hist_det0_nt2 = "bda13fd7237b1bf2bcb8358e75ad659e",
     reg_hist_det1_nt1 = "42562c42ec0b77ec51ab7435eb1aefc8",
-    cls_hist = "340c32f55abb8a82b31e86f60ce16b6b",
+    cls_hist = "719991df65460831fa19ff46ac215275",
     loss2_L1 = "dba097acec118e611d815c3a827a9584",
     loss2_L2 = "85225648f3adc80ba30457e17b58bdcb",
-    loss2_logit = "8c8fbbb8303640dc4c427d88cd5c6485",
-    loss2_exponential = "8263e3e7d9a5db55c64b9b3c4e5f7f92",
+    loss2_logit = "c56a5e11d60f2288d70c6586237aa3b0",
+    loss2_exponential = "ebaa3818d1b9abba9b33c91b39faa9ce",
     loss3_L1 = "843dbb7d5d6bafd8eeea3bf3c67487d2",
     loss3_L2 = "c536e887470b870fe87d9304f9f55493",
-    loss3_logit = "83019105a5aa663469cebfba59029d89",
-    loss3_exponential = "b1b1c73482efca7888f794065fee677b",
+    loss3_logit = "f6da6785a4c565d8371f8bfa9398fc3d",
+    loss3_exponential = "12266abb766d52f722abbbc3bbeb7e59",
     pur_reg_components = "8bf9507e8cb6ea716ec1315c58a5de81",
     pur_reg_names = "42622b48bfaa4d8b2d366e09a55c3899",
     pur_reg_predict = "e6f84aa30eedf710528f06e78974a9ea",
@@ -60,12 +63,12 @@ test_that("small fits are bitwise as pinned", {
     pur_reg_cap2_components = "d45cc7f004f85b3d25ece7b214ee21d4",
     pur_reg_cap2_names = "42622b48bfaa4d8b2d366e09a55c3899",
     pur_reg_cap2_predict = "0cc253088d2a69994515aeead292a80c",
-    pur_cls_components = "0cf62d6b6e8cadfae3d42cba3e5fb765",
+    pur_cls_components = "2dcb042bbbd344fabc8473fd487bbdf5",
     pur_cls_names = "794adb1c06fff2ed25b9afca6d9c654c",
-    pur_cls_predict = "14baa1eda3d55e5a1563d44b1537c377",
-    pur_cls3_components = "633193af37d38ee9f17b0d5df2cb3c23",
+    pur_cls_predict = "c8429429261b27535c675d9a4eb377ce",
+    pur_cls3_components = "1a1bfe68d5da0fdac396ce5e54c3cfa3",
     pur_cls3_names = "28a2f21ce2d5f51e5a2bb15a2eee3069",
-    pur_cls3_predict = "fef1a31c437d2efebbb14fdf3afc5866"
+    pur_cls3_predict = "64d1b3fb84cfffba12f77636c788645c"
   )
   got <- port_fit_digests()
   expect_identical(names(got), names(expected))
diff --git a/tests/testthat/test-no-split-family.R b/tests/testthat/test-no-split-family.R
index 7d145ba..66bec45 100644
--- a/tests/testthat/test-no-split-family.R
+++ b/tests/testthat/test-no-split-family.R
@@ -30,6 +30,8 @@ test_that("a family without a split predicts the mean of its rows", {
 })
 
 test_that("a classification family without a split predicts its rows' fit", {
+  # the null leaf's own fit, before the level (test-classif-level.R)
+  rlang::local_options(randomPlantedForest.level = FALSE)
   dat <- no_split_data()
   dat[c("x1", "x2")] <- 1
   shares <- function(y) as.numeric(prop.table(table(y)))
@@ -128,6 +130,8 @@ test_that("min_leaf_rows above half the rows leaves every family without a split
 })
 
 test_that("rpf()'s default call (splits chosen out of bag): a family without a split predicts its rows' fit", {
+  # the null leaf's own fit, before the level (test-classif-level.R)
+  rlang::local_options(randomPlantedForest.level = FALSE)
   # The adaptive fit rebuilds the forest from the families' records at the chosen number of splits; a family without
   # a split there gets its rows' fit, as a fit with that number of splits gives it. Checked on the forest itself
   # (`$fit`, kept with keep_forest = TRUE) and on the pruned model.
@@ -168,3 +172,72 @@ test_that("rpf()'s default call (splits chosen out of bag): a family without a s
     check_cls(fit, dat[[outcome]], paste(outcome, "min_leaf_rows"))
   }
 })
+
+test_that("the pruned model's intercept includes the families that made no split", {
+  # one row differs in the only predictor: the families whose bootstrap sample misses it cannot split
+  set.seed(1)
+  n <- 60
+  d <- data.frame(x1 = c(1, rep(0, n - 1)), y = rnorm(n, 5))
+  fit <- rpf(y ~ x1, data = d, n_families = 20, n_splits = 3, prune = FALSE, max_interaction = 1)
+  families <- fit$fit$get_model()
+  no_split <- vapply(families, function(f) identical(as.integer(unlist(f$variables)), 0L), logical(1))
+  expect_true(any(no_split) && !all(no_split))
+  # rule D at tau = 0 is the forest
+  forest <- prune_components(fit, d, rule = "D", tau = 0, max_panels = Inf)
+  expect_equal(predict(forest, d)$.pred, predict(fit, d)$.pred, tolerance = 1e-12)
+})
+
+test_that("the OOB path tracks a family without a split at its rows' fit, not at 0", {
+  # the path itself: the families' null leaves without the level (test-classif-level.R)
+  rlang::local_options(randomPlantedForest.level = FALSE)
+  dat <- no_split_data()
+  # 100 rows: no threshold keeps 51 rows on each side, so no family splits
+  path_fit <- function(y, n_families, ...) {
+    set.seed(2)
+    rpf(
+      stats::as.formula(paste(y, "~ x1 + x2")),
+      data = dat,
+      n_families = n_families,
+      n_splits = 5,
+      min_leaf_rows = 51,
+      prune = FALSE,
+      oob_path = c(1, 5),
+      ...
+    )
+  }
+  cases <- list(
+    list(y = "y"),
+    list(y = "cl", loss = "L2"),
+    list(y = "cl", loss = "logit"),
+    list(y = "cl3", loss = "logit"),
+    list(y = "cl", loss = "exponential"),
+    list(y = "cl3", loss = "exponential")
+  )
+  for (case in cases) {
+    label <- paste(case, collapse = " ")
+    # one family: every OOB row is at the family's fit, at every checkpoint
+    fit <- do.call(path_fit, c(list(case$y, 1), case[-1]))
+    families <- fit$fit$get_model()
+    expect_true(all(vapply(families, function(f) length(f$variables) == 1, TRUE)), label = label)
+    null_value <- families[[1]]$values[[1]][[1]]
+    V <- length(null_value)
+    path <- fit$fit$oob_path()
+    ok <- path$count[, 1] + path$count[, 2] > 0
+    expect_true(any(ok), label = label)
+    for (c in 1:2) {
+      cols <- (c - 1) * V + seq_len(V)
+      g <- path$sum0[ok, cols, drop = FALSE] + path$sum1[ok, cols, drop = FALSE]
+      expect_equal(g, matrix(null_value, sum(ok), V, byrow = TRUE), tolerance = 1e-12, label = label)
+    }
+  }
+
+  # regression: the OOB error is that of the families' bootstrap means (at 0 it was the mean of y, about 10.5)
+  fit <- path_fit("y", 10)
+  expect_lt(max(fit$oob_path$oob_rmse_raw), 2 * stats::sd(dat$y))
+  # the adaptive fit tracks them the same way: its first checkpoint's raw criterion is that of the families' fits
+  # (at the score 0 it was log(2))
+  set.seed(2)
+  ad <- rpf(cl ~ x1 + x2, data = dat, loss = "logit", n_families = 10, min_leaf_rows = 51, prune = FALSE)
+  p <- mean(dat$cl == "b")
+  expect_lt(ad$oob_path$oob_loss_raw[1], -(p * log(p) + (1 - p) * log(1 - p)) + 0.05)
+})
diff --git a/tests/testthat/test-oob-path.R b/tests/testthat/test-oob-path.R
index 84a5786..f0186a6 100644
--- a/tests/testthat/test-oob-path.R
+++ b/tests/testthat/test-oob-path.R
@@ -59,6 +59,8 @@ leaf_list <- function(fit) {
 }
 
 test_that("the prefix property: the k-split forest is the first k splits of a longer fit", {
+  # the growth, before the level a classification fit adds after its last split
+  rlang::local_options(randomPlantedForest.level = FALSE)
   # With delete_leaves = FALSE a split only appends leaves, so every tree of the k-split forest must be the first
   # leaves of that tree in the longer fit: every structure, regression and classification.
   d <- path_sim()
diff --git a/vignettes/articles/classification.Rmd b/vignettes/articles/classification.Rmd
index dfa611a..57ffab4 100644
--- a/vignettes/articles/classification.Rmd
+++ b/vignettes/articles/classification.Rmd
@@ -124,8 +124,10 @@ Two things follow:
   A round in which no candidate gains makes no split (it still counts towards `n_splits`).
 
 The table shows the score of a point in a pure region as one tree family grows: the two classes are separated by `x1 > 0.5`, and `x2` is noise.
+The fits leave out the forest's level (below), so that the table shows the steps alone.
 
 ```{r pure-region}
+old <- options(randomPlantedForest.level = FALSE)
 set.seed(1)
 n <- 200
 sep <- data.frame(x1 = runif(n), x2 = runif(n))
@@ -144,6 +146,7 @@ pure <- data.frame(
   t(sapply(rounds, function(k) score_after("logit", k)))
 )
 names(pure) <- c("rounds", "exponential score", "exponential leaves", "logit score", "logit leaves")
+options(old)
 pure
 ```
 
@@ -157,3 +160,10 @@ So a fixed, large `n_splits` does not make the model overconfident.
 With `exact_step = TRUE` (the default), the chosen split's two leaves then get the exact minimisers, by Newton's method on their rows; `exact_step = FALSE` keeps the one Newton step, which is a little faster.
 The search is the same either way, and `exact_step` does not change `"exponential"` fits.
 
+### The level
+
+The tree families are averaged on the score scale.
+Where they disagree, the probability of their mean score is not the mean of their probabilities: for a class rarer than 1/2 it is lower, and the more splits the families make, the more they disagree.
+So after the fit, a forest with `"logit"` or `"exponential"` gets one constant per score column, its level: the intercept that makes every class's mean out-of-bag probability equal to its frequency in the training rows (the maximum-likelihood intercept given the forest's out-of-bag scores).
+The level is part of the forest's intercept, so the components and the purified forest include it, and the pruning starts from the forest with it.
+`$level` holds it, `print()` shows it, and `options(randomPlantedForest.level = FALSE)` fits without it.