diff --git a/DESCRIPTION b/DESCRIPTION index e4242be0..7c0eebcb 100644 --- a/DESCRIPTION +++ b/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", diff --git a/NEWS.md b/NEWS.md index d1a61d5c..7894c2f4 100644 --- a/NEWS.md +++ b/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 diff --git a/R/gg_ale_rfsrc.R b/R/gg_ale_rfsrc.R index 57a1818c..73f977b3 100644 --- a/R/gg_ale_rfsrc.R +++ b/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 diff --git a/R/gg_auct.R b/R/gg_auct.R index d4c95664..fe027162 100644 --- a/R/gg_auct.R +++ b/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) diff --git a/R/gg_vimp.R b/R/gg_vimp.R index f1ebf07d..f9209fe9 100644 --- a/R/gg_vimp.R +++ b/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), ] diff --git a/R/plot.gg_auct.R b/R/plot.gg_auct.R index be4c9ddf..060a4698 100644 --- a/R/plot.gg_auct.R +++ b/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)") +} diff --git a/R/plot.gg_partial_varpro.R b/R/plot.gg_partial_varpro.R index 02bbd67b..18b9a83c 100644 --- a/R/plot.gg_partial_varpro.R +++ b/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) { diff --git a/R/utils.R b/R/utils.R index 39cabaef..f58017e8 100644 --- a/R/utils.R +++ b/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) } diff --git a/man/gg_auct.Rd b/man/gg_auct.Rd index d17e82e6..8cca634c 100644 --- a/man/gg_auct.Rd +++ b/man/gg_auct.Rd @@ -40,8 +40,10 @@ used when \code{auct_fit} is supplied, though the value is still validated.} A \code{data.frame} of class \code{c("gg_auct", "data.frame")} with columns \code{time}, \code{auc}, \code{se}, \code{lower}, \code{upper}, \code{marker} (CI columns \code{NA} when no bootstrap), an \code{iauc} attribute (a list with \code{uno}, \code{std}, \code{uno.se}, -\code{std.se}, \code{conf.level}), and a \code{provenance} attribute derived from -\code{object} (source, family, ntree, n). +\code{std.se}, \code{conf.level}), a \code{method} attribute (\code{"cumulative"} or +\code{"incident"}, taken from the \code{auct.rhf} fit, which is what +\code{\link[=plot.gg_auct]{plot.gg_auct()}} names on its y axis), and a \code{provenance} attribute +derived from \code{object} (source, family, ntree, n). } \description{ Extracts the time-dependent AUC curve from \code{\link[randomForestRHF:auct.rhf]{randomForestRHF::auct.rhf()}} diff --git a/man/plot.gg_partial_varpro.Rd b/man/plot.gg_partial_varpro.Rd index cc0a5d7d..bcae1c63 100644 --- a/man/plot.gg_partial_varpro.Rd +++ b/man/plot.gg_partial_varpro.Rd @@ -29,8 +29,10 @@ \item{x}{A \code{\link{gg_partial_varpro}} object.} \item{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.} \item{labels}{Optional variable labels for the facet strips. One of: a named character vector (\code{c(bpd_last = "BP Diastole")}); a labelled data frame, diff --git a/tests/testthat/_snaps/snapshots/gg-auct-chf.svg b/tests/testthat/_snaps/snapshots/gg-auct-chf.svg index d6fc7b02..dbd9e8d3 100644 --- a/tests/testthat/_snaps/snapshots/gg-auct-chf.svg +++ b/tests/testthat/_snaps/snapshots/gg-auct-chf.svg @@ -71,7 +71,7 @@ 2000 3000 Time -AUC(t) +Cumulative/dynamic AUC(t) Time-varying AUC (cumhaz) iAUC (Uno) = 0.818 | iAUC (standardized) = 0.661 diff --git a/tests/testthat/_snaps/snapshots/gg-partial-varpro-both.svg b/tests/testthat/_snaps/snapshots/gg-partial-varpro-both.svg index b0e68bbd..caa2f3d1 100644 --- a/tests/testthat/_snaps/snapshots/gg-partial-varpro-both.svg +++ b/tests/testthat/_snaps/snapshots/gg-partial-varpro-both.svg @@ -56,28 +56,26 @@ - - - - + + + - - - - + + + + - - - + + @@ -103,26 +101,23 @@ 60 70 80 --0.50 --0.25 -0.00 -0.25 - - - - +-0.50 +-0.25 +0.00 +0.25 + + + + Partial Effect - -Effect type - - - - - - -causal -nonparametric -parametric + +Effect type + + + + +nonparametric +parametric @@ -143,30 +138,22 @@ - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + @@ -195,32 +182,24 @@ Partial Effect - -Effect type - - - - - - - - - - - - - - - - - - - - - -causal -nonparametric -parametric + +Effect type + + + + + + + + + + + + + + +nonparametric +parametric gg-partial-varpro-both diff --git a/tests/testthat/_snaps/snapshots/gg-partial-varpro-categorical.svg b/tests/testthat/_snaps/snapshots/gg-partial-varpro-categorical.svg index aa48e362..24e202d3 100644 --- a/tests/testthat/_snaps/snapshots/gg-partial-varpro-categorical.svg +++ b/tests/testthat/_snaps/snapshots/gg-partial-varpro-categorical.svg @@ -27,30 +27,22 @@ - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + @@ -80,32 +72,24 @@ Partial Effect - -Effect type - - - - - - - - - - - - - - - - - - - - - -causal -nonparametric -parametric + +Effect type + + + + + + + + + + + + + + +nonparametric +parametric gg-partial-varpro-categorical diff --git a/tests/testthat/_snaps/snapshots/gg-partial-varpro-continuous.svg b/tests/testthat/_snaps/snapshots/gg-partial-varpro-continuous.svg index 750a252d..84ced623 100644 --- a/tests/testthat/_snaps/snapshots/gg-partial-varpro-continuous.svg +++ b/tests/testthat/_snaps/snapshots/gg-partial-varpro-continuous.svg @@ -27,9 +27,8 @@ - - - + + @@ -56,26 +55,23 @@ 60 70 80 --0.50 --0.25 -0.00 -0.25 - - - - +-0.50 +-0.25 +0.00 +0.25 + + + + Partial Effect - -Effect type - - - - - - -causal -nonparametric -parametric + +Effect type + + + + +nonparametric +parametric gg-partial-varpro-continuous diff --git a/tests/testthat/_snaps/snapshots/gg-partial-varpro-mortality.svg b/tests/testthat/_snaps/snapshots/gg-partial-varpro-mortality.svg index fb2bfaee..4353d75b 100644 --- a/tests/testthat/_snaps/snapshots/gg-partial-varpro-mortality.svg +++ b/tests/testthat/_snaps/snapshots/gg-partial-varpro-mortality.svg @@ -27,9 +27,8 @@ - - - + + @@ -56,26 +55,23 @@ 60 70 80 --0.50 --0.25 -0.00 -0.25 - - - - +-0.50 +-0.25 +0.00 +0.25 + + + + Ensemble mortality (expected events) - -Effect type - - - - - - -causal -nonparametric -parametric + +Effect type + + + + +nonparametric +parametric gg-partial-varpro-mortality diff --git a/tests/testthat/test_gg_ale_rfsrc.R b/tests/testthat/test_gg_ale_rfsrc.R index 1c5092ca..d82754f8 100644 --- a/tests/testthat/test_gg_ale_rfsrc.R +++ b/tests/testthat/test_gg_ale_rfsrc.R @@ -324,3 +324,51 @@ test_that("gg_ale_rfsrc does not count NA toward the cat_limit cardinality", { expect_gt(nrow(g$categorical), 0L) expect_equal(nrow(g$continuous), 0L) }) + +test_that("categorical ALE steps average over both adjacent levels", { + ## pred = z * level code, so the step k -> k+1 changes the prediction by z. + ## Averaging over cases at level k AND k+1 (Apley and Zhu) gives the mean z + ## of those two levels; the upper level alone gives only its own mean. + d <- data.frame(g = factor(c("a", "a", "a", "b", "c", "c")), + z = c(1, 1, 1, 10, 2, 2)) + pred_fun <- function(nd) nd$z * as.integer(nd$g) + out <- ggRandomForests:::.ale_categorical("g", d, pred_fun, + model_levels = c("a", "b", "c")) + step_ab <- mean(d$z[d$g %in% c("a", "b")]) # 13 / 4 + step_bc <- mean(d$z[d$g %in% c("b", "c")]) # 14 / 3 + expect_equal(diff(out$yhat), c(step_ab, step_bc)) +}) + +test_that("ALE predicts once per variable and once per interaction surface", { + ## Each bin (and each interaction cell corner) used to be its own predict() + ## call: 2,500 of them for a default interaction surface. Stack the frames and + ## predict once; the values must not change. + set.seed(9) + d <- data.frame(x1 = stats::runif(200), x2 = stats::runif(200)) + calls <- 0L + pred_fun <- function(nd) { + calls <<- calls + 1L + 2 * nd$x1 + nd$x2 + 3 * nd$x1 * nd$x2 + } + cont <- ggRandomForests:::.ale_continuous("x1", d, pred_fun, n_eval = 10) + expect_equal(calls, 1L) + ## Linear-plus-interaction: the first-order ALE slope for x1 is 2 + 3 E[x2|bin]. + expect_true(all(diff(cont$yhat) > 0)) + calls <- 0L + int <- ggRandomForests:::.ale_interaction(pred_fun, d, "x1", "x2", + cat_limit = 10, n_eval = 5) + expect_equal(calls, 1L) + expect_true(all(is.finite(int$ale))) +}) + +test_that("ALE keeps an NA prediction visible instead of zeroing its bin", { + ## Only an EMPTY bin contributes no local effect. A bin whose predictions + ## include NA must stay NA, as the per-bin loop left it, rather than be + ## reported as a flat step. + d <- data.frame(x1 = seq(0, 1, length.out = 40)) + pred_fun <- function(nd) ifelse(nd$x1 > 0.9, NA_real_, nd$x1) + cont <- ggRandomForests:::.ale_continuous("x1", d, pred_fun, n_eval = 4) + expect_true(anyNA(cont$yhat)) + expect_equal(ggRandomForests:::.ale_group_mean(c(1, NA, 3), c(1, 1, 3), 3), + c(NA, 0, 3)) +}) diff --git a/tests/testthat/test_gg_auct.R b/tests/testthat/test_gg_auct.R index bf8c644b..a216f16c 100644 --- a/tests/testthat/test_gg_auct.R +++ b/tests/testthat/test_gg_auct.R @@ -70,9 +70,26 @@ test_that("gg_auct passes ... through to auct.rhf", { expect_true(any(is.finite(gg$upper))) }) -test_that("gg_auct ignores method and ... when auct_fit is supplied", { +test_that("gg_auct uses auct_fit as given, and warns when method disagrees", { o <- .rhf_pbc() fit <- .auct_pbc_noboot() - gg <- gg_auct(o, marker = "chf", method = "incident", auct_fit = fit) + ## The fit decides the curve. A method that contradicts it used to be + ## dropped in silence, leaving a cumulative curve behind an incident request. + expect_warning( + gg <- gg_auct(o, marker = "chf", method = "incident", auct_fit = fit), + "auct_fit was computed with method = 'cumulative'") expect_equal(attr(gg, "iauc")$uno, fit$iAUC.uno) + expect_equal(attr(gg, "method"), "cumulative") + ## Leaving method at its default raises nothing. + expect_no_warning(gg_auct(o, marker = "chf", auct_fit = fit)) +}) + +test_that("gg_auct records the method and plot.gg_auct names it", { + o <- .rhf_pbc() + set.seed(20260828L) + inc <- gg_auct(o, marker = "haz", method = "incident") + expect_equal(attr(inc, "method"), "incident") + expect_match(plot(inc)$labels$y, "Incident/dynamic", fixed = TRUE) + cum <- gg_auct(o, marker = "chf", auct_fit = .auct_pbc_noboot()) + expect_match(plot(cum)$labels$y, "Cumulative/dynamic", fixed = TRUE) }) diff --git a/tests/testthat/test_gg_vimp.R b/tests/testthat/test_gg_vimp.R index ea43b77f..687beecb 100644 --- a/tests/testthat/test_gg_vimp.R +++ b/tests/testthat/test_gg_vimp.R @@ -536,6 +536,19 @@ test_that("gg_vimp: randomForest classification nvar counts variables, not rows" expect_setequal(unique(as.character(gg$vars)), top2) }) +test_that("gg_vimp: rfsrc classification nvar keeps the top nvar, not the first", { + set.seed(1) + rf <- randomForestSRC::rfsrc(Species ~ ., iris, ntree = 100, + importance = TRUE) + gg <- as.data.frame(gg_vimp(rf, nvar = 2)) + ## rfsrc's importance matrix is in predictor order, so trimming it unsorted + ## returned the two least important variables. Rank by the overall "all" + ## column, which is what which.outcome = 0 selects. + top2 <- names(sort(rf$importance[, "all"], decreasing = TRUE))[1:2] + expect_setequal(unique(as.character(gg$vars)), top2) + expect_equal(nrow(gg), 2L * ncol(rf$importance)) +}) + test_that("gg_vimp: randomForest classification importance=FALSE falls back to node purity", { skip_on_cran() diff --git a/tests/testthat/test_plot_gg_auct.R b/tests/testthat/test_plot_gg_auct.R index 45d61e16..caf6b15d 100644 --- a/tests/testthat/test_plot_gg_auct.R +++ b/tests/testthat/test_plot_gg_auct.R @@ -3,7 +3,7 @@ test_that("plot.gg_auct builds an AUC(t) ggplot with a 0.5 reference", { p <- plot(gg) expect_s3_class(p, "ggplot") expect_gt(nrow(ggplot2::layer_data(p)), 0) - expect_equal(p$labels$y, "AUC(t)") + expect_equal(p$labels$y, "Cumulative/dynamic AUC(t)") }) test_that("plot.gg_auct adds a ribbon when bootstrap CI is present (no warning)", { diff --git a/tests/testthat/test_varpro_importance_order.R b/tests/testthat/test_varpro_importance_order.R index 5dd6eac3..339a56af 100644 --- a/tests/testthat/test_varpro_importance_order.R +++ b/tests/testthat/test_varpro_importance_order.R @@ -59,3 +59,13 @@ test_that(".varpro_importance_order keeps list order among unranked names", { local_mocked_bindings(get.topvars = function(...) c("y"), .package = "varPro") expect_equal(.varpro_importance_order(pd, fake), c("y", "z", "x")) }) + +test_that(".varpro_rank_of does not let a real column borrow a longer name's rank", { + ## x1 is its own column of the feature matrix. When get.topvars() ranks x10 + ## but not x1, the digit-suffix fallback matched '^x1[0-9]+$' to x10 and put + ## x1 at x10's rank. The fallback is for expanded factors, whose own name is + ## absent from the feature matrix, so it applies only to those. + xn <- c("x1", "x10", "sex0", "sex1") + expect_equal(unname(.varpro_rank_of("x1", c("x10"), xnames = xn)), Inf) + expect_equal(unname(.varpro_rank_of("sex", c("age", "sex1"), xnames = xn)), 2) +}) diff --git a/vignettes/varpro.qmd b/vignettes/varpro.qmd index 6cb9cd6b..cb78e83a 100644 --- a/vignettes/varpro.qmd +++ b/vignettes/varpro.qmd @@ -211,12 +211,14 @@ gg_pd <- if (is.null(.vp$pd_boston)) { } else { .vp$pd_boston } -plot(gg_pd) +plot(gg_pd, type = c("parametric", "nonparametric", "causal")) ``` Each panel is a single predictor. The three curves correspond to the three estimators varPro carries (parametric, non-parametric, and causal); -read them as a sensitivity analysis. When all three agree, you have a +read them as a sensitivity analysis. `plot()` draws only the first two by +default, because `causal` is a contrast that starts at zero rather than a +level, so we ask for all three here. When all three agree, you have a stable signal; when the causal curve diverges from the others, that's a hint that the variable's observed relationship with the response may be partly driven by its correlation with other predictors in the rule