fix: rc7 review findings across the 4.0.0 RC series

DESCRIPTION

@@ -2,7 +2,7 @@ Package: ggRandomForests
 Type: Package
 Title: Visually Exploring Random Forests
 Version: 4.0.0
-Date: 2026-08-31
+Date: 2026-09-29
 Authors@R: person("John", "Ehrlinger",
   role = c("aut", "cre"),
   email = "john.ehrlinger@gmail.com",

NEWS.md

@@ -3,6 +3,18 @@ Version: 4.0.0
 
 ggRandomForests v4.0.0 (development)
 ====================================
+* `gg_vimp(nvar = )` on a multi-class `rfsrc` forest now keeps the `nvar`
+  most important variables. The importance matrix arrives in predictor order
+  and was trimmed before it was sorted, so `nvar = 2` on `iris` returned
+  `Sepal.Length` and `Sepal.Width`, the two least important. Variables are now
+  ranked by the overall (`all`) column first, as `which.outcome = 0` selects.
+  The `randomForest` method already did this.
+* `plot.gg_partial_varpro()` now draws only the `parametric` and
+  `nonparametric` curves by default. `causal` is a centered contrast that
+  starts at zero, so on a shared axis with the level curves it squeezed them
+  into a thin band; on a regression fit the levels sat near the response mean
+  and the contrast near 0. Ask for it with `type = "causal"` (or include it in
+  `type`), on an unbounded scale as before.
 * `gg_partial_varpro(scale = "prob")` now restores each subject's level
   before averaging, so the curve is the expected proportion it is documented
   as. `partialpro()` fits each subject's curve separately but returns every
@@ -64,7 +76,11 @@ ggRandomForests v4.0.0 (development)
   those regions. Supplying `xvar2.name` returns the second-order (interaction)
   surface for a pair of continuous predictors, which is zero everywhere when
   the two act additively. Survival forests are not supported, the same
-  limitation `gg_shap()` carries.
+  limitation `gg_shap()` carries. For a categorical predictor each step from
+  one level to the next averages the prediction change over the observations
+  at both levels, as Apley and Zhu define it. All the modified frames for a
+  variable, or for an interaction surface, go to the forest in one `predict()`
+  call.
 
   `plot()` and `autoplot()` methods draw the first-order curves as lines and
   bars faceted by variable, matching `plot.gg_partial_rfsrc()`'s layout so the
@@ -181,8 +197,12 @@ ggRandomForests v4.0.0 (development)
   frame `time / auc / se / lower / upper / marker` with an `iauc`
   attribute (Uno + standardized integrated AUC); `plot.gg_auct()` draws
   AUC(t) with a bootstrap CI ribbon when available and a 0.5 reference
-  line. `gg_auct.rhf(object, marker, auct_fit = NULL)` computes
-  `auct.rhf()` internally or reuses a cached fit.
+  line. `gg_auct.rhf(object, marker, auct_fit = NULL, method)` computes
+  `auct.rhf()` internally or reuses a cached fit. `method` chooses the
+  cumulative/dynamic (default) or incident/dynamic AUC; the result records it
+  as `attr(, "method")` and the plot names it on the y axis, since the two
+  curves otherwise look alike. A cached `auct_fit` decides the method, and a
+  `method` that contradicts it now warns rather than being dropped.
 * `gg_rhf_importance()` / `plot.gg_rhf_importance()`: tidy wrapper and point
   matrix for time-localized variable priority from
   `randomForestRHF::importance.rhf()` (RHF Phase 3). It returns

R/gg_ale_rfsrc.R

@@ -265,6 +265,27 @@ gg_ale_rfsrc <- function(rf_model,
   fJ - avg
 }
 
+## Predict several frames of the same columns with one predict() call and hand
+## the predictions back per frame. Each call to a forest's predict() pays a
+## fixed dispatch cost, which dominated ALE when every bin (and every corner of
+## every interaction cell) was predicted separately.
+.ale_predict_stacked <- function(pred_fun, frames) {
+  sizes <- vapply(frames, nrow, integer(1))
+  pred  <- pred_fun(do.call(rbind, frames))
+  split(pred, factor(rep(seq_along(frames), sizes), levels = seq_along(frames)))
+}
+
+## Mean of 'value' within each group 1..n_group; a group with no members is 0,
+## matching a bin that contributes no local effect. Only EMPTY groups are
+## zeroed: an occupied group whose predictions include NA stays NA, so a
+## prediction failure is visible rather than reported as a flat step.
+.ale_group_mean <- function(value, group, n_group) {
+  out <- tapply(value, factor(group, levels = seq_len(n_group)), mean)
+  out <- as.numeric(out)
+  out[tabulate(group, nbins = n_group) == 0L] <- 0
+  out
+}
+
 ## First-order ALE for one continuous predictor.
 .ale_continuous <- function(xname, newx, pred_fun, n_eval) {
   xval <- newx[[xname]]
@@ -276,18 +297,16 @@ gg_ale_rfsrc <- function(rf_model,
   n_bin     <- length(edges) - 1L
   bin   <- .ale_bin_index(xval, edges)
 
-  delta <- numeric(n_bin)
-  nk    <- numeric(n_bin)
-  for (k in seq_len(n_bin)) {
-    idx <- which(bin == k)
-    nk[k] <- length(idx)
-    if (nk[k] == 0L) next
-    dd_lo <- dd_all[idx, , drop = FALSE]
-    dd_hi <- dd_lo
-    dd_lo[[xname]] <- edges[k]
-    dd_hi[[xname]] <- edges[k + 1L]
-    delta[k] <- mean(pred_fun(dd_hi) - pred_fun(dd_lo))
-  }
+  ## Every observation sits in exactly one bin, so both edge frames can be
+  ## built for all rows at once and predicted in a single call, rather than two
+  ## predict() calls per bin.
+  dd_lo <- dd_all
+  dd_hi <- dd_all
+  dd_lo[[xname]] <- edges[bin]
+  dd_hi[[xname]] <- edges[bin + 1L]
+  pr    <- .ale_predict_stacked(pred_fun, list(dd_hi, dd_lo))
+  nk    <- tabulate(bin, nbins = n_bin)
+  delta <- .ale_group_mean(pr[[1]] - pr[[2]], bin, n_bin)
 
   ale <- .ale_accumulate(delta, nk)
   data.frame(x = edges, yhat = ale, name = xname, type = "continuous")
@@ -360,20 +379,21 @@ gg_ale_rfsrc <- function(rf_model,
   ## weighs on where the curve is centered.
   n_level <- tabulate(code, nbins = m)
 
-  n_bin <- m - 1L
-  delta <- numeric(n_bin)
-  for (k in seq_len(n_bin)) {
-    ## Bin k's members are observations at the UPPER level of the step, the
-    ## same convention as the continuous case (bin k = values up through
-    ## edge_k).
-    idx <- which(code == k + 1L)
-    if (length(idx) == 0L) next
-    dd_lo <- dd_all[idx, , drop = FALSE]
-    dd_hi <- dd_lo
-    dd_lo[[xname]] <- .ale_impose_level(dd_all[[xname]], fvalues[k])
-    dd_hi[[xname]] <- .ale_impose_level(dd_all[[xname]], fvalues[k + 1L])
-    delta[k] <- mean(pred_fun(dd_hi) - pred_fun(dd_lo))
-  }
+  ## Step k -> k+1 is averaged over the observations at BOTH levels, each moved
+  ## from level k to level k+1 (Apley and Zhu's categorical ALE, as ALEPlot
+  ## computes it). A continuous bin has observations on either side of its
+  ## edges; a level does not, so taking only the upper level's observations
+  ## would estimate a step from a handful of rows whenever that level is rare.
+  ## An observation at an interior level therefore serves two steps.
+  n_bin   <- m - 1L
+  members <- lapply(seq_len(n_bin), function(k) which(code %in% c(k, k + 1L)))
+  step    <- rep(seq_len(n_bin), lengths(members))
+  dd_lo   <- dd_all[unlist(members), , drop = FALSE]
+  dd_hi   <- dd_lo
+  dd_lo[[xname]] <- .ale_impose_level(dd_all[[xname]], fvalues[step])
+  dd_hi[[xname]] <- .ale_impose_level(dd_all[[xname]], fvalues[step + 1L])
+  pr    <- .ale_predict_stacked(pred_fun, list(dd_hi, dd_lo))
+  delta <- .ale_group_mean(pr[[1]] - pr[[2]], step, n_bin)
 
   ale <- .ale_accumulate_categorical(delta, n_level)
   data.frame(x = factor(flabels, levels = flabels), yhat = ale, name = xname,
@@ -439,30 +459,24 @@ gg_ale_rfsrc <- function(rf_model,
   b1 <- .ale_bin_index(x1, e1)
   b2 <- .ale_bin_index(x2, e2)
 
-  delta <- matrix(0, n_bin1, n_bin2)
-  cnt   <- matrix(0, n_bin1, n_bin2)
-  for (k in seq_len(n_bin1)) {
-    for (l in seq_len(n_bin2)) {
-      idx <- which(b1 == k & b2 == l)
-      cnt[k, l] <- length(idx)
-      if (cnt[k, l] == 0L) next
-      dd    <- dd_all[idx, , drop = FALSE]
-      dd_hh <- dd
-      dd_hh[[xname1]] <- e1[k + 1L]
-      dd_hh[[xname2]] <- e2[l + 1L]
-      dd_hl <- dd
-      dd_hl[[xname1]] <- e1[k + 1L]
-      dd_hl[[xname2]] <- e2[l]
-      dd_lh <- dd
-      dd_lh[[xname1]] <- e1[k]
-      dd_lh[[xname2]] <- e2[l + 1L]
-      dd_ll <- dd
-      dd_ll[[xname1]] <- e1[k]
-      dd_ll[[xname2]] <- e2[l]
-      delta[k, l] <- mean(pred_fun(dd_hh) - pred_fun(dd_hl) -
-                            pred_fun(dd_lh) + pred_fun(dd_ll))
-    }
+  ## Every observation sits in exactly one cell, so the four corner frames are
+  ## built for all rows at once and predicted in a single call; a default
+  ## 25 x 25 surface used to make up to 2,500 predict() calls. Cell (k, l) is
+  ## column-major index (l - 1) * n_bin1 + k, which is how delta is filled.
+  corner <- function(i1, i2) {
+    dd <- dd_all
+    dd[[xname1]] <- e1[i1]
+    dd[[xname2]] <- e2[i2]
+    dd
   }
+  pr <- .ale_predict_stacked(pred_fun, list(
+    corner(b1 + 1L, b2 + 1L), corner(b1 + 1L, b2),
+    corner(b1, b2 + 1L),      corner(b1, b2)))
+  cell  <- (b2 - 1L) * n_bin1 + b1
+  n_cel <- n_bin1 * n_bin2
+  delta <- matrix(.ale_group_mean(pr[[1]] - pr[[2]] - pr[[3]] + pr[[4]],
+                                  cell, n_cel), n_bin1, n_bin2)
+  cnt   <- matrix(tabulate(cell, nbins = n_cel), n_bin1, n_bin2)
 
   ## Double cumulative sum, down columns then across rows. Written as loops
   ## rather than nested apply(): apply() drops the dimension when an axis has a

R/gg_auct.R

@@ -25,8 +25,10 @@
 #' @return A `data.frame` of class `c("gg_auct", "data.frame")` with columns
 #'   `time`, `auc`, `se`, `lower`, `upper`, `marker` (CI columns `NA` when no
 #'   bootstrap), an `iauc` attribute (a list with `uno`, `std`, `uno.se`,
-#'   `std.se`, `conf.level`), and a `provenance` attribute derived from
-#'   `object` (source, family, ntree, n).
+#'   `std.se`, `conf.level`), a `method` attribute (`"cumulative"` or
+#'   `"incident"`, taken from the `auct.rhf` fit, which is what
+#'   [plot.gg_auct()] names on its y axis), and a `provenance` attribute
+#'   derived from `object` (source, family, ntree, n).
 #'
 #' @references
 #' Ishwaran H, Hsich EM, Kogalur UB, Lee DKK (2026). Random Hazard Forests.
@@ -67,6 +69,7 @@ gg_auct <- function(object, ...) {
 #' @export
 gg_auct.rhf <- function(object, marker = c("chf", "haz"), auct_fit = NULL,
                         method = c("cumulative", "incident"), ...) {
+  method_given <- !missing(method)
   marker <- match.arg(marker)
   method <- match.arg(method)
 
@@ -82,6 +85,15 @@ gg_auct.rhf <- function(object, marker = c("chf", "haz"), auct_fit = NULL,
     stop("auct_fit must be an 'auct.rhf' object from ",
          "randomForestRHF::auct.rhf().", call. = FALSE)
   }
+  ## A supplied fit decides the curve. Say so when the caller asked for the
+  ## other method, rather than returning a cumulative curve for an incident
+  ## request with nothing to show for it.
+  fit_method <- auct_fit$method %||% method
+  if (method_given && !identical(fit_method, method)) {
+    warning("gg_auct: auct_fit was computed with method = '", fit_method,
+            "', so method = '", method, "' is ignored. Recompute auct_fit, ",
+            "or leave auct_fit = NULL.", call. = FALSE)
+  }
 
   abt  <- auct_fit$AUC.by.time
   boot <- auct_fit$boot
@@ -103,6 +115,7 @@ gg_auct.rhf <- function(object, marker = c("chf", "haz"), auct_fit = NULL,
     std.se     = if (!is.null(boot)) boot$iAUC.std.se else NA_real_,
     conf.level = if (!is.null(boot)) boot$conf.level  else NA_real_
   )
+  attr(gg_dta, "method") <- fit_method
   class(gg_dta) <- c("gg_auct", class(gg_dta))
   gg_dta <- .set_provenance(gg_dta, object)
   invisible(gg_dta)

R/gg_vimp.R

@@ -298,8 +298,12 @@ gg_vimp.rfsrc <- function(object, nvar, ...) {
       colnames(gg_v)[1] <- which_col
       gg_dta <- gg_v
     } else {
-      # No specific class requested: attach variable names and pivot.
+      # No specific class requested: attach variable names and pivot. rfsrc's
+      # importance matrix is in predictor order, so rank by the overall
+      # measure (column 1, "all") before the nvar trim below keeps the top
+      # nvar rather than the first nvar.
       gg_dta$vars <- rownames(gg_dta)
+      gg_dta <- gg_dta[order(gg_dta[[1]], decreasing = TRUE), , drop = FALSE]
     }
 
     gg_dta <- gg_dta[seq_len(nvar), ]

R/plot.gg_auct.R

@@ -63,8 +63,19 @@ plot.gg_auct <- function(x, ...) {
     ggplot2::geom_line() +
     ggplot2::geom_hline(yintercept = 0.5, linetype = "dashed",
                         colour = "grey50") +
-    ggplot2::labs(x = "Time", y = "AUC(t)",
+    ggplot2::labs(x = "Time", y = .auct_ylabel(attr(x, "method")),
                   title = sprintf("Time-varying AUC (%s)", x$marker[1]),
                   caption = caption) +
     ggplot2::theme_bw()
 }
+
+## Name the AUC(t) estimand on the axis: a cumulative/dynamic and an
+## incident/dynamic curve otherwise render identically. A gg_auct object saved
+## before the method was recorded falls back to the bare label.
+#' @keywords internal
+.auct_ylabel <- function(method) {
+  switch(method %||% "",
+    cumulative = "Cumulative/dynamic AUC(t)",
+    incident   = "Incident/dynamic AUC(t)",
+    "AUC(t)")
+}

R/plot.gg_partial_varpro.R

@@ -82,8 +82,10 @@
 #'
 #' @param x A \code{\link{gg_partial_varpro}} object.
 #' @param type Character vector; one or more of \code{"parametric"},
-#'   \code{"nonparametric"}, \code{"causal"}.  Defaults to all three.
-#'   Ignored for path-C objects.
+#'   \code{"nonparametric"}, \code{"causal"}.  Defaults to the two level curves,
+#'   \code{"parametric"} and \code{"nonparametric"}.  \code{"causal"} is a
+#'   centered contrast rather than a level, so it has to be asked for, and it is
+#'   shown only on an unbounded scale.  Ignored for path-C objects.
 #' @param labels Optional variable labels for the facet strips.  One of: a named
 #'   character vector (\code{c(bpd_last = "BP Diastole")}); a labelled data frame,
 #'   whose \code{attr(col, "label")} values are read; or a two-column
@@ -512,6 +514,11 @@ plot.gg_partial_varpro <- function(x, # nolint: cyclocomp_linter
 ## back to the level curves if causal was the only requested type.
 #' @keywords internal
 .partial_varpro_plot_type <- function(type, type_user, prov) {
+  ## By default draw the level curves only. The causal contrast starts at zero
+  ## by construction, so on a shared axis with the levels it compresses them
+  ## into a thin band (a regression on sale price put the levels near 200,000
+  ## and the contrast near 0). Asking for it keeps it, where the scale allows.
+  if (!type_user) type <- setdiff(type, "causal")
   if (is.null(prov) || !.is_bounded_scale(prov$scale %||% "generic"))
     return(type)
   if (type_user && "causal" %in% type) {

R/utils.R

@@ -194,7 +194,7 @@ shift <- function(x, shift_by = 1) {
     return(nms)
   }
 
-  rank_key <- .varpro_rank_of(nms, ranked)
+  rank_key <- .varpro_rank_of(nms, ranked, xnames = colnames(object$x))
   ## seq_along() as the tiebreaker keeps incoming order stable among names that
   ## share a rank (in practice, all the unranked ones at Inf).
   nms[order(rank_key, seq_along(nms))]
@@ -203,12 +203,14 @@ shift <- function(x, shift_by = 1) {
 ## Position of each name in the ranked vector.  Exact match wins; failing that a
 ## one-hot level (name followed by digits, e.g. sex0/sex1) is accepted and the
 ## best -- lowest -- position across levels is taken.  Requiring digits keeps
-## 'age' from being captured by 'age_group'.
+## 'age' from being captured by 'age_group'.  The fallback is only for a name
+## that is not itself a column of the feature matrix ('xnames'): a real column
+## such as x1 must not borrow the rank of x10.
 #' @keywords internal
-.varpro_rank_of <- function(nms, ranked) {
+.varpro_rank_of <- function(nms, ranked, xnames = NULL) {
   vapply(nms, function(nm) {
     hit <- which(ranked == nm)
-    if (length(hit) == 0L) {
+    if (length(hit) == 0L && !(nm %in% xnames)) {
       pat <- paste0("^", .escape_regex(nm), "[0-9]+$")
       hit <- grep(pat, ranked)
     }