Package 'ggRandomForests'

Description

ggRandomForests provides ggplot2 (Wickham 2009) diagnostic and exploration figures for fits from randomForestSRC::rfsrc() (rfsrc class; randomForestSRC >= 3.4.0), randomForestRHF::rhf() (rhf class; randomForestRHF >= 2.0.0), varPro::varpro() (varpro class; varPro >= 3.1.0), or randomForest.

randomForestSRC gives a unified treatment of Breiman's (2001) random forests across data settings: regression and classification forests when the response is numeric or categorical, survival and competing-risk forests (Ishwaran et al. 2008) for right-censored data.

Details

The package is built on one decision: keep the data step and the figure step apart. A gg_* function pulls a tidy data object out of the forest; its plot() method turns that object into a figure. Two things follow.

The data object stands on its own. It carries everything its plot needs, so you can save it, inspect it, or come back to it later without keeping the original forest (which can be large) in memory.

You are never locked into the default figure. Each plot() method returns a single plottable object: a ggplot you extend with +, or a patchwork composite for the multi-panel methods. Add layers, swap scales, apply a theme, or ignore the default entirely and build the figure from the tidy data yourself. Every gg_* object also carries print() and summary() methods: print() shows a short header rather than dumping every row, and summary() returns a diagnostics object.

Forest diagnostics

• gg_rfsrc: predicted versus observed values.

• gg_error: OOB error against the number of trees.

• gg_vimp: variable importance ranking (Ishwaran et al. 2010).

• gg_variable: marginal variable dependence.

• gg_roc: ROC curves for classification forests (see also calc_roc and calc_auc).

• gg_survival: Kaplan-Meier / Nelson-Aalen estimates.

• gg_brier: time-resolved Brier score and CRPS for survival forests.

Partial dependence

• gg_partial: tidies the output of randomForestSRC::plot.variable(partial = TRUE).

• gg_partial_rfsrc: computes partial dependence from the fitted forest directly, via partial.rfsrc.

SHAP explanations

• gg_shap: SHAP values for regression and classification forests, wrapping kernelshap, with shap_importance, shap_beeswarm and shap_dependence figures.

varPro rule-based variable selection

• gg_varpro: variable importance from a varpro fit.

• gg_partial_varpro (alias gg_partialpro): partial effects on an interpretable scale.

• gg_beta_varpro: per-region lasso coefficients.

• gg_isopro: isolation-forest outlier scores.

• gg_ivarpro: individual (per-observation) importance.

Unsupervised varPro

• gg_udependent: variable dependency graph.

• gg_beta_uvarpro: entropy-based variable ranking.

• gg_sdependent: signal-variable detection.

Random Hazard Forests These functions are for time-to-event data whose predictors can change during follow-up. They show how event risk, time-varying discrimination, and variable priority change over time, and let you inspect a saved tree-size tuning search.

• gg_rhf: per-case event-risk curves over time.

• gg_auct: AUC curves over follow-up time.

• gg_rhf_importance: variable priority by time window.

• gg_tune_rhf: saved tree-size tuning paths by OOB risk or OOB iAUC.

See vignette("rhf", package = "ggRandomForests") for a longitudinal workflow that uses all four RHF families.

varPro is a required dependency (Imports), so the varPro families are always available. randomForestRHF and kernelshap are in Suggests: their functions check for the package and fail with a clear message when it is not installed.

References

Breiman, L. (2001). Random forests, Machine Learning, 45:5-32.

Ishwaran H, Kogalur U (2026). Fast Unified Random Forests for Survival, Regression, and Classification (RF-SRC). R package version 3.6.2. https://cran.r-project.org/package=randomForestSRC

Ishwaran H. and Kogalur U.B. (2007). Random survival forests for R. R News 7(2), 25–31.

Ishwaran H., Kogalur U.B., Blackstone E.H. and Lauer M.S. (2008). Random survival forests. Ann. Appl. Statist. 2(3), 841–860.

Ishwaran, H., U. B. Kogalur, E. Z. Gorodeski, A. J. Minn, and M. S. Lauer (2010). High-dimensional variable selection for survival data. J. Amer. Statist. Assoc. 105, 205-217.

Ishwaran, H. (2007). Variable importance in binary regression trees and forests. Electronic J. Statist., 1, 519-537.

Lu M, Ishwaran H (2024). Model-independent variable selection via the rule-based variable priority framework. arXiv:2409.09003. https://arxiv.org/abs/2409.09003

Ishwaran H, Kogalur U (2026). Model-Independent Variable Selection via the Rule-Based Variable Priority. R package version 3.2.0. https://CRAN.R-project.org/package=varPro

Ishwaran H, Hsich EM, Kogalur UB, Lee DKK (2026). Random Hazard Forests. arXiv:2608.21597. doi:10.48550/arXiv.2608.21597.

Ishwaran H, Kogalur UB (2026). randomForestRHF: Random Hazard Forests. R package version 2.0.3. https://CRAN.R-project.org/package=randomForestRHF

Wickham, H. ggplot2: elegant graphics for data analysis. Springer New York, 2009.

Description

These let you call ggplot2::autoplot() on any gg_* object ggRandomForests returns. Each is a thin wrapper around the matching plot.gg_*() S3 method, and ... passes straight through, so every argument those plot methods take is still available here.

Usage

## S3 method for class 'gg_error'
autoplot(object, ...)

## S3 method for class 'gg_vimp'
autoplot(object, ...)

## S3 method for class 'gg_rfsrc'
autoplot(object, ...)

## S3 method for class 'gg_variable'
autoplot(object, ...)

## S3 method for class 'gg_partial'
autoplot(object, ...)

## S3 method for class 'gg_partial_rfsrc'
autoplot(object, ...)

## S3 method for class 'gg_partialpro'
autoplot(object, ...)

## S3 method for class 'gg_partial_varpro'
autoplot(object, ...)

## S3 method for class 'gg_roc'
autoplot(object, ...)

## S3 method for class 'gg_survival'
autoplot(object, ...)

## S3 method for class 'gg_brier'
autoplot(object, ...)

## S3 method for class 'gg_rhf'
autoplot(object, ...)

## S3 method for class 'gg_tune_rhf'
autoplot(object, ...)

## S3 method for class 'gg_auct'
autoplot(object, ...)

## S3 method for class 'gg_rhf_importance'
autoplot(object, ...)

## S3 method for class 'gg_varpro'
autoplot(object, ...)

## S3 method for class 'gg_udependent'
autoplot(object, ...)

## S3 method for class 'gg_isopro'
autoplot(object, ...)

## S3 method for class 'gg_shap'
autoplot(object, ...)

## S3 method for class 'gg_ale_rfsrc'
autoplot(object, ...)

## S3 method for class 'gg_ale_interaction'
autoplot(object, ...)

Arguments

Details

The following gg_* classes are supported:

gg_error

OOB error vs. number of trees

gg_vimp

Variable importance ranking

gg_rfsrc

Predicted vs. observed values

gg_variable

Marginal dependence

gg_partial

Partial dependence (via plot.variable)

gg_partial_rfsrc

Partial dependence (via partial.rfsrc)

gg_partial_varpro

Partial dependence (via varPro)

gg_partialpro

Partial dependence via varPro (deprecated alias)

gg_varpro

Variable importance from varPro

gg_roc

ROC curve

gg_survival

Survival / cumulative hazard curves

gg_brier

Time-resolved Brier score and CRPS

gg_rhf

Random Hazard Forest hazard curves

gg_tune_rhf

Random Hazard Forest tuning path

gg_auct

Random Hazard Forest time-varying AUC

gg_rhf_importance

Random Hazard Forest variable priority

gg_shap

SHAP explanations

Value

A ggplot object.

Examples

library(ggplot2)
set.seed(42)
rf <- randomForestSRC::rfsrc(Ozone ~ ., data = na.omit(airquality),
                              ntree = 50, importance = TRUE,
                              tree.err = TRUE)
autoplot(gg_error(rf))
autoplot(gg_vimp(rf))

Description

Area Under the ROC Curve calculator

Usage

calc_auc(x)

Arguments

Details

calc_auc uses the trapezoidal rule to calculate the area under the ROC curve.

This is a helper function for the gg_roc functions.

Value

AUC. 50% is random guessing, higher is better.

See Also

calc_roc gg_roc

plot.gg_roc

Examples

##
## Taken from the gg_roc example
rfsrc_iris <- randomForestSRC::rfsrc(Species ~ ., data = iris, ntree = 100)

gg_dta <- gg_roc(rfsrc_iris, which_outcome = 1)

calc_auc(gg_dta)

gg_dta <- gg_roc(rfsrc_iris, which_outcome = 2)

calc_auc(gg_dta)

## randomForest tests
rf_iris <- randomForest::randomForest(Species ~ ., data = iris)
gg_dta <- gg_roc(rfsrc_iris, which_outcome = 2)

calc_auc(gg_dta)

Description

Receiver Operator Characteristic calculator

Usage

## S3 method for class 'rfsrc'
calc_roc(object, dta, which_outcome = "all", oob = TRUE, ...)

Arguments

Details

For a randomForestSRC prediction and the actual response value, calculate the specificity (1-False Positive Rate) and sensitivity (True Positive Rate) of a predictor.

This is a helper function for the gg_roc functions, and not intended for use by the end user.

Value

A gg_roc data.frame with columns sens (sensitivity), spec (specificity), and pct (the probability threshold), with one row per unique prediction value. Suitable for passing to calc_auc or plot.gg_roc.

See Also

calc_auc gg_roc

plot.gg_roc

Examples

## Taken from the gg_roc example
rfsrc_iris <- randomForestSRC::rfsrc(Species ~ ., data = iris, ntree = 100)

gg_dta <- calc_roc(rfsrc_iris, rfsrc_iris$yvar,
  which_outcome = 1, oob = TRUE
)
gg_dta <- calc_roc(rfsrc_iris, rfsrc_iris$yvar,
  which_outcome = 1, oob = FALSE
)

rf_iris <- randomForest::randomForest(Species ~ ., data = iris)
# randomForest stores the response in $y (rfsrc uses $yvar); pass the
# original training factor so calc_roc has the class labels.
gg_dta <- calc_roc(rf_iris, iris$Species,
  which_outcome = 1
)
gg_dta <- calc_roc(rf_iris, iris$Species,
  which_outcome = 2
)

Description

A partial dependence curve (gg_partial_rfsrc) marginalizes the forest's prediction by averaging over the joint distribution of the other predictors, a computation that is misleading when predictors are correlated, because it evaluates the forest at combinations of predictor values that never occur together in the data. Accumulated Local Effects (Apley and Zhu, 2020) avoid this by only ever perturbing a predictor within small local neighborhoods of its own observed values, then accumulating those local effects into a global curve.

Usage

gg_ale_rfsrc(
  rf_model,
  xvar.names,
  xvar2.name = NULL,
  newx = NULL,
  cat_limit = 10,
  n_eval = 25,
  which.class = 1
)

Arguments

Details

For a continuous predictor, gg_ale_rfsrc bins the observed values into n_eval quantile-based intervals. Within each bin, every observation's predictor value is replaced first with the bin's lower edge and then with its upper edge (all other predictors held at that observation's own values), and the average change in prediction is the bin's local effect. These local effects are accumulated (cumulatively summed) across bins and centered to have a weighted mean of zero, giving a curve that is directly comparable to gg_partial_rfsrc's output but immune to extrapolation into implausible predictor combinations.

For a categorical predictor, the model's own factor level order is used as the "grid": the local effect for the step from level $k$ to level $k+1$ is estimated only from observations actually at level $k+1$, comparing their prediction at that level against the counterfactual of level $k$. Relevel the predictor before fitting the forest to control this ordering.

Supplying xvar2.name switches to second-order (interaction) ALE between xvar.names and xvar2.name, isolating the part of their joint effect that is not explained by either variable's own main effect. It is the ALE analogue of an interaction term, computed on a 2-D grid of bins using the same local-perturbation idea, with the main effects removed via a row/column/grand weighted-mean decomposition (the same device used to isolate an interaction term in a two-way ANOVA). A purely additive forest returns an all-zero surface.

Value

For first-order ALE (xvar2.name = NULL), a named list with two elements, classed "gg_ale_rfsrc":

continuous

A data.frame with columns x (the bin edges, numeric), yhat (the centered ALE value at each edge), and name (variable name), for all continuous predictors.

categorical

The same columns but x kept as a factor in the model's level order, for low-cardinality predictors.

For second-order ALE (xvar2.name supplied), a data.frame classed "gg_ale_interaction" with columns x (grid values of xvar.names), y (grid values of xvar2.name), ale (the interaction surface value), name1, and name2.

References

Apley, D. W. and Zhu, J. (2020). Visualizing the effects of predictor variables in black box supervised learning models. Journal of the Royal Statistical Society Series B, 82(4), 1059-1086.

See Also

gg_partial_rfsrc, plot.gg_ale_rfsrc

Examples

## ------------------------------------------------------------
## regression, first-order ALE
## ------------------------------------------------------------
airq.obj <- randomForestSRC::rfsrc(Ozone ~ ., data = na.omit(airquality),
                                    ntree = 100)
ale_dta <- gg_ale_rfsrc(airq.obj, xvar.names = c("Wind", "Temp"))
plot(ale_dta)


## ------------------------------------------------------------
## second-order (interaction) ALE between two continuous predictors
## ------------------------------------------------------------
ale_int <- gg_ale_rfsrc(airq.obj, xvar.names = "Wind",
                         xvar2.name = "Temp", n_eval = 15)
plot(ale_int)

Description

Extracts the time-dependent AUC curve from randomForestRHF::auct.rhf() into a tidy long data frame, one row per time point, with bootstrap confidence bounds when available and the integrated AUC (iAUC) summary attached as an attribute.

Usage

gg_auct(object, ...)

## S3 method for class 'rhf'
gg_auct(
  object,
  marker = c("chf", "haz"),
  auct_fit = NULL,
  method = c("cumulative", "incident"),
  ...
)

Arguments

Value

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), 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).

Note

The two definitions answer different questions rather than better and worse versions of the same one. Cumulative/dynamic AUC ranks accumulated risk through a horizon, comparing subjects who have failed by that horizon against subjects still event-free at it. Incident/dynamic AUC ranks local failures within the risk set at each time. Pick the one that matches the question you are asking, and read the two curves as separate estimands rather than as a check on each other.

References

Ishwaran H, Hsich EM, Kogalur UB, Lee DKK (2026). Random Hazard Forests. arXiv:2608.21597. doi:10.48550/arXiv.2608.21597.

Ishwaran H, Kogalur UB (2026). randomForestRHF: Random Hazard Forests. R package version 2.0.3. https://CRAN.R-project.org/package=randomForestRHF.

See Also

plot.gg_auct(), randomForestRHF::auct.rhf()

Examples

if (requireNamespace("randomForestRHF", quietly = TRUE)) {
  data(pbc, package = "randomForestSRC")
  d <- randomForestRHF::convert.counting(
    survival::Surv(days, status) ~ ., na.omit(pbc))
  o <- randomForestRHF::rhf("Surv(id, start, stop, event) ~ .", d, ntree = 30)
  plot(gg_auct(o, marker = "chf"))
}

Description

Tidy wrapper around varPro::get.beta.entropy() for a uvarpro object. Where gg_beta_varpro() refines the supervised release-rule contrast, gg_beta_uvarpro() does the unsupervised analogue: uvarpro() builds entropy regions with no response, and get.beta.entropy() fits a cross-validated lasso within each region to ask how strongly every other variable explains the released variable. Averaging the absolute lasso coefficients per variable gives one number per variable: an unsupervised, lasso-flavoured importance.

Usage

gg_beta_uvarpro(object, ..., cutoff = NULL, beta_fit = NULL)

Arguments

Details

get.beta.entropy(o) returns a (released-variable x variable) numeric matrix of absolute lasso coefficients. The column mean (na.rm = TRUE) is the per-variable importance reported here, matching the canonical sort(colMeans(beta, na.rm = TRUE), decreasing = TRUE) idiom in the varPro::uvarpro() help ("iowa housing - illustrates lasso importance").

Because get.beta.entropy() is expensive (a cross-validated glmnet per region), the beta_fit argument accepts a pre-computed matrix so you can iterate on the cutoff without re-fitting. The pairing mirrors the beta_fit argument of gg_beta_varpro().

Value

A gg_beta_uvarpro object (a data.frame), one row per variable, most-important first, with columns:

variable

factor; levels reversed so the most-important variable lands at the top after coord_flip() (the gg_vimp convention).

beta_mean

⁠mean(|lasso beta|)⁠ over the released regions (colMeans(beta, na.rm = TRUE)).

n_released

number of regions contributing a non-NA coefficient for the variable.

selected

logical; beta_mean >= cutoff.

The provenance attribute records source, family ("unsupv"), cutoff, n_var, n_released_regions, and precomputed.

See Also

gg_beta_varpro() (supervised analogue), gg_udependent(), varPro::get.beta.entropy(), varPro::uvarpro().

Examples

if (requireNamespace("varPro", quietly = TRUE)) {
  set.seed(1)
  o <- varPro::uvarpro(mtcars, ntree = 50)
  gg <- gg_beta_uvarpro(o)
  plot(gg)
}

Description

Tidy wrapper around varPro::beta.varpro() for the regression or classification family. Aggregates the per-rule lasso coefficient $\hat{\beta}$ by variable into the mean absolute value $\mathrm{mean}(|\hat{\beta}|)$ and flags variables above a scalar cutoff. Optional beta_fit argument lets callers compute the expensive beta.varpro() step once and reuse the result.

Usage

gg_beta_varpro(object, ..., cutoff = NULL, beta_fit = NULL, which_class = NULL)

Arguments

Value

A data.frame of class c("gg_beta_varpro", "data.frame"). For a regression fit: one row per released variable, sorted by beta_mean descending. For a classification fit: long-format with an extra class column, one row per (variable, class) pair; variable is a factor whose levels run in ascending ⁠mean(|sum-of-class-beta|)⁠, most important last, so every facet / panel shares the same row order with the top variable at the top. which_class (or the binary default last-factor-level) collapses the output to a single class.

What this is doing

Think of the varPro release-rule mechanism as asking: "given a region of the feature space that the forest carved out, what changes when I remove the constraint on this one variable and let observations leave?" The standard importance answer (from gg_varpro()) measures that change as a z-scored contrast between local estimators: no synthetic data, no permutation. beta.varpro() asks the same question with a different ruler: for each rule (a tree-branch pair), it fits a one-predictor lasso regression of the response on the released variable's values, restricted to the OOB observations inside the rule's region. The wrapper aggregates those per-rule coefficients into one number per variable.

The key distinction from gg_vimp(), which measures Breiman-Cutler permutation importance by perturbing a variable's values and watching OOB error climb, is that neither gg_varpro() nor gg_beta_varpro() touches the data synthetically: all contrasts are between real subsets defined by the forest's rules.

What imp actually is (pedantic, because the column name is misleading)

The imp column on beta.varpro()'s ⁠$results⁠ is not a variable-importance score in the conventional sense. It is a regularised regression coefficient. Specifically:

• Per rule, glmnet fits a one-predictor lasso of the response on the released variable inside the rule's OOB region. use.cv = TRUE selects $\lambda$ by 10-fold CV (default nfolds = 10); use.1se = TRUE (default) picks lambda.1se. use.cv = FALSE uses the full $\lambda$ path.

• imp is the absolute fitted coefficient $|\hat{\beta}|$ at the chosen $\lambda$. There is no direction here: varPro::beta.varpro() wraps every coefficient it returns in abs(), so the sign is discarded upstream and never reaches us. Read imp as strength of local association, not its direction. (For a signed local estimator, see gg_ivarpro(), where the sign does survive.) Magnitude depends on the predictor's units (raw x, no standardization); a predictor in millimetres has a smaller $|\hat{\beta}|$ than the same predictor in metres.

• Lasso shrinkage can drive $\hat{\beta}$ to exactly zero. Those zeros are data, not missingness, and are kept in the aggregation. Convergence failures land as NA_real_ and are dropped.

• The per-variable aggregate is beta_mean, $\mathrm{mean}(|\hat{\beta}|)$, across the rules where this variable was released. It is not a permutation importance, not a split-strength importance, and not directly comparable on the same numeric axis to gg_varpro()'s z-scores. Disagreement with gg_varpro is often diagnostic, not a bug.

In code form: imp_r is the glmnet coefficient $\hat{\beta}$ fit on y ~ x_v restricted to rule r, with $\lambda$ chosen by use.cv / use.1se.

What's in the output

One row per released variable. Columns:

• variable: predictor name.

• beta_mean: mean of $|\hat{\beta}|$ across that variable's rules.

• n_rules: count of rules contributing (zero-beta rules included; only NA failures excluded).

• selected: beta_mean >= cutoff.

Provenance attribute carries source, family, ntree, cutoff, cutoff_default, use.cv, n_rules_total, n_rules_nonzero, precomputed, and xvar.names.

What you use this for

Picking variables when local effects matter more than aggregate split-strength contribution. Compare side-by-side with gg_varpro(): a variable that scores high here but low in gg_varpro is one whose local linear effect inside many rules is real even though its release-rule contrast is modest.

Caching

beta.varpro() is the expensive call (per-rule glmnet / cv.glmnet, often minutes on real data). Compute it once and reuse:

v <- varPro::varpro(mpg ~ ., data = mtcars, ntree = 200)
b <- varPro::beta.varpro(v, use.cv = TRUE)        # expensive, once
gg_a <- gg_beta_varpro(v, beta_fit = b)            # cheap
gg_b <- gg_beta_varpro(v, beta_fit = b, cutoff = 0.5)

Provenance carries precomputed = TRUE when beta_fit was supplied.

Classification

For a varpro classification fit (object$family == "class", binary or multi-class), the returned frame is long-format with an extra class column: one row per (variable, class) pair. The beta_mean column aggregates the per-class lasso coefficient $\hat{\beta}$ stored in beta.varpro()'s ⁠imp.<k>⁠ columns (one per class level). Same pedantic-beta semantics as regression, applied independently to each class.

Binary default: which_class = NULL resolves to the last factor level of the response, the positive-class convention used by glm and gg_roc. For a 30-day-mortality outcome with levels c("no", "yes"), that means the wrapper shows you "yes" (the event) by default.

Multi-class default: which_class = NULL returns all K classes; the plot method renders facet_wrap(~ class) with one cutoff line per facet.

which_class = "<name>" filters to a single class regardless of K. Errors if the name isn't in the response levels.

Per-class cutoffs: cutoff = NULL resolves to each class's mean(beta_mean). A scalar broadcasts. A named numeric vector overrides per class; missing names fall back to that class's mean.

Example (30-day mortality, binary):

fit <- varPro::varpro(event_30d ~ ., data = clinical, ntree = 200)
gg  <- gg_beta_varpro(fit)   # default: "yes" panel
plot(gg)

Reproducibility

Byte-for-byte agreement between cached (beta_fit = b) and uncached (beta_fit = NULL) outputs requires that b was computed by beta.varpro(object, ...) on the same object; set.seed() alone is not sufficient, because beta.varpro's internal cv.glmnet fits can pick slightly different folds across separate calls. Reuse beta_fit when reproducibility matters.

Note

Multivariate regression (⁠regr+⁠) and survival fits are not supported. The unsupported-family path reports the family it received.

See Also

gg_varpro(), gg_vimp(), plot.gg_beta_varpro(), varPro::beta.varpro().

Examples

if (requireNamespace("varPro", quietly = TRUE)) {
  set.seed(1)
  v <- varPro::varpro(mpg ~ ., data = mtcars, ntree = 50)
  b <- varPro::beta.varpro(v)
  gg <- gg_beta_varpro(v, beta_fit = b)
  plot(gg)
}

Description

The Brier score asks a familiar question of any probabilistic forecast: how far did the predicted probability sit from what actually happened? For a survival forest the forecast is the predicted survival probability at a given moment, and the "what happened" is whether the subject was still alive at that moment. The score is computed at every event time, so you get a curve rather than a single number, and lower is better everywhere. A perfectly calibrated forest that predicts 0 for every subject who died and 1 for every subject who survived would score 0; a forest that predicts 0.5 for everyone scores roughly 0.25 regardless of the true outcome. That is the "uninformative" reference, not a ceiling; a forest can score worse than it.

Usage

gg_brier(object, ...)

Arguments

Details

This function extracts the time-resolved Brier score for a survival forest grown with randomForestSRC, both overall and broken down by mortality-risk quartile (lowest-risk to highest-risk subjects). It also returns, in the crps column, the running continuous ranked probability score (CRPS): the Brier score integrated over time and divided by elapsed time, so each value summarizes calibration up to that point on the time axis.

Because subjects are right-censored, a plain Brier score is biased: censored subjects contribute no outcome information yet still inflate the denominator. The score here uses inverse-probability-of-censoring weighting (IPCW), which up-weights uncensored observations to compensate. The censoring distribution is estimated either by Kaplan-Meier (cens.model = "km", the default) or by a separate censoring forest (cens.model = "rfsrc") when the censoring mechanism is itself covariate-dependent.

Internally, this wraps get.brier.survival and rebuilds the quartile decomposition and running CRPS from the returned brier.matx and mort components, following the approach in the internal plot.survival of randomForestSRC.

Value

A gg_brier data.frame with columns

time

event time grid (object$time.interest).

brier

overall Brier score at each time.

bs.q25, bs.q50, bs.q75, bs.q100

Brier score within each mortality-risk quartile (lowest to highest risk).

bs.lower, bs.upper

15th and 85th percentile of per-subject Brier contributions at each time. Used by plot.gg_brier(by_quartile = TRUE) to draw an envelope around the overall curve.

crps

running CRPS (overall) at each time, normalized by elapsed time.

crps.q25, crps.q50, crps.q75, crps.q100

running CRPS within each mortality-risk quartile.

crps.lower, crps.upper

running CRPS of the 15th / 85th per-subject Brier percentile, normalized by elapsed time.

The integrated CRPS is attached as attr(., "crps_integrated"). It is get.brier.survival()$crps, the raw area under the Brier curve, not normalized by time, so it is in the units of the time axis and grows with follow-up. The time-normalized value, get.brier.survival()$crps.std (the raw integral divided by max(.$time)), is attached as attr(., "crps_std") and reads on the Brier scale. The last value of the crps column instead divides by the time elapsed since the first event time.

Note

Brier score / CRPS is randomForestSRC survival-only; there is no randomForest method.

References

Graf E., Schmoor C., Sauerbrei W., Schumacher M. (1999). Assessment and comparison of prognostic classification schemes for survival data. Statistics in Medicine, 18(17-18):2529-2545.

Gerds T.A., Schumacher M. (2006). Consistent estimation of the expected Brier score in general survival models with right-censored event times. Biometrical Journal, 48(6):1029-1040.

See Also

plot.gg_brier, get.brier.survival, gg_error

Examples

library(survival)   # Surv() must be on the search path for rfsrc()
data(pbc, package = "randomForestSRC")
rfsrc_pbc <- randomForestSRC::rfsrc(
  Surv(days, status) ~ ., data = pbc, nsplit = 10, ntree = 100
)
gg_dta <- gg_brier(rfsrc_pbc)
plot(gg_dta)
plot(gg_dta, type = "crps")
plot(gg_dta, envelope = TRUE)   # overall line + 15-85% envelope

# Multi-model comparison: stack gg_brier outputs and plot with ggplot2.
rf2 <- randomForestSRC::rfsrc(
  Surv(days, status) ~ ., data = pbc, nsplit = 10, mtry = 4, ntree = 100
)
compare_dta <- dplyr::bind_rows(
  dplyr::mutate(gg_brier(rfsrc_pbc), model = "default"),
  dplyr::mutate(gg_brier(rf2),       model = "mtry=4")
)
ggplot2::ggplot(compare_dta,
  ggplot2::aes(x = time, y = brier, color = model)) +
  ggplot2::geom_line()

Description

Extract the cumulative out-of-bag (OOB) or in-bag training error rate from randomForestSRC and randomForest fits as a function of the number of grown trees.

Usage

gg_error(object, ...)

Arguments

Details

You have to ask rfsrc for the trajectory. How many trees you get a curve over is decided by randomForestSRC::rfsrc(), not here. Its block.size argument controls how often the error is recorded, and it defaults to NULL unless you request importance; NULL records the error at the final tree only. So a default fit gives gg_error() a single point, not a curve, and tree.err = TRUE on its own does not change that. Grow the forest with block.size = 1 for an error recorded at every tree, or a larger block.size for every nth. If a plot.gg_error() figure comes back as one dot, this is why.

For randomForestSRC objects the function reshapes the rfsrc$err.rate matrix, drops the trees where rfsrc recorded nothing, and annotates the rest with the tree index required by plot.gg_error. When supplied a randomForest object, the method inspects either the $mse or $err.rate component and, when training = TRUE is requested, reconstructs the original training set via the model call to compute an in-bag error curve using per-tree predictions. Training curves are only available when the forest was stored (keep.forest = TRUE) and the original data can be recovered.

Value

A gg_error data.frame containing at least the cumulative OOB error columns and an ntree counter. When training = TRUE is honored an additional train column is included.

References

Breiman L. (2001). Random forests, Machine Learning, 45:5-32.

Ishwaran H. and Kogalur U.B. (2007). Random survival forests for R, Rnews, 7(2):25-31.

Ishwaran H, Kogalur U (2026). Fast Unified Random Forests for Survival, Regression, and Classification (RF-SRC). R package version 3.6.2. https://cran.r-project.org/package=randomForestSRC

See Also

plot.gg_error, gg_vimp, gg_variable, rfsrc, randomForest, plot.rfsrc

Examples

## Examples from RFSRC package...
## ------------------------------------------------------------
## classification example
## ------------------------------------------------------------
## ------------- iris data
## You can build a randomForest
rfsrc_iris <- randomForestSRC::rfsrc(Species ~ ., data = iris,
  ntree = 100, tree.err = TRUE, block.size = 1)

# Get a data.frame containing error rates
gg_dta <- gg_error(rfsrc_iris)

# Plot the gg_error object
plot(gg_dta)

## RandomForest example
rf_iris <- randomForest::randomForest(Species ~ .,
  data = iris,
  tree.err = TRUE
)
gg_dta <- gg_error(rf_iris)
plot(gg_dta)

gg_dta <- gg_error(rf_iris, training = TRUE)
plot(gg_dta)
## ------------------------------------------------------------
## Regression example
## ------------------------------------------------------------

## ------------- airq data
rfsrc_airq <- randomForestSRC::rfsrc(Ozone ~ .,
  data = airquality, ntree = 100,
  na.action = "na.impute", tree.err = TRUE, block.size = 1
)

# Get a data.frame containing error rates
gg_dta <- gg_error(rfsrc_airq)

# Plot the gg_error object
plot(gg_dta)


## ------------- Boston data
if (requireNamespace("MASS", quietly = TRUE)) {
  data(Boston, package = "MASS")
  Boston$chas <- as.logical(Boston$chas)
  rfsrc_boston <- randomForestSRC::rfsrc(medv ~ .,
    data = Boston,
    ntree = 100,
    forest = TRUE,
    importance = TRUE,
    tree.err = TRUE,
    save.memory = TRUE
  )

  # Get a data.frame containing error rates
  gg_dta <- gg_error(rfsrc_boston)

  # Plot the gg_error object
  plot(gg_dta)
}


## ------------- mtcars data
rfsrc_mtcars <- randomForestSRC::rfsrc(mpg ~ ., data = mtcars,
  ntree = 100, tree.err = TRUE, block.size = 1)

# Get a data.frame containing error rates
gg_dta<- gg_error(rfsrc_mtcars)

# Plot the gg_error object
plot(gg_dta)


## ------------------------------------------------------------
## Survival example
## ------------------------------------------------------------
## ------------- veteran data
## randomized trial of two treatment regimens for lung cancer
data(veteran, package = "randomForestSRC")
rfsrc_veteran <- randomForestSRC::rfsrc(Surv(time, status) ~ ., data = veteran,
                       ntree = 100, tree.err = TRUE, block.size = 1)

gg_dta <- gg_error(rfsrc_veteran)
plot(gg_dta)

## ------------- pbc data
data(pbc, package = "randomForestSRC", envir = environment())
# Recode: binary columns to logical, low-cardinality columns to factors.
for (ind in seq_len(dim(pbc)[2])) {
  if (!is.factor(pbc[, ind])) {
    if (length(unique(pbc[which(!is.na(pbc[, ind])), ind])) <= 2) {
      if (sum(range(pbc[, ind], na.rm = TRUE) == c(0, 1)) == 2) {
        pbc[, ind] <- as.logical(pbc[, ind])
      }
    }
  } else {
    if (length(unique(pbc[which(!is.na(pbc[, ind])), ind])) <= 2) {
      if (sum(sort(unique(pbc[, ind])) == c(0, 1)) == 2) {
        pbc[, ind] <- as.logical(pbc[, ind])
      }
      if (sum(sort(unique(pbc[, ind])) == c(FALSE, TRUE)) == 2) {
        pbc[, ind] <- as.logical(pbc[, ind])
      }
    }
  }
  if (!is.logical(pbc[, ind]) &&
    length(unique(pbc[which(!is.na(pbc[, ind])), ind])) <= 5) {
    pbc[, ind] <- factor(pbc[, ind])
  }
}
# Convert age from days to years
pbc$age <- pbc$age / 364.24
pbc$years <- pbc$days / 364.24
pbc <- pbc[, -which(colnames(pbc) == "days")]
pbc$treatment <- as.numeric(pbc$treatment)
pbc$treatment[which(pbc$treatment == 1)] <- "DPCA"
pbc$treatment[which(pbc$treatment == 2)] <- "placebo"
pbc$treatment <- factor(pbc$treatment)
# Split into a training set (assigned treatment) and a test set (not assigned)
dta_train <- pbc[-which(is.na(pbc$treatment)), ]
pbc_test <- pbc[which(is.na(pbc$treatment)), ]

#========
# build the forest:
rfsrc_pbc <- randomForestSRC::rfsrc(
  Surv(years, status) ~ .,
 dta_train,
 ntree = 100,
 nsplit = 10,
 na.action = "na.impute",
 tree.err = TRUE,
 forest = TRUE,
 importance = TRUE,
 save.memory = TRUE
)


gg_dta <- gg_error(rfsrc_pbc)
plot(gg_dta)

Description

Pulls per-observation anomaly scores out of a isopro fit so you can plot them, sort them, or write them to disk without having to know the internal shape of the fit.

Usage

gg_isopro(object, ..., newdata = NULL)

Arguments

Value

A data.frame of class c("gg_isopro", "data.frame"), one row per observation. Columns:

obs

Integer; observation index 1..n, in the same order as the rows of the data passed to isopro.

case.depth

Numeric; mean isolation depth across the forest. Lower means the observation was isolated quickly, so more anomalous.

howbad

Numeric in [0, 1]; the case.depth values pushed through their own empirical CDF and flipped so higher means more anomalous. This is the score the plot method draws by default.

A provenance attribute records source = "varPro::isopro", the observation count n, and the number of trees ntree.

What isopro is doing

An isolation forest (Liu, Ting and Zhou 2008) is a random forest grown on very small subsamples of the data and asked to split until each observation lands in its own terminal node. The intuition is geometric: a typical observation sits in the dense middle of the feature cloud and takes many splits to isolate, while an unusual observation sits out near an edge and gets cut off after only a few. So the depth at which an observation is isolated is a proxy for how typical it is: shallow depth means anomalous, deep depth means ordinary. Average a single observation's depth across many trees and the noise washes out, leaving a stable per-observation rank.

isopro supports three flavours of isolation forest, which differ in how the splits are chosen:

"rnd"

The original Liu/Ting/Zhou method: each tree node picks a variable at random and a split point uniformly at random in the variable's range. Fast, no model, surprisingly effective.

"unsupv"

Unsupervised splitting from randomForestSRC: at each split a random subset of the features serves as pseudo-responses (ytry of them, about the square root of the number of features by default) and the split is chosen to separate those pseudo-responses. More structured than "rnd"; sometimes more accurate, especially when the anomalies follow a coherent direction.

"auto"

An auto-encoder formulation that grows a multivariate forest predicting each feature from the others. Most expressive, slowest, best suited to low-dimensional data.

No method is universally best. The varPro authors recommend trying at least two and comparing the score distributions; the plot method here colors per-method curves automatically when you stack the results.

What's in the output

The fit gives back two parallel per-observation vectors: case.depth is the raw mean isolation depth (units of "splits", lower = more anomalous) and howbad is the same information on a [0, 1] scale, the empirical CDF of case.depth. On the fit those two run the same way: varPro::isopro()'s howbad is lower = more anomalous, like the depth it came from.

We flip it. gg_isopro() returns 1 - object$howbad, so the column you get is higher = more anomalous and reads the way a score should. That is our transformation, not the fit's, and it means gg$howbad is 1 - fit$howbad rather than a copy of it, worth knowing if you compare the two side by side. Both columns are kept so you can plot in either space and have the raw depth on hand for diagnostics; howbad is the canonical score and is what the plot method uses by default.

What you use this for

This is screening, not inference. Reach for it when you want to:

• flag observations that may be data-entry errors, out-of-range measurements, or distinct subpopulations before fitting a primary model;

• check whether a held-out cohort sits inside the training distribution before scoring with a model trained elsewhere;

• give the analyst a ranked list of "look at these first" cases for a manual review;

• score a held-out cohort or a fresh batch of incoming data against a fitted model and compare the test scores to the training distribution.

The score is a rank, not a probability of being an outlier: two observations with howbad = 0.92 are both unusual, not "92% likely to be anomalous". Pick a cutoff by looking at where the elbow rises; plot.gg_isopro can annotate either a score (threshold) or a top-percent (top_n_pct) for you.

Scoring new data

Pass a data.frame as newdata and the extractor calls predict.isopro twice: once with quantiles = FALSE to get the raw mean case depth per row, and once with quantiles = TRUE to get the per-row quantile of that depth against the training-data depth distribution.

varPro's predict.isopro returns quantiles where smaller is more anomalous, which is the opposite polarity of the wrapper's howbad (where higher is more anomalous). The wrapper exposes both conventions so nothing is hidden:

• case.depth carries varPro's native polarity, lower = more anomalous. This is the unmodified output of predict(object, newdata, quantiles = FALSE). Use it to cross-reference against raw varPro output.

• howbad is the flipped, wrapper-convention version. The relationship is howbad = 1 - predict(object, newdata, quantiles = TRUE).

To overlay training and test scores in one plot, bind the two extractor calls with a method label column (the same column plot.gg_isopro uses to color rnd / unsupv / auto comparisons):

gg_train <- gg_isopro(fit)
gg_test  <- gg_isopro(fit, newdata = test_df)
gg_both  <- rbind(cbind(gg_train, method = "train"),
                  cbind(gg_test,  method = "test"))
class(gg_both) <- c("gg_isopro", "data.frame")
plot(gg_both)

Comparing methods

To compare methods ("rnd", "unsupv", "auto"), call gg_isopro on each fit and dplyr::bind_rows() the results with a method label column. The plot method auto-detects method and colors the curves.

References

Liu, F. T., Ting, K. M., and Zhou, Z. H. (2008). Isolation Forest. Eighth IEEE International Conference on Data Mining, 413-422.

Ishwaran H, Kogalur U (2026). varPro: Model-Independent Variable Selection via the Rule-Based Variable Priority. R package version 3.2.0. https://CRAN.R-project.org/package=varPro

See Also

plot.gg_isopro, isopro

Examples

if (requireNamespace("varPro", quietly = TRUE)) {
  set.seed(1)
  fit <- varPro::isopro(data = iris[, 1:4], method = "rnd",
                        sampsize = 32, ntree = 50)
  gg <- gg_isopro(fit)
  plot(gg)
}

Description

Tidy wrapper around varPro::ivarpro() for the regression or classification family. Returns one row per (observation, variable) pair where the local-importance cell is non-NA; classification adds a class column. which_obs collapses to a per-observation profile; which_class collapses to a single class. Optional ivarpro_fit argument lets callers cache the expensive ivarpro() call.

Usage

gg_ivarpro(
  object,
  ...,
  which_obs = NULL,
  which_class = NULL,
  cutoff = NULL,
  ivarpro_fit = NULL
)

Arguments

Value

A data.frame of class c("gg_ivarpro", "data.frame"). Regression: columns obs / variable / local_imp / selected. Classification: long-format with an extra class column. variable is a factor whose levels run in ascending ⁠mean(|local_imp|)⁠ across all rows, most important last (the unified ranking axis shared across facets / panels).

What this is doing

The varPro framework builds importance from release rules: for a given rule region, it compares a local estimator inside that region to what the estimator becomes after the constraint on the tested variable is removed ("released"). That contrast is summed over many rules and trees to get a global z-score: the quantity gg_varpro() shows. What ivarpro() adds is a per-observation view of the same mechanism.

Concretely: ivarpro() walks the forest's rules and, for each (observation, variable) pair, computes a scaled per-rule contribution to predicting that observation. Per-rule LOO removes the observation from its own rule before scoring, so the contribution is not inflated by the observation having helped define the region. Per-region scaling (scale = "local", default) standardizes the contribution by the rule's local response standard deviation so values are comparable across rules of different size. Aggregating those per-rule scores into one number per (obs, variable) pair gives the local_imp cell.

No permutation, no synthetic data: the contrast is always between real subsets of the observed data, defined by the forest's own rules. This is the same no-synthetic-features property that distinguishes gg_varpro() from gg_vimp()'s Breiman-Cutler permutation importance.

What local_imp actually is (pedantic)

local_imp[i, v] is the scaled aggregated rule contribution of variable v to predicting observation i, NOT a permutation importance and NOT a SHAP value. Sign carries direction of the local response shift inside the rule's region. Magnitude is on the response scale when scale = "global", or unit-free when scale = "local" (the default). The matrix is heavily sparse: an observation contributes only to rules that retain it as OOB; on real data, per-variable NA fractions of 50-95% are common. Comparison with gg_varpro() (aggregate split-strength) and gg_beta_varpro() (per-rule lasso beta) is diagnostic: a variable that's important globally but has low per-observation contribution for a specific case is interesting; the inverse (high local but low global) flags a regime-specific signal.

What's in the output

Long-format tidy frame. Regression has columns obs, variable, local_imp, selected. Classification adds a class column (factor in response-level order). variable is a factor whose levels run in ascending ⁠mean(|local_imp|)⁠ across all rows, so the most important variable is the last level and plots at the top; for classification that aggregate is across all (obs, class) so every facet / panel shows variables in the same row order. NA cells are filtered out; the source matrix is sparse, and the tidy frame only carries the cells where local importance is defined.

Provenance attribute carries source, family, ntree, cutoff (named numeric vector: length 1 named "regr" for regression, length K named with class levels for classification), cutoff_default, use.loo, scale, n_train, n_obs, n_var, precomputed, xvar.names, class_levels (classification only), which_obs, which_class.

What you use this for

Per-observation interpretation ("which variables drive this prediction?"), variable-selection diagnostics via the aggregate distribution view, and side-by-side comparison against gg_varpro() / gg_beta_varpro() to spot variables that matter locally but not globally (or vice versa).

Caching

ivarpro() is the most expensive call in varPro (per-rule leave-one-out + per-region scaling, often minutes on real data). Compute it once and reuse:

v   <- varPro::varpro(medv ~ ., data = Boston, ntree = 200)
iv  <- varPro::ivarpro(v, scale = "local")              # expensive, once
gg_aggregate <- gg_ivarpro(v, ivarpro_fit = iv)          # cheap
gg_case1     <- gg_ivarpro(v, ivarpro_fit = iv, which_obs = 1L)

Provenance carries precomputed = TRUE when ivarpro_fit was supplied.

Classification

For a classification fit, ivarpro() returns a list of K matrices (one per class) for multi-class, or a flat data.frame for binary. The binary frame explains one class, named by its target attribute and by default the first factor level. Because the two class probabilities sum to one, the wrapper derives the other class exactly, as the negated local slopes (or the same values under use.abs = TRUE), so both levels are available. The wrapper stacks per-class frames into a long-format frame with a class column. which_class = NULL returns all classes (binary defaults to the last factor level, the positive-class convention used by glm and gg_roc); which_class = "<name>" filters to a single class. cutoff polymorphism mirrors gg_beta_varpro(): NULL is per-class mean(|local_imp|), a scalar broadcasts, a named numeric vector overrides per class with fallback to that class's mean.

Reproducibility

Byte-for-byte agreement between cached (ivarpro_fit = iv) and uncached (ivarpro_fit = NULL) outputs requires reusing the same ivarpro() result. set.seed() alone is not sufficient because per-rule LOO subsampling can drift across separate calls. Reuse ivarpro_fit when reproducibility matters.

Note

Multivariate regression (⁠regr+⁠) and survival fits are not supported. The unsupported-family path reports the family it received.

See Also

gg_varpro(), gg_vimp(), gg_beta_varpro(), varPro::ivarpro().

Examples

if (requireNamespace("varPro", quietly = TRUE) &&
    requireNamespace("MASS", quietly = TRUE)) {
  set.seed(1)
  v <- varPro::varpro(medv ~ ., data = MASS::Boston, ntree = 50)
  iv <- varPro::ivarpro(v)
  gg <- gg_ivarpro(v, ivarpro_fit = iv)
  plot(gg)
}

Description

A partial dependence curve answers a what-if question about a forest: hold every other predictor at its observed value, sweep one of them across its range, and watch how the ensemble prediction moves. Marginalized over the joint distribution of the other variables, the resulting curve isolates the average effect of the swept predictor alone.

Usage

gg_partial(part_dta, nvars = NULL, cat_limit = 10, model = NULL)

Arguments

Details

gg_partial handles the bookkeeping step after you've already called randomForestSRC::plot.variable(partial = TRUE): it takes the list that function returns and separates the variables into two tidy data frames, one for continuous predictors (plotted as lines) and one for categorical predictors (plotted as box plots). The split is controlled by cat_limit: variables with more unique x-values than this threshold are treated as continuous; all others are categorical.

If you'd rather skip the plot.variable step and pass the fitted forest directly, see gg_partial_rfsrc, which calls partial.rfsrc for you.

Value

A named list with two elements:

continuous

data.frame with columns x, yhat, name (and optionally model) for continuous variables

categorical

data.frame with the same columns but with x as a factor, for low-cardinality / categorical variables. Unlike continuous, it is not averaged: plot.variable() returns one prediction per training observation per level, and plot.gg_partial draws their spread as box plots. Average yhat within x for the partial dependence value.

Note

Partial-dependence extraction is randomForestSRC-only; there is no randomForest method (the randomForest package provides no comparable partial-dependence interface).

For survival forests, randomForestSRC::plot.variable defaults to surv.type = "mort", so yhat is mortality (the expected number of events) and not a survival probability. It is therefore not on $[0, 1]$ and is not directly comparable with the survival probabilities returned by gg_variable. For a comparable quantity, ask for it explicitly: randomForestSRC::plot.variable(rf, partial = TRUE, surv.type = "surv"). The label describing the plotted quantity is recorded on the returned object as attr(x, "ylabel") and is used as the y-axis title by plot.gg_partial.

gg_partial_rfsrc defaults to partial.type = "surv" and so reports survival probabilities. The two entry points therefore report different quantities by default.

See Also

gg_partial_rfsrc gg_partialpro

Examples

## Build a small regression forest on the airquality dataset
set.seed(42)
airq <- na.omit(airquality)
rf <- randomForestSRC::rfsrc(Ozone ~ ., data = airq, ntree = 50)

## Compute partial dependence via plot.variable (show.plots = FALSE to
## suppress the base-graphics output, we only want the data)
pv <- randomForestSRC::plot.variable(rf, partial = TRUE,
                                      show.plots = FALSE)

## Split into continuous and categorical data frames
result <- gg_partial(pv)
head(result$continuous)

## Label this model for later comparison with a second forest
result_labeled <- gg_partial(pv, model = "airq_model")
unique(result_labeled$continuous$model)

Description

A partial dependence curve marginalizes the forest's prediction over all other predictors: for each evaluation point of the target variable, the forest scores every training observation with that value substituted in, then averages the result. What you get is the average effect of the target variable after "integrating out" the rest, a curve that would be flat if the variable carried no signal.

Usage

gg_partial_rfsrc(
  rf_model,
  xvar.names = NULL,
  xvar2.name = NULL,
  newx = NULL,
  partial.time = NULL,
  partial.type = c("surv", "chf", "mort"),
  cat_limit = 10,
  n_eval = 25
)

Arguments

Details

This function builds those curves for one or more predictors by calling partial.rfsrc and then tidy-stacking the results into separate data frames for continuous and categorical variables. Unlike gg_partial (which wraps plot.variable), you pass the fitted rfsrc object directly, with no intermediate plot.variable step.

For survival forests, the marginalized quantity depends on partial.type: survival probability ("surv"), cumulative hazard function ("chf"), or expected mortality ("mort"). You can request the curve at one or more time horizons via partial.time; the resulting data have a time column so the plot layers them as separate colored lines. Mortality is the exception: it is summed over every event time, so partial.time has no effect on it and the data carry no time column.

Value

A named list with two elements:

continuous

A data.frame with columns x (numeric), yhat, name (variable name), and optionally grp (the level of xvar2.name) and time (survival forests with partial.type "surv" or "chf") for all continuous predictors.

categorical

A data.frame with the same columns but x kept as a factor (levels in the model's level order, not alphabetical), for low-cardinality predictors. Unlike continuous, it is not averaged: there is one row per training observation per level (and per time point), so plot.gg_partial_rfsrc draws it as box plots. Average yhat within x for the partial dependence value.

Survival forests and partial.time

partial.rfsrc expects every value in partial.time to be an exact member of the model's time.interest vector, the unique observed event times stored in the fitted object. Pass an arbitrary time, even a plausible one such as c(1, 3) for a study measured in years, and you get a C-level prediction error from inside partial.rfsrc.

gg_partial_rfsrc takes care of this: every element of partial.time is silently snapped to its nearest time.interest value before the call. To target a specific follow-up horizon, find the closest grid point yourself and pass it explicitly:

ti  <- rf_model$time.interest
t1  <- ti[which.min(abs(ti - 1))]   # nearest to 1 year
pd  <- gg_partial_rfsrc(rf_model, xvar.names = "x", partial.time = t1)

The snap is silent, and it matters most before the first event: every request earlier than the first event time returns that first event's curve, so partials at, say, 10, 30 and 60 days come back identical in a sample whose first death is later. Check the time column of the result, and choose horizons after the first event.

Logical predictor columns

partial.rfsrc does not handle logical predictor columns correctly in survival forests (randomForestSRC <= 3.5.1). If your training data contains binary 0/1 columns, convert them to factor rather than logical before fitting the model.

Note

For survival forests this function defaults to partial.type = "surv", so yhat is a survival probability on $[0, 1]$. gg_partial wraps randomForestSRC::plot.variable instead, which defaults to surv.type = "mort" and therefore yields mortality (an expected event count). The two entry points consequently report different quantities by default; check the y-axis label, which each records on the returned object.

See Also

gg_partial, partial.rfsrc, get.partial.plot.data

Examples

## ------------------------------------------------------------
##
## regression
##
## ------------------------------------------------------------

airq.obj <- randomForestSRC::rfsrc(Ozone ~ ., data = airquality, ntree = 100)

## partial effect for wind
prt_dta <- gg_partial_rfsrc(airq.obj,
                       xvar.names = c("Wind"))

Description

varPro::partialpro returns one list, with continuous and categorical predictors mixed together. This function splits that list into two tidy data frames, one for each kind, and resolves the y-axis label the plot method will use.

Deprecated. gg_partialpro() has been superseded by gg_partial_varpro() and is now a thin alias for it. It will be removed in a future release; use gg_partial_varpro() directly.

Usage

gg_partial_varpro(
  part_dta = NULL,
  object = NULL,
  scale = c("auto", "prob", "prob_typical", "odds", "logodds", "rmst", "surv",
    "mortality", "chf"),
  time = NULL,
  nvars = NULL,
  cat_limit = 10,
  model = NULL,
  ...
)

gg_partialpro(
  part_dta,
  object = NULL,
  scale = c("auto", "prob", "prob_typical", "odds", "logodds", "rmst", "surv",
    "mortality", "chf"),
  time = NULL,
  nvars = NULL,
  cat_limit = 10,
  model = NULL,
  ...
)

Arguments

Details

Which variables you actually get: varpro screens twice before anything reaches partial dependence. The split-weight screen decides which predictors are worth guiding the trees with, and what survives it lands in object$xvar.names; varPro::get.topvars then ranks a shorter list out of that. So the design matrix, the reachable set, and the default list are three different sizes; a fit on 45 predictors might carry 26 in object$xvar.names and 15 in get.topvars. partialpro can only reach the middle one. How much gets screened off depends on the data and the fit, so check rather than assume: length(object$xvar.names) against ncol(object$x) tells you where you stand.

This bites when you bring a variable list in from somewhere else, say the top names off an rfsrc VIMP ranking. partialpro intersects your xvar.names with what it can reach and keeps the overlap, so a request for twelve variables can come back with ten. It is the intermittent kind of trap: a top-10 list may come back whole while a top-12 list loses two. partialpro warns and names what it dropped, but only after the isolation-forest work is done. A quick setdiff(my_names, object$xvar.names) answers the same question before you spend the computation.

For a complete view, fit with both screens off: varPro::varpro(..., sparse = FALSE, split.weight = FALSE). split.weight = FALSE is the one that lifts the ceiling. It puts every predictor in object$xvar.names, so partial dependence can reach them all, and it leaves a strong variable's curve where it was. sparse = FALSE does the smaller thing, deepening varPro::get.topvars so the reported ranking shows its tail rather than the screened top. Both are varPro's own arguments, and its defaults (both TRUE) go the other way, toward the screened set; keep the defaults when that sparser set is what you want. nvar is not the knob here; it only caps how much gets reported.

Which rows you actually get: varPro has no imputation. Every entry point opens by growing a one-node stump through randomForestSRC::rfsrc to settle the family and hand back cleaned data, and that call takes rfsrc's default na.action = "na.omit". Any case with a missing value, in a predictor or in the outcome, is deleted before the fit. What you are told about it depends on the varPro version. Before 3.2.2 the deletion is silent: no warning, no message, and na.action = "na.impute" passed to varPro::varpro lands in ..., never reaches the stump, and is discarded without remark. From 3.2.2, varpro warns with the input, omitted and retained counts, records them in object$model.info$observations, and stops with an error on na.action (or any other argument it does not recognise) instead of ignoring it. The rows are deleted either way; 3.2.2 makes the deletion visible.

The loss compounds across predictors rather than adding up. At 5% missing per column, independently, retention is $0.95^p$: 60% of rows at 10 predictors, 36% at 20, 8% at 50. On a wide clinical frame a little missingness everywhere can delete most of the cohort. Before 3.2.2 nothing in the fit records it: object$rf$n is the count after deletion and no original is kept, so neither you nor this package can recover the number from the object. On any version, check before you fit: nrow(dta) against sum(complete.cases(dta)).

Imputing first, without inventing outcomes: varPro::roughfix does mean and modal fill, randomForestSRC::impute does the forest-based job, and varPro's own help pages use the latter. Both fill every column they are handed, the outcome included, so where the outcome is itself missing they manufacture it, and the release rules are then fit partly to invented responses. Put the observed outcome back and drop the cases that never had one:

imp   <- randomForestSRC::impute(y ~ ., data = dta)
imp$y <- dta$y                    # restore the observed outcome
imp   <- imp[!is.na(imp$y), ]     # drop cases with no real outcome

That is von Hippel's impute-then-delete: keep the outcome in the imputation model, since leaving it out attenuates the associations varPro is looking for, then discard the cases whose outcome was imputed, which carry no information about the response. Only the deletion step depends on having missing outcomes; the imputation matters whenever a predictor is missing. Two cautions. Imputing with the outcome stamps its signal into the filled predictor cells, which is right for a single fit but crosses fold boundaries in varPro::cv.varpro; impute inside the folds if the selection error has to mean anything. And a completed frame is one dataset, not many, so the curves here carry no uncertainty from the imputation; read them as conditional on it.

nvars selects by importance, not by list position: when object is supplied, the variables in part_dta are first reordered by their rank in varPro::get.topvars(object), and nvars then keeps the top nvars of that ordering, not the first nvars elements as they arrived. get.topvars() typically ranks far fewer variables than part_dta contains; any variable it does not rank keeps its incoming order and is appended after the ranked block, so nothing is ever dropped from consideration, only reordered. The returned continuous$name / categorical$name column is a factor whose levels follow this same importance order (raw list order when no object was supplied), so facet panels and legend entries sort by importance rather than alphabetically.

Scale detection: with scale = "auto" and an object in hand, the scale resolves to "surv" for a survival forest, "prob" for a classification forest and "generic" for a regression forest. We do not read an RMST horizon $\tau$ off the varpro object. Before varPro 3.2.2 there is none to read; from 3.2.2 one is recorded only when varpro itself was fit with an rmst horizon. So RMST output requires you to pass scale = "rmst", time = tau explicitly.

RMST partial dependence (scale = "rmst"): varPro::partialpro has no time argument, so its default survival learner returns ensemble mortality at every horizon. Before varPro 3.2.2 a horizon passed through ... is silently dropped, and multi-horizon plots built that way differ only by Monte-Carlo noise, not by $\tau$; from 3.2.2 partialpro stops with an error instead. To get a genuine RMST($\tau$) curve, scale = "rmst" supplies partialpro a learner that returns $\mathrm{RMST}(\tau)=\int_0^\tau S(t)\,dt$ from the survival forest, so the curve actually depends on $\tau$. This path recomputes from object, so it needs object (a survival fit) with part_dta = NULL; a precomputed part_dta can only be relabeled, and gg_partial_varpro warns when you try. A $\tau$ beyond the model's largest event time is truncated there (with a warning), since $S(t)$ cannot be extrapolated.

Classification scale (scale = "prob"/"odds"/"logodds"): varPro::partialpro returns classification effects as log-odds of the target class. scale = "prob" (the classification default) back-transforms to probability $P(Y = \mathrm{target})$, "odds" to the odds, and "logodds" keeps the raw scale. The back-transform is applied per observation before averaging, so the curve is the mean predicted probability, not the probability of the mean log-odds.

Restoring each subject's level (scale = "prob"): partialpro fits each subject's curve separately but returns every row at the cohort-mean intercept, keeping only the subject's own slope. On the log-odds scale that sets every subject to the average log-odds, and averaging per-subject probabilities then no longer gives the expected proportion. So for "prob", gg_partial_varpro shifts each subject's curve to pass through that subject's own out-of-bag log-odds at its observed value of the variable, and then back-transforms and averages. The shape partialpro fitted is unchanged. This needs object, the classification fit, and is recorded as anchored in the provenance. Without object, or when ... passes partialpro a custom learner or newdata, the curve is left as partialpro returned it (with a warning when object is missing). A precomputed part_dta is anchored whenever object is supplied, which assumes it came from partialpro(object) with its default learner; a part_dta built with your own learner should be passed with scale = "logodds" instead. A subject that was never out of bag takes its in-bag prediction as its anchor. Binary variables are never shifted, because partialpro already returns per-subject levels for them. The causal contrast is shown only on "logodds" (see plot.gg_partial_varpro).

Two probability scales, and why they disagree (scale = "prob" vs "prob_typical"): partialpro returns a matrix of per-subject log-odds, one row per observation and one column per grid point. Collapsing it to a curve takes an average and a back-transform, and the order of those two steps is a modelling choice, not a detail:

$\mathrm{prob}(x) = \frac{1}{n}\sum_i \mathrm{logit}^{-1}(z_i(x)) \qquad \mathrm{prob\_typical}(x) = \mathrm{logit}^{-1}\!\left( \frac{1}{n}\sum_i z_i(x)\right)$

"prob" (the classification default) transforms per observation and then averages, so the curve is the mean predicted probability, the expected proportion of this cohort, and the standard partial dependence quantity on a probability scale. "prob_typical" averages on the log-odds scale and then transforms once, giving the probability for a subject sitting at the mean log-odds.

These are different estimands and they do not agree. $\mathrm{logit}^{-1}$ is concave above zero and convex below it, so by Jensen's inequality "prob" is pulled toward $0.5$ relative to "prob_typical", at both ends of the curve. The gap widens with the spread of per-subject log-odds, and on a heterogeneous cohort it is not small. In a simulation where a second variable spreads the subjects' log-odds to an SD near 4, one point reads $0.13$ under "prob_typical" and $0.35$ under "prob", against a true partial dependence of $0.38$. Before the level restoration described above, "prob" read $0.14$ there: with every subject at the mean log-odds, the two scales nearly coincide. "prob_typical" uses partialpro's values as returned.

Which to report is a question about the claim, not about the code. If the sentence is "what fraction of these patients would wean", that is "prob". If it is "what would we predict for a typical patient", that is "prob_typical", while remembering that the mean-log-odds subject need not resemble anyone in the data. A figure captioned as a percentage of patients wants "prob".

The distinction applies only to the continuous frame. The categorical frame keeps its values unaveraged so the plot method can draw boxplots, so there is no averaging order to choose and the two scales return the same numbers there.

Survival probability (scale = "surv"): scale = "surv" (the survival default) computes $S(\tau \mid x)$ through partialpro (the same UVT engine as mortality and RMST), bounded in $[0, 1]$. partialpro smooths each case's $S(\tau)$ with an unbounded polynomial, so where survival is near 1 the averaged curve can run slightly past it; the curve is clamped to $[0, 1]$ after averaging, and the categorical frame per value. The curve is still the average of the smoothed per-case $S(\tau)$, with only impossible values changed. When time is not supplied, $\tau$ defaults to the median follow-up time of the fit, a data-driven horizon that is always in the model's own time units, so it cannot be mis-specified the way a hand-typed $\tau$ can. The resolved $\tau$ is reported in a message and the axis label; pass time = tau to choose another. scale = "mortality" keeps the unbounded ensemble-mortality score as an explicit opt-in.

Ensemble mortality (scale = "mortality"): here the y-axis is ensemble mortality, the expected number of events a subject would see if they were exposed to the study-average cumulative hazard. It is the same quantity as the rfsrc predicted value for survival forests (Ishwaran, Kogalur, Blackstone & Lauer, 2008 https://doi.org/10.1214/08-AOAS169). This is an unbounded relative-risk score, not a survival probability and not $1 - S(t)$; don't read it as one. For a bounded, time-anchored survival summary, use scale = "rmst", time = tau (restricted mean survival time, in the time units of the outcome) or scale = "surv" / "chf".

Arguments are documented on gg_partial_varpro; this alias shares its formals and forwards every argument unchanged.

Value

A named list of class "gg_partial_varpro" with elements:

continuous

data.frame with columns variable, parametric, nonparametric, causal, name (and optionally model).

categorical

data.frame with the same columns, one row per observation per category level.

A "provenance" attribute carries source, family, ntree, n, scale, rmst_tau, xvar.names, and path.

A gg_partial_varpro object (see gg_partial_varpro).

What partialpro is doing

A partial dependence curve answers the question, "if I hold a single variable at a grid of values and average out everything else, how does the model's prediction move?" That is the same question rfsrc partial dependence answers. What varPro::partialpro adds is two wrinkles that are worth understanding before you read the curves.

First, partialpro filters the partial grid through an isolation forest (Unlimited Virtual Twins, or UVT) so that unlikely combinations of the focal variable with the rest of the data are downweighted. The rfsrc version, by contrast, averages over the full marginal grid regardless of plausibility. So when a covariate is highly correlated with others, the two methods can disagree, and partialpro's curve is the one restricted to the data manifold.

That isolation forest is worth one note. partialpro grows it with varPro::isopro, whose method defaults to "unsupv", a forest with no outcome. randomForestSRC then passes a zero-length yvar.wt into its native code and decrements that pointer (entry.c:184), which is undefined behaviour and is reported by UBSAN builds. It is benign in practice (the pointer is formed but never dereferenced), and it is an upstream issue rather than one this package can fix (ggRandomForests is pure R); it is tracked as kogalur/randomForestSRC issue #477. Passing method = "rnd" through to partialpro avoids the unsupervised grow entirely, which is what this package's own tests and examples do.

Second, partialpro fits a local polynomial model to the predicted values rather than just plotting their mean. That gives three parallel curves per variable, stored as yhat.par, yhat.nonpar, and yhat.causal, which the plot method overlays so you can see whether a smooth parametric story and the raw forest predictions are telling you the same thing.

Interpretation of the y-axis depends on the outcome (per varPro::partialpro): response scale for regression, log-odds of the target class for classification, and either ensemble mortality (default) or RMST (if the original varpro call set rmst) for survival.

What's in the output

We split partialpro's mixed list into two tidy data frames so the plot method does not have to. A variable with more than cat_limit distinct grid points goes into $continuous, one row per grid point with the column means of yhat.par, yhat.nonpar, and yhat.causal stored as parametric, nonparametric, and causal. A variable at or below cat_limit goes into $categorical, one row per observation per category level, carrying the same three columns unaveraged so the plot method can draw boxplots. Path C (scale %in% c("surv","chf")) takes a different route: we hand the underlying rfsrc forest to gg_partial_rfsrc so you get a survival-probability or cumulative-hazard curve on the usual rfsrc scale instead.

What you use this for

• read the marginal shape of a relationship the varpro model found important: monotone, threshold, U-shape, flat;

• compare the three partialpro estimators on the same variable and flag the ones where parametric and nonparametric disagree, those are the candidates for closer inspection;

• report a survival partial dependence on the probability or cumulative-hazard scale (scale = "surv" or "chf") rather than the unbounded mortality scale.

A varpro partial dependence curve is a description of the model, not a causal effect. The causal column is varpro's local estimator, not a structural causal claim about the data-generating process.

Factors, bands and small data

varPro one-hot encodes a factor, so a factor stage becomes the terms stage1, stage2, ..., and partialpro() looks terms up by those names (pass xvar.names = "stage4", not "stage"). Each level is partialled on its own, every case set to 0 and to 1 on that level. partialpro() leaves the factor's other level columns as they were, so for a factor with three or more levels some of those cases hold two levels at once; this is upstream behaviour in varPro::partialpro(). Two-level factors, stored as a single 0/1 column, are not affected. To check a multi-level factor, use gg_partial_rfsrc, which sets the factor itself to each level.

The spread across cases is not a confidence interval. It shows how much cases differ under one fitted forest; varPro's is narrower still because partialpro() gives every case the cohort-mean intercept. To see how uncertain the curve is, refit on subsamples. The varpro-small-data vignette shows both.

References

von Hippel PT (2007). Regression with missing Ys: An improved strategy for analyzing multiply imputed data. Sociological Methodology, 37(1), 83–117. doi:10.1111/j.1467-9531.2007.00180.x.

Ishwaran H, Kogalur UB, Blackstone EH, Lauer MS (2008). Random survival forests. The Annals of Applied Statistics, 2(3), 841–860. doi:10.1214/08-AOAS169.

Ishwaran H, Blackstone EH (2025). Harnessing the power of virtual (digital) twins: Graphical causal tools for understanding patient and hospital differences. Computational and Structural Biotechnology Journal, 28, 312.

See Also

plot.gg_partial_varpro, gg_varpro, gg_vimp, gg_partialpro (deprecated), gg_partial_rfsrc, varpro_feature_names

gg_partial_varpro

Examples

set.seed(42)
n_obs <- 30; n_pts <- 15
mock_data <- list(
  age = list(
    xvirtual    = seq(30, 80, length.out = n_pts),
    xorg        = sample(seq(30, 80, by = 5), n_obs, replace = TRUE),
    yhat.par    = matrix(rnorm(n_obs * n_pts), nrow = n_obs),
    yhat.nonpar = matrix(rnorm(n_obs * n_pts), nrow = n_obs),
    yhat.causal = matrix(rnorm(n_obs * n_pts), nrow = n_obs)
  ),
  sex = list(
    xvirtual    = c(0, 1),
    xorg        = sample(c(0, 1), n_obs, replace = TRUE),
    yhat.par    = matrix(rnorm(n_obs * 2), nrow = n_obs),
    yhat.nonpar = matrix(rnorm(n_obs * 2), nrow = n_obs),
    yhat.causal = matrix(rnorm(n_obs * 2), nrow = n_obs)
  )
)
## The two probability scales differ by the ORDER of averaging and
## back-transform, and disagree whenever subjects are heterogeneous.
## Mock data has no fit, so "prob" cannot restore subject levels and warns.
pa <- suppressWarnings(gg_partial_varpro(mock_data, scale = "prob"))
pt <- gg_partial_varpro(mock_data, scale = "prob_typical")
head(data.frame(prob = pa$continuous$parametric,
                prob_typical = pt$continuous$parametric))

result <- gg_partial_varpro(mock_data, scale = "logodds")
head(result$continuous)
head(result$categorical)


## The object-driven path: hand gg_partial_varpro() the varpro fit and let
## it call partialpro() for you.  This is where the reachability ceiling
## shows up (see Details).
set.seed(42)
vp <- varPro::varpro(mpg ~ ., data = mtcars, ntree = 50)

## Three different sizes.  partialpro() can only reach the second.
ncol(vp$x)                    # predictors in the data
length(vp$xvar.names)         # what the fit reaches
length(varPro::get.topvars(vp))   # the default when xvar.names is absent

## Say these came from an rfsrc VIMP ranking.  Check what the fit cannot
## reach before you spend the computation -- this is the habit worth having.
wanted <- c("wt", "hp", "qsec", "vs")
setdiff(wanted, vp$xvar.names)

## Ask anyway and partialpro() warns, naming what it dropped.
## (method = "rnd" is passed through to partialpro(); see the note on
## isolation-forest method in Details.)
pd <- gg_partial_varpro(object = vp, xvar.names = wanted,
                        method = "rnd")

## Refitting without the split-weight screen reaches every predictor.
vp_all <- varPro::varpro(mpg ~ ., data = mtcars, ntree = 50,
                         split.weight = FALSE)
length(vp_all$xvar.names)
setdiff(wanted, vp_all$xvar.names)   # empty; nothing to drop

Description

Every tree in a random forest makes its own prediction, and the forest's "ensemble" prediction is the average across all trees. The out-of-bag (OOB) variant averages only over the trees that did not include a given observation in their bootstrap sample, a built-in cross-validation estimate that requires no held-out test set. gg_rfsrc pulls those ensemble predictions out of the fitted forest and arranges them for plotting with plot.gg_rfsrc.

Usage

## S3 method for class 'rfsrc'
gg_rfsrc(object, oob = TRUE, by, ...)

Arguments

Details

The structure of the returned data depends on the forest family. For a regression forest you get a scatter of OOB predicted values against the observed response. For classification you get predicted class probabilities alongside the observed class label. For a survival forest you get the ensemble survival function (or cumulative hazard, or mortality, controlled by surv_type) at each unique event time (one curve per observation), which together trace the range of predicted risk in the cohort. Pass conf.int to add pointwise bootstrap confidence bands around the mean survival curve, or by to stratify all of the above by a predictor group.

For survival forests, use the surv_type argument ("surv", "chf", or "mortality") to select the predicted quantity. Bootstrap confidence bands are requested by passing conf.int (e.g. conf.int = 0.95); the number of resamples is controlled by bs.sample.

Value

A gg_rfsrc object (a classed data.frame) whose structure depends on the forest family:

regression

Columns yhat and the response name; optionally a group column when by is supplied.

classification

One column per class with predicted probabilities; a y column with observed class labels; optionally group.

survival (no CI / grouping)

Long-format with columns variable (event time), value (survival probability), obs_id, and event.

survival (with conf.int or by)

Wide-format with pointwise bootstrap CI columns (lower, upper, median, mean) per time point; a group column when by is supplied.

The object carries class attributes for the forest family so that plot.gg_rfsrc dispatches correctly.

See Also

plot.gg_rfsrc, rfsrc, gg_survival

Examples

## ------------------------------------------------------------
## classification example (small, runs on CRAN)
## ------------------------------------------------------------
## -------- iris data
set.seed(42)
rfsrc_iris <- randomForestSRC::rfsrc(Species ~ ., data = iris, ntree = 50)
gg_dta <- gg_rfsrc(rfsrc_iris)
plot(gg_dta)


## ------------------------------------------------------------
## Additional regression / survival examples are guarded with
## \donttest because the cumulative example time exceeds the
## 10-second CRAN budget. Run locally with `R CMD check --run-donttest`
## (or `devtools::check(run_dont_test = TRUE)`) to exercise them.
## ------------------------------------------------------------

## -------- air quality data (regression)
rfsrc_airq <- randomForestSRC::rfsrc(Ozone ~ ., data = airquality,
                    na.action = "na.impute", ntree = 50)
plot(gg_rfsrc(rfsrc_airq))

## -------- Boston data (rfsrc + randomForest)
if (requireNamespace("MASS", quietly = TRUE)) {
  data(Boston, package = "MASS")
  Boston$chas <- as.logical(Boston$chas)
  rfsrc_boston <- randomForestSRC::rfsrc(medv ~ ., data = Boston, ntree = 50,
                        forest = TRUE, importance = TRUE,
                        tree.err = TRUE, save.memory = TRUE)
  plot(gg_rfsrc(rfsrc_boston))

  rf_boston <- randomForest::randomForest(medv ~ ., data = Boston,
                                          ntree = 50)
  plot(gg_rfsrc(rf_boston))
}

## -------- mtcars data
rfsrc_mtcars <- randomForestSRC::rfsrc(mpg ~ ., data = mtcars, ntree = 50)
plot(gg_rfsrc(rfsrc_mtcars))

## -------- veteran data (survival; with CI and group-by)
data(veteran, package = "randomForestSRC")
rfsrc_veteran <- randomForestSRC::rfsrc(Surv(time, status) ~ ., data = veteran,
                       ntree = 50)
plot(gg_rfsrc(rfsrc_veteran))
plot(gg_rfsrc(rfsrc_veteran, conf.int = .95))
plot(gg_rfsrc(rfsrc_veteran, by = "trt"))

Description

Extracts case-specific ensemble hazard and cumulative-hazard estimates from a fitted randomForestRHF::rhf() object into a tidy long data frame, one row per (case, time) pair on the forest's time.interest grid.

Usage

gg_rhf(object, ...)

## S3 method for class 'rhf'
gg_rhf(object, source = c("oob", "inbag"), ...)

Arguments

Value

A data.frame of class c("gg_rhf", "data.frame") with columns id, time, hazard, chf, source, an integer ntime attribute (number of grid points), and a provenance attribute.

The frame always has ntime rows per case. Both hazard and chf may be NA, and they are masked on different rules. From randomForestRHF 2.0.0 the pointwise hazard is defined only where a grid point falls inside one of the case's supplied ⁠(start, stop]⁠ intervals, and is NA in gaps and after the final stop time. From 2.0.3 the cumulative hazard is NA after each case's final stop, but stays flat through an internal gap rather than going missing there. So on a fit whose cases each have a single interval the two masks look the same, and with time-dependent covariates they do not. Earlier versions carried chf across the whole grid. The values are passed through unchanged; use na.rm = TRUE when summarising either column, or see plot.gg_rhf(), which drops the masked cells before drawing.

References

Ishwaran H, Hsich EM, Kogalur UB, Lee DKK (2026). Random Hazard Forests. arXiv:2608.21597. doi:10.48550/arXiv.2608.21597.

Ishwaran H, Kogalur UB (2026). randomForestRHF: Random Hazard Forests. R package version 2.0.3. https://CRAN.R-project.org/package=randomForestRHF.

See Also

plot.gg_rhf(), randomForestRHF::rhf()

Examples

if (requireNamespace("randomForestRHF", quietly = TRUE)) {
  data(pbc, package = "randomForestSRC")
  d <- randomForestRHF::convert.counting(
    survival::Surv(days, status) ~ ., na.omit(pbc))
  o <- randomForestRHF::rhf("Surv(id, start, stop, event) ~ .", d, ntree = 30)
  gg <- gg_rhf(o)
  plot(gg, idx = c(1, 5, 10))
}

Description

Extracts the variable-priority matrix from randomForestRHF::importance.rhf() into a tidy data frame, one row per variable and time window. The score measures how much the fitted integrated hazard changes when rules involving a variable are released.

Usage

gg_rhf_importance(object, ...)

## S3 method for class 'rhf'
gg_rhf_importance(
  object,
  importance_fit = NULL,
  cache = NULL,
  time.index = NULL,
  ...
)

Arguments

Details

Variable priority is time-localized. A large value means that releasing rules involving that variable changed the log integrated hazard more in that window. It is a ranking score, not a z-score, and this function does not apply a significance cutoff.

Calculating the upstream result can be expensive. For an analysis you will revisit, calculate importance_fit once and supply it here. The extractor accepts cache, time.index, and additional calculation arguments only when importance_fit is NULL.

Variables are ordered for plotting by their 90th percentile (q90) priority across windows. This changes the factor levels, but does not change the upstream row order or priority values.

Value

A data.frame of class c("gg_rhf_importance", "data.frame") with columns:

variable

Variable name, ordered by q90 priority for plotting.

time_window

Upstream time-window label.

time

Evaluation time at the end of the window.

time_index

Index of time on the RHF time grid.

start, stop, midpoint

Window boundaries and midpoint.

n_risk

Number of observations at risk.

n_rules

Number of rules contributing to the window.

priority

RHF variable-priority score.

A provenance attribute records the source forest, upstream settings, whether importance_fit was supplied, and the installed randomForestRHF version.

References

Ishwaran H, Hsich EM, Kogalur UB, Lee DKK (2026). Random Hazard Forests. arXiv:2608.21597. doi:10.48550/arXiv.2608.21597.

Ishwaran H, Kogalur UB (2026). randomForestRHF: Random Hazard Forests. R package version 2.0.3. https://CRAN.R-project.org/package=randomForestRHF.

See Also

plot.gg_rhf_importance(), randomForestRHF::importance.rhf()

Examples

if (requireNamespace("randomForestRHF", quietly = TRUE)) {
  data(pbc, package = "randomForestSRC")
  d <- randomForestRHF::convert.counting(
    survival::Surv(days, status) ~ ., na.omit(pbc))
  o <- randomForestRHF::rhf(
    "Surv(id, start, stop, event) ~ .", d, ntree = 30)

  priority_fit <- randomForestRHF::importance.rhf(o)
  priority <- gg_rhf_importance(o, importance_fit = priority_fit)
  plot(priority)
}

Description

A classifier does not hand you a class; it hands you a predicted probability, and you pick a threshold. Slide that threshold from 0 to 1 and the trade-off between catching the positives and crying wolf shifts the whole way. The ROC curve traces that trade-off. For one class of a classification rfsrc or randomForest forest, gg_roc walks every threshold and records sensitivity (the true positive rate) against specificity (1 minus the false positive rate).

Usage

## S3 method for class 'rfsrc'
gg_roc(object, which_outcome, oob = TRUE, per_class = FALSE, ...)

Arguments

Value

A gg_roc data.frame, one row per unique prediction threshold, with columns:

sens

Sensitivity (true positive rate) at the threshold.

spec

Specificity (true negative rate) at the threshold.

pct

The probability threshold used for that row.

Pass it to calc_auc for the area under the curve.

Defaults across entry points

plot(gg_roc(x)) and a direct plot.gg_roc(x) call on a raw multi-class forest do not draw the same figure: gg_roc(x) returns a single curve, while plot.gg_roc given a raw forest overlays one curve per class. Plain plot(x) on a forest reaches neither, as it dispatches to the forest's own method in randomForestSRC or randomForest. Pass which_outcome explicitly when the distinction matters; issue #72 tracks reconciling the two.

See Also

plot.gg_roc, calc_roc, calc_auc, rfsrc, randomForest

Examples

## ------------------------------------------------------------
## classification example
## ------------------------------------------------------------
## -------- iris data
rfsrc_iris <- randomForestSRC::rfsrc(Species ~ ., data = iris, ntree = 100)

# ROC for setosa
gg_dta <- gg_roc(rfsrc_iris, which_outcome = 1)
plot(gg_dta)

# ROC for versicolor
gg_dta <- gg_roc(rfsrc_iris, which_outcome = 2)
plot(gg_dta)

# ROC for virginica
gg_dta <- gg_roc(rfsrc_iris, which_outcome = 3)
plot(gg_dta)

## -------- iris data
rf_iris <- randomForest::randomForest(Species ~ ., data = iris)

# ROC for setosa
gg_dta <- gg_roc(rf_iris, which_outcome = 1)
plot(gg_dta)

# ROC for versicolor
gg_dta <- gg_roc(rf_iris, which_outcome = 2)
plot(gg_dta)

# ROC for virginica
gg_dta <- gg_roc(rf_iris, which_outcome = 3)
plot(gg_dta)

Description

Tidy wrapper around varPro::sdependent() for a uvarpro object. Where gg_udependent() draws the cross-variable dependency graph, gg_sdependent() surfaces sdependent()'s signal-variable detection: a ranked table of the per-variable signal score and graph degree, with the variables flagged as "signal" (those whose dependency structure clears the detection threshold).

Usage

gg_sdependent(
  object,
  ...,
  threshold = 0.25,
  q.signal = 0.75,
  directed = TRUE,
  min.degree = NULL,
  beta_fit = NULL
)

Arguments

Details

sdependent() runs on the varPro::get.beta.entropy() lasso-coefficient matrix and returns, with plot = FALSE, a list of imp.score (per-variable signal score), degree (node degree in the dependency graph), and signal.vars (the detected signal set). This wrapper tidies that into one row per candidate variable, ranked by imp.score. Because the entropy matrix is the expensive part, beta_fit accepts a precomputed varPro::get.beta.entropy() matrix (shared with gg_beta_uvarpro() and gg_udependent()).

Value

A gg_sdependent object (a data.frame), one row per candidate variable, most-signal first, with columns:

variable

factor; levels reversed so the top variable lands at the top after coord_flip().

imp_score

sdependent() per-variable signal score.

degree

node degree in the dependency graph.

signal

logical; variable is in sdependent()$signal.vars.

The provenance attribute records source, family ("unsupv"), threshold, q.signal, directed, n_signal, and n_var.

See Also

gg_udependent() (the dependency graph), gg_beta_uvarpro() (lasso importance), varPro::sdependent(), varPro::uvarpro().

Examples

if (requireNamespace("varPro", quietly = TRUE)) {
  set.seed(1)
  o <- varPro::uvarpro(mtcars, ntree = 50)
  gg <- gg_sdependent(o)
  plot(gg)
}

Description

gg_shap computes SHAP values for a rfsrc or randomForest regression or classification forest by wrapping kernelshap, and reshapes them into a tidy data set with one row per (observation, variable).

Usage

gg_shap(object, newdata, bg_n = 50, which.class = 1, ...)

Arguments

Value

A gg_shap object: a data.frame with columns id (observation index), vars (variable name, an ordered factor ranked by mean absolute SHAP), shap (the signed SHAP contribution), value (numeric feature value, NA for categorical features), and value_label (feature value as character). The background-sample mean prediction is stored in the "baseline" attribute.

See Also

kernelshap

Examples

if (requireNamespace("kernelshap", quietly = TRUE)) {
  rf <- randomForestSRC::rfsrc(Ozone ~ ., data = na.omit(airquality),
                               ntree = 50)
  gg_dta <- gg_shap(rf, bg_n = 20)
}

Description

Nonparametric survival estimates.

Usage

gg_survival(
  object = NULL,
  interval = NULL,
  censor = NULL,
  by = NULL,
  data = NULL,
  type = c("kaplan", "nelson"),
  ...
)

## S3 method for class 'rfsrc'
gg_survival(
  object,
  interval = NULL,
  censor = NULL,
  by = NULL,
  data = NULL,
  type = c("kaplan", "nelson"),
  ...
)

## Default S3 method:
gg_survival(
  object = NULL,
  interval = NULL,
  censor = NULL,
  by = NULL,
  data = NULL,
  type = c("kaplan", "nelson"),
  ...
)

Arguments

Details

Placing the forest's ensemble survival curve over the marginal Kaplan-Meier baseline is a calibration check, not a test of signal. The ensemble curve averages the per-subject curves, so it should track the Kaplan-Meier estimate whether or not the predictors carry information; a forest grown on pure noise tracks it about as closely as one grown on real predictors. A clear gap points to miscalibration, most often at late times when few subjects remain at risk. To see what the predictors carry, stratify with by and compare the groups. gg_survival computes the nonparametric baseline (the Kaplan-Meier or Nelson-Aalen estimate) so you can place it on the same canvas as the forest predictions from gg_rfsrc.

gg_survival is an S3 generic that dispatches on the class of its first argument:

rfsrc

Extracts the outcome columns from the fitted forest's $yvar slot (time in column 1, event indicator in column 2) and delegates to kaplan. Use by to stratify on a predictor from $xvar: you get one Kaplan-Meier curve per group, ready to compare against the forest's group-specific ensemble curves.

default

Accepts raw survival columns directly via interval, censor, and data. Delegates to kaplan (the default) or nelson depending on type.

Value

A gg_survival data.frame with columns time, surv, cum_haz, lower, upper, n.risk, and optionally groups when by is supplied.

Note

Survival estimation is randomForestSRC-only; randomForest has no survival forest, so no randomForest method exists.

See Also

kaplan nelson

plot.gg_survival

Examples

## -------- pbc data (default method, raw data columns)
data(pbc, package = "randomForestSRC")
pbc$time <- pbc$days / 364.25

gg_dta <- gg_survival(interval = "time", censor = "status", data = pbc)
plot(gg_dta, error = "none")

# Stratified
gg_dta <- gg_survival(
  interval = "time", censor = "status",
  data = pbc, by = "treatment"
)
plot(gg_dta)

Description

Extracts an already calculated tree-size tuning path from randomForestRHF::tune.treesize.rhf() into a data frame for inspection and plotting. The expensive step is upstream tuning. Calculate and retain that result once, then supply it to gg_tune_rhf() when you need the saved search path or its plot. gg_tune_rhf() only prepares that path; it never tunes a forest.

Usage

gg_tune_rhf(tune_fit, ...)

## S3 method for class 'tune.treesize.rhf'
gg_tune_rhf(tune_fit, ...)

Arguments

Details

The returned path preserves the row order in tune_fit$path. Its columns are treesize (evaluated forest size), metric ("OOB risk" or "OOB iAUC"), value (the observed metric), se (the supplied bootstrap iAUC standard error, or NA_real_), and selected (whether that size is the upstream best.size). Upstream tuning minimizes OOB risk or maximizes OOB iAUC.

Provenance is stored in the provenance attribute: best_size is the selected tree size; best_err is OOB risk for risk tuning and 1 - iAUC for iAUC tuning; perf identifies the criterion; method is the upstream search method; bounds gives its tree-size range; n_evaluations counts the evaluated sizes; and randomForestRHF_version records the installed upstream package version. The optional fitted forest is not copied into the tidy result.

Value

A data.frame with class c("gg_tune_rhf", "data.frame") and columns treesize, metric, value, se, and selected. The provenance attribute contains the upstream settings described in gg_tune_rhf.

References

Ishwaran H, Hsich EM, Kogalur UB, Lee DKK (2026). Random Hazard Forests. arXiv:2608.21597. doi:10.48550/arXiv.2608.21597.

Ishwaran H, Kogalur UB (2026). randomForestRHF: Random Hazard Forests. R package version 2.0.3. https://CRAN.R-project.org/package=randomForestRHF.

See Also

plot(), randomForestRHF::tune.treesize.rhf()

Examples

if (requireNamespace("randomForestRHF", quietly = TRUE)) {
  ## Calculate this expensive result once and retain it for reuse.
  simulated <- randomForestRHF::hazard.simulation(1, n = 100, nrecords = 3)
  tune_fit <- randomForestRHF::tune.iAUC.rhf(
    "Surv(id, start, stop, event) ~ .",
    simulated$dta,
    ntree = 12L,
    lower = 2L,
    upper = 5L,
    verbose = FALSE,
    forest = FALSE
  )
  tuning <- gg_tune_rhf(tune_fit)
  plot(tuning)
}

Description

A uvarpro fit records how strongly each variable depends on the others. This function pulls those cross-variable dependency scores from the fit with get.beta.entropy and sdependent, and returns them as a tidy list that plot.gg_udependent can draw as a network.

Usage

gg_udependent(
  object,
  threshold = 0.25,
  q.signal = 0.75,
  directed = TRUE,
  min.degree = NULL,
  ...,
  pre.filter = TRUE,
  beta.mat = NULL
)

Arguments

Value

A named list of class "gg_udependent" with elements:

$edges

Data frame: variable_from, variable_to, weight (raw cross-importance value).

$nodes

Data frame: variable (factor, levels by descending degree), degree (integer; out-degree when directed = TRUE, total degree when directed = FALSE), selected (logical, TRUE if in sdependent's signal set).

$graph

igraph object. NULL if no dependencies detected.

A "provenance" attribute carries threshold, q.signal, directed, min.degree, pre.filter (NA when beta.mat was supplied, since it was not used), beta.mat.supplied, xvar.names, and n.

What cross-variable dependency is doing

UVarPro (Zhou, Lu and Ishwaran, 2026) extends the varpro framework to the unsupervised setting: grow a forest without a response, then use the same region-release contrasts varpro uses for supervised importance to ask, "which variables explain the structure in the data?" The lasso-driven variant frames each region-release contrast as a classification task (does an observation belong to the region or to its release?) and fits a lasso logistic regression with the other variables as predictors. The coefficient on variable $j$ in the model for variable $i$'s region-release contrast is the entry $I[i, j]$ of the matrix varPro::get.beta.entropy() returns.

Read that entry as "how much does knowing $j$ help separate $i$'s region from its release". A large $I[i, j]$ says $j$ carries information about the structure varpro picked up in $i$. varPro::sdependent thresholds that matrix at a user-chosen cut and returns the set of "signal" variables: the nodes with high enough out-degree to be worth keeping. We pass the threshold through to sdependent and use the same matrix to weight the edges of the resulting graph.

The graph is directed by default because $I[i, j]$ and $I[j, i]$ are separate lasso coefficients and need not agree; setting directed = FALSE collapses each pair by taking the larger of the two, which is appropriate when you only want to see that two variables are dependent, not which way the dependency reads.

Mixed-type data

The matrix get.beta.entropy() returns is usually not square. Its rows are the variables that produced region releases; its columns are every one-hot-encoded predictor column, so factor levels that never released, and variables the pre-filter dropped, appear only as columns. We pad it to a square matrix over the union of the names, filling with 0 (no release, no dependency). We also zero the entries between levels of the same factor (cyl4 and cyl8, say): the levels are mutually exclusive by construction, so that edge says nothing about the data.

What's in the output

$edges has one row per surviving edge with the raw weight I[i, j] (or, for undirected graphs, the max of the two directions). $nodes has one row per surviving variable with its degree (out-degree for directed, total degree for undirected) and a selected flag for membership in the sdependent signal set. $graph is the same information packaged as an igraph object, with weight, degree, and selected attached so plot.gg_udependent can render it without recomputing anything.

What you use this for

• screen a wide unsupervised dataset for the small set of variables UVarPro thinks are carrying the signal: the nodes with high degree, or those flagged selected = TRUE;

• spot clusters of mutually dependent variables (hubs and the spokes around them) that may be measuring the same underlying construct;

• compare two datasets, or two preprocessing pipelines, by looking at how their dependency graphs change.

An edge in this graph is a statistical dependency in the unsupervised decomposition of the data. It is not a causal arrow. A high $I[i, j]$ says $j$ predicts $i$'s region membership, not that $j$ causes $i$.

References

Zhou, L., Lu, M. and Ishwaran, H. (2026). Variable priority for unsupervised variable selection. Pattern Recognition, 172:112727.

See Also

plot.gg_udependent

Examples

set.seed(42)
uv <- varPro::uvarpro(iris[, -5], ntree = 50)
gg <- gg_udependent(uv)
print(gg)

Description

plot.variable generates a data.frame containing the marginal variable dependence or the partial variable dependence. The gg_variable function creates a data.frame of containing the full set of covariate data (predictor variables) and the predicted response for each observation. Marginal dependence figures are created using the plot.gg_variable function. A randomForest fit does not keep the model frame, so for those objects gg_variable rebuilds it from the stored call. That lets the same predictors be paired with the in-sample predictions.

A few optional arguments tune the extraction: time (one survival time, or a vector of them), time_labels (labels for multiple survival horizons), and oob, which switches between out-of-bag and in-bag predictions when the forest carries both.

Usage

gg_variable(object, ...)

Arguments

Details

The marginal variable dependence is determined by comparing relation between the predicted response from the randomForest and a covariate of interest.

The gg_variable function operates on a rfsrc object, the output from the plot.variable function, or on a fitted randomForest object via the formula interface.

Value

A gg_variable object: a data.frame pairing every training predictor column with the OOB (or in-bag) predicted response. For survival forests, each requested time horizon adds a column named by time_labels. The object carries a "family" class attribute ("regr", "class", or "surv") that plot.gg_variable uses for dispatch.

See Also

plot.gg_variable, plot.variable

Examples

## ------------------------------------------------------------
## classification (small, runs on CRAN)
## ------------------------------------------------------------
## -------- iris data
set.seed(42)
rfsrc_iris <- randomForestSRC::rfsrc(Species ~ ., data = iris, ntree = 50)

gg_dta <- gg_variable(rfsrc_iris)
plot(gg_dta, xvar = "Sepal.Width")


## ------------------------------------------------------------
## Additional classification / regression / survival examples are
## guarded with \donttest because the cumulative example time exceeds
## the 10-second CRAN budget. Run locally with `R CMD check
## --run-donttest` (or `devtools::check(run_dont_test = TRUE)`) to
## exercise them.
## ------------------------------------------------------------
plot(gg_dta, xvar = "Sepal.Length")
plot(gg_dta, xvar = rfsrc_iris$xvar.names, panel = TRUE)

## ------------------------------------------------------------
## regression
## ------------------------------------------------------------

## -------- air quality data
rfsrc_airq <- randomForestSRC::rfsrc(Ozone ~ ., data = airquality, ntree = 50)
gg_dta <- gg_variable(rfsrc_airq)

# an ordinal variable
gg_dta[, "Month"] <- factor(gg_dta[, "Month"])

plot(gg_dta, xvar = "Wind")
plot(gg_dta, xvar = "Temp")
plot(gg_dta, xvar = "Solar.R")
plot(gg_dta, xvar = c("Solar.R", "Wind", "Temp", "Day"), panel = TRUE)
plot(gg_dta, xvar = "Month", notch = TRUE)

## -------- motor trend cars data
rfsrc_mtcars <- randomForestSRC::rfsrc(mpg ~ ., data = mtcars, ntree = 50)

gg_dta <- gg_variable(rfsrc_mtcars)

# mtcars$cyl is an ordinal variable
gg_dta$cyl <- factor(gg_dta$cyl)
gg_dta$am <- factor(gg_dta$am)
gg_dta$vs <- factor(gg_dta$vs)
gg_dta$gear <- factor(gg_dta$gear)
gg_dta$carb <- factor(gg_dta$carb)

plot(gg_dta, xvar = "cyl")
plot(gg_dta, xvar = "disp")
plot(gg_dta, xvar = "hp")
plot(gg_dta, xvar = "wt")
plot(gg_dta, xvar = c("disp", "hp", "drat", "wt", "qsec"), panel = TRUE)
plot(gg_dta,
  xvar = c("cyl", "vs", "am", "gear", "carb"), panel = TRUE,
  notch = TRUE
)

## -------- Boston data
if (requireNamespace("MASS", quietly = TRUE)) {
  data(Boston, package = "MASS")
  rf_boston <- randomForest::randomForest(medv ~ ., data = Boston)
  gg_dta <- gg_variable(rf_boston)
  plot(gg_dta)
  plot(gg_dta, panel = TRUE)
}

## ------------------------------------------------------------
## survival examples
## ------------------------------------------------------------

## -------- veteran data
data(veteran, package = "randomForestSRC")
rfsrc_veteran <- randomForestSRC::rfsrc(Surv(time, status) ~ ., veteran,
  nsplit = 10,
  ntree = 50
)

# get the 90-day survival time.
gg_dta <- gg_variable(rfsrc_veteran, time = 90)

# Generate variable dependence plots for age and diagtime
plot(gg_dta, xvar = "age")
plot(gg_dta, xvar = "diagtime")

# Generate coplots
plot(gg_dta, xvar = c("age", "diagtime"), panel = TRUE, se = FALSE)

# Compare survival at 30, 90, and 365 days simultaneously
gg_dta <- gg_variable(rfsrc_veteran, time = c(30, 90, 365))
plot(gg_dta, xvar = "age")

Description

Pulls the per-tree importance scores out of a fitted varpro object and summarizes them into a data structure the plot method can draw as a boxplot. The box hinges are the 15th and 85th percentiles and the whiskers run to the 5th and 95th, not the usual Tukey 1.5 IQR whiskers. For a classification forest you can also keep the class-conditional importances.

Usage

gg_varpro(
  object,
  local.std = TRUE,
  cutoff = 0.79,
  faithful = FALSE,
  conditional = FALSE,
  nvar = NULL,
  trim = 0.1,
  max.tree = object$max.tree,
  max.rules.tree = object$max.rules.tree,
  ...
)

Arguments

Value

A named list of class "gg_varpro" with elements:

$imp

Summary data frame: variable (factor whose levels run least- to most-important by median per-tree z, so the most-important variable sits at the top of the plot after coord_flip), z (aggregate z-score, NA when it is undefined), selected (logical, z > cutoff), coverage (share of the scored trees in which a release rule for the variable returned an importance value) and tested (coverage >= trim).

$imp.tree

NULL when faithful = FALSE; otherwise an ntree x p matrix of per-tree importance values.

$stats

Per-variable summary: variable, median, q05, q15, q85, q95 of the per-tree values (on z's scale when local.std = TRUE, raw when FALSE), mean (their plain mean), and coverage and tested as in $imp.

$conditional

NULL when conditional = FALSE; otherwise a data frame with columns variable, class, z (one row per variable x class combination). class is a factor of the outcome's own level names, also for a binary outcome, which varPro recodes to 0/1 internally.

A "provenance" attribute carries family, local.std, cutoff, faithful, conditional, xvar.names, n, trim, max.tree and max.rules.tree.

What varpro is doing

Permutation importance asks "what happens to OOB accuracy when I scramble this variable?" That works, but it leans on artificial data (the permuted column) and the answer can be unstable when variables are correlated. The varpro framework (Lu and Ishwaran, 2024) replaces permutation with release rules. The forest is grown with guided splitting; from a subset of trees varpro samples a collection of decision-rule branches; for each variable it then compares the response inside the rule's region to the response after the rule's constraint on that variable is "released". The size of that change, aggregated over many rules and trees, is the variable's importance. No synthetic covariates, no permutation: the contrast is between two real subsets of the data.

Because varpro builds importance from rules sampled over trees, every tree contributes its own importance value for each variable. Those are the per-tree scores we summarize here, and the aggregate z comes from the same rules on the same trees, so the boxes and z always agree. varPro::importance() takes a winsorised mean over trees (trimmed by trim), counting a tree with no usable release rule for the variable as 0. With local.std = TRUE (the default, as in varPro) the release rules are re-sampled and scored by their local importance, the standard deviation is fixed at 1, and z is that winsorised mean, so the per-tree values already sit on z's scale. With local.std = FALSE the fit-time rules are used and z is the winsorised mean divided by the winsorised standard deviation. cutoff = 0.79 is varpro's default selection threshold, but the scale of z depends on how many variables share each tree, so the same cutoff does not mean the same thing across settings.

A variable scored in fewer than about trim of the trees (a tree scores it when one of its release rules for the variable returns an importance value) is untested: its z is 0 (or undefined with local.std = FALSE) because every non-zero tree is winsorised away, not because it was tested and found unimportant. coverage and tested in the output say which variables those are, and the plot draws them as open boxes. varpro_coverage counts coverage the same way and adds a per-variable roll-up and a noise-based cutoff.

For a classification forest, varpro also returns a class-conditional z table: the same importance computed restricting attention to rules relevant to each class. conditional = TRUE keeps that table so the plot method can show which variables matter for which class rather than only in aggregate.

What's in the output

$imp is the one-row-per-variable summary: aggregate z from varPro's importance aggregation, a selected flag for z > cutoff, and the variable's coverage and tested flag. $stats holds the box quantiles (5/15/50/85/95 percentiles, plus the raw mean) computed from the per-tree matrix; these are what the boxplot draws. $imp.tree is the per-tree matrix itself, kept only when faithful = TRUE so the plot method can scatter individual tree values over the box. $conditional is the tidy class x variable z table, present only when conditional = TRUE and the family is classification.

What you use this for

• rank candidate variables by importance and pick a working set above varpro's z cutoff;

• see, via the boxplot's spread and the per-tree points (faithful = TRUE), how stable each variable's importance is across trees: a high median with a wide box is a different story from a high median with a tight box;

• for a classification forest, ask which variables drive which class (conditional = TRUE) rather than just which variables drive the model overall.

The z-score is a standardized ranking statistic, not a p-value or a probability. Two variables with the same z are "similarly important by this method", not "equally likely to be true signal". For a data-driven cutoff rather than the 0.79 default, see varPro::cv.varpro.

Reading the plot

Rows are ordered by the median tree, but the highlight follows z, varPro's winsorised mean across the trees, against cutoff. The two can disagree: a variable a few trees find very important and many find unimportant has a long box and a low median.

A tree with no usable release rule for a variable contributes a zero, which pulls the lower whisker and box toward 0 as coverage falls; at high coverage the box need not reach 0. varPro winsorises the per-tree values (trim = 0.1: values beyond the 10th and 90th percentiles are pulled in to them), so a variable scored in fewer than about trim of the trees comes back with z exactly 0. That variable is untested, not unimportant; it is drawn as an open box, and $imp$coverage and $imp$tested say so.

By default varPro scores a random sample of 150 trees, so on a large forest the plot can change from run to run. Pass max.tree = object$rf$ntree to score every tree. On small or sparse data the 0.79 cutoff decides a great deal; see varpro_coverage and vignette("varpro-small-data", package = "ggRandomForests").

References

Lu, M. and Ishwaran, H. (2024). Model-independent variable selection via the rule-based variable priority framework. arXiv preprint arXiv:2409.09003.

See Also

plot.gg_varpro, varpro_coverage, gg_vimp

Examples

set.seed(42)
vp <- varPro::varpro(mpg ~ ., data = mtcars, ntree = 50)
gg <- gg_varpro(vp)
print(gg)
plot(gg)

Description

gg_vimp extracts the variable importance (VIMP) information from a rfsrc or randomForest object and reshapes it into a tidy data set.

Usage

gg_vimp(object, nvar, ...)

Arguments

Details

gg_vimp() reports whatever importance the forest stored; it computes nothing itself. Usually that is permutation (Breiman-Cutler) variable importance: the forest permutes a variable's observed values across the out-of-bag (OOB) cases, runs those perturbed cases down the already-grown trees, and measures how much the OOB prediction error climbs. That perturbation is synthetic (the variable's link to the response is broken on purpose) so a large increase means the variable was carrying genuine signal; near-zero or negative values mean it added noise or nothing at all.

A randomForest fit needs importance = TRUE to give you this. randomForest::randomForest() defaults to importance = FALSE, and that fit stores only a node-impurity measure (IncNodePurity, RSS, for regression; MeanDecreaseGini for classification), which is not a permutation quantity and is not comparable to one. It is the only importance the forest kept, so it is what gg_vimp() reports, in the vimp column, same as any other. Nothing marks the difference in the plot. So gg_vimp(randomForest(y ~ ., data)) ranks by node purity; pass importance = TRUE and you get permutation VIMP (%IncMSE for regression), and colnames(object$importance) tells you which one you have. randomForestSRC::rfsrc() has no such trap: its importance argument yields permutation VIMP.

When a randomForest fit carries both measures, gg_vimp() reports the permutation one and leaves node purity out of the ranking; the two run on different scales and mean different things, so putting them in one ordering would be meaningless. Read randomForest::importance(object) if you want both. A classification fit names that pair MeanDecreaseAccuracy and MeanDecreaseGini, and stores a permutation column per class besides. Those per-class columns are all permutation measures on one scale, so gg_vimp() keeps them together, names each in the set column, and drops only the Gini one.

gg_varpro() takes the opposite route, comparing local estimators on real observed data through varPro's release rules, with no permutation and no synthetic features. The two approaches answer "which variables matter?" by opposite mechanisms, so a variable can rank differently under each, and that disagreement is itself informative: it often signals interaction structure or non-monotone effects that one mechanism surfaces and the other obscures.

For survival forests, VIMP is measured against the ensemble cumulative hazard function (CHF); the error metric is one minus the concordance index (C-statistic). Variables with non-positive VIMP are flagged in the positive column and colored differently by plot.gg_vimp.

Value

gg_vimp object. A data.frame of VIMP measures in rank order, with columns vars, set (the importance measure), vimp and positive. A classification forest with class-wise permutation importance has one set per class plus the overall measure; a randomForest classification fit grown with importance = FALSE stores only node impurity, so it has a single set. When randomForest objects lack stored importance values a warning is issued and NA placeholders (columns vimp, vars and positive) are returned so plots remain reproducible.

References

Ishwaran H. (2007). Variable importance in binary regression trees and forests, Electronic J. Statist., 1:519-537.

See Also

plot.gg_vimp rfsrc

vimp gg_varpro

Examples

## ------------------------------------------------------------
## classification example
## ------------------------------------------------------------
## -------- iris data
rfsrc_iris <- randomForestSRC::rfsrc(Species ~ .,
  data = iris,
  ntree = 50,
  importance = TRUE
)
gg_dta <- gg_vimp(rfsrc_iris)
plot(gg_dta)

## ------------------------------------------------------------
## regression example
## ------------------------------------------------------------

## -------- air quality data
rfsrc_airq <- randomForestSRC::rfsrc(Ozone ~ ., airquality,
  ntree = 50,
  importance = TRUE
)
gg_dta <- gg_vimp(rfsrc_airq)
plot(gg_dta)


## -------- Boston data
if (requireNamespace("MASS", quietly = TRUE)) {
  data(Boston, package = "MASS")
  rfsrc_boston <- randomForestSRC::rfsrc(medv ~ ., Boston,
    ntree = 50,
    importance = TRUE
  )
  gg_dta <- gg_vimp(rfsrc_boston)
  plot(gg_dta)

  ## importance = TRUE for permutation VIMP; without it randomForest stores
  ## only IncNodePurity, which is what you would be ranking (see Details).
  rf_boston <- randomForest::randomForest(medv ~ ., Boston,
    importance = TRUE
  )
  gg_dta <- gg_vimp(rf_boston)
  plot(gg_dta)
}


## -------- mtcars data
rfsrc_mtcars <- randomForestSRC::rfsrc(mpg ~ .,
  data = mtcars,
  ntree = 50,
  importance = TRUE
)
gg_dta <- gg_vimp(rfsrc_mtcars)
plot(gg_dta)

## ------------------------------------------------------------
## survival example
## ------------------------------------------------------------

## -------- veteran data
data(veteran, package = "randomForestSRC")
rfsrc_veteran <- randomForestSRC::rfsrc(Surv(time, status) ~ .,
  data = veteran,
  ntree = 50,
  importance = TRUE
)

gg_dta <- gg_vimp(rfsrc_veteran)
plot(gg_dta)

## -------- pbc data
data(pbc, package = "randomForestSRC", envir = environment())
# Recode: binary columns to logical, low-cardinality columns to factors.
for (ind in seq_len(dim(pbc)[2])) {
  if (!is.factor(pbc[, ind])) {
    if (length(unique(pbc[which(!is.na(pbc[, ind])), ind])) <= 2) {
      if (sum(range(pbc[, ind], na.rm = TRUE) == c(0, 1)) == 2) {
        pbc[, ind] <- as.logical(pbc[, ind])
      }
    }
  } else {
    if (length(unique(pbc[which(!is.na(pbc[, ind])), ind])) <= 2) {
      if (sum(sort(unique(pbc[, ind])) == c(0, 1)) == 2) {
        pbc[, ind] <- as.logical(pbc[, ind])
      }
      if (sum(sort(unique(pbc[, ind])) == c(FALSE, TRUE)) == 2) {
        pbc[, ind] <- as.logical(pbc[, ind])
      }
    }
  }
  if (!is.logical(pbc[, ind]) &&
    length(unique(pbc[which(!is.na(pbc[, ind])), ind])) <= 5) {
    pbc[, ind] <- factor(pbc[, ind])
  }
}
# Convert age from days to years
pbc$age <- pbc$age / 364.24
pbc$years <- pbc$days / 364.24
pbc <- pbc[, -which(colnames(pbc) == "days")]
pbc$treatment <- as.numeric(pbc$treatment)
pbc$treatment[which(pbc$treatment == 1)] <- "DPCA"
pbc$treatment[which(pbc$treatment == 2)] <- "placebo"
pbc$treatment <- factor(pbc$treatment)
# Split into a training set (assigned treatment) and a test set (not assigned)
dta_train <- pbc[-which(is.na(pbc$treatment)), ]
pbc_test <- pbc[which(is.na(pbc$treatment)), ]

# ========
# build the forest:
rfsrc_pbc <- randomForestSRC::rfsrc(
  Surv(years, status) ~ .,
  dta_train,
  ntree = 50,
  nsplit = 10,
  na.action = "na.impute",
  forest = TRUE,
  importance = TRUE,
  save.memory = TRUE
)

gg_dta <- gg_vimp(rfsrc_pbc)
plot(gg_dta)

# Restrict to only the top 10.
gg_dta <- gg_vimp(rfsrc_pbc, nvar = 10)
plot(gg_dta)

Description

nonparametric Kaplan-Meier estimates

Usage

kaplan(interval, censor, data, by = NULL, ...)

Arguments

Value

gg_survival object

See Also

gg_survival nelson plot.gg_survival

Examples

# These get run through the gg_survival examples.
data(pbc, package = "randomForestSRC")
pbc$time <- pbc$days / 364.25

# This is the same as gg_survival
gg_dta <- kaplan(
  interval = "time", censor = "status",
  data = pbc
)

plot(gg_dta, error = "none")
plot(gg_dta)

# Stratified on treatment variable.
gg_dta <- gg_survival(
  interval = "time", censor = "status",
  data = pbc, by = "treatment"
)

plot(gg_dta, error = "none")
plot(gg_dta)

Description

cum_haz is the Nelson-Aalen estimate of the cumulative hazard: at each event time the number of events is divided by the number at risk, and the ratios are summed. The surv column, its standard error and its confidence limits are the Kaplan-Meier estimates, as kaplan returns them, and hazard, density, life and proplife are derived from that surv. The two functions therefore differ only in cum_haz, where kaplan reports $-\log S(t)$. The two agree closely while the risk set is large and diverge in the tail; when the last observation is an event, $-\log S(t)$ is infinite there and the Nelson-Aalen sum is not.

Usage

nelson(interval, censor, data, by = NULL, weight = NULL, ...)

Arguments