Package 'randomForestSRC'

Fast Unified Random Forests for Survival, Regression, and Classification (RF-SRC)

Description

Fast OpenMP-parallel implementation of Breiman's random forests (Breiman, 2001) for regression, classification, survival analysis (Ishwaran, 2008), competing risks (Ishwaran, 2012), multivariate outcomes (Segal and Xiao, 2011), unsupervised learning (Mantero and Ishwaran, 2020), quantile regression (Meinshausen, 2006; Zhang et al., 2019; Greenwald and Khanna, 2001), and imbalanced q-classification (O'Brien and Ishwaran, 2019).

Supports deterministic and randomized splitting rules (Geurts et al., 2006; Ishwaran, 2015) across all families. Variable importance (VIMP), holdout VIMP, and confidence regions (Ishwaran and Lu, 2019) can be computed for single and grouped variables. Includes minimal depth variable selection (Ishwaran et al., 2010, 2011) and a fast interface for missing data imputation using multiple forest-based methods (Tang and Ishwaran, 2017).

Tree structures can be visualized in Safari or Chrome for any family; see get.tree.

Package Overview

This package contains many useful functions. Users are encouraged to read the help files in full for detailed guidance. Below is a brief overview of key functions to help navigate the package.

1. rfsrc

   The main entry point to the package. Builds a random forest using user-supplied training data. The returned object is of class (rfsrc, grow).

2. rfsrc.fast

   A computationally efficient version of rfsrc using subsampling.

3. quantreg.rfsrc, quantreg

   Univariate and multivariate quantile regression forests for training and testing. Includes methods such as the Greenwald-Khanna (2001) algorithm, ideal for large data due to its memory efficiency.

4. predict.rfsrc, predict

   Predicts outcomes by dropping test data down the trained forest. Returns an object of class (rfsrc, predict).

5. sidClustering.rfsrc, sidClustering

   Unsupervised clustering using SID (Staggered Interaction Data). Also includes Breiman's artificial two-class method (Breiman, 2003).

6. vimp, subsample, holdout.vimp

   Functions for variable selection and importance assessment:

   1. vimp: Computes variable importance (VIMP) by perturbing each variable (e.g., via permutation). Can also be computed directly in rfsrc and predict.rfsrc.

   2. subsample: Computes confidence intervals for VIMP using subsampling.

   3. holdout.vimp: Measures the effect of removing a variable from the model.

   4. VarPro (VarPro package): For advanced model-independent variable selection using rule-based variable priority. Supports regression, classification, survival, and unsupervised data. See https://www.varprotools.org.

7. imbalanced.rfsrc, imbalanced

   Implements q-classification and G-mean-based VIMP for class-imbalanced data.

8. impute.rfsrc, impute

   A fast interface for missing data imputation. While rfsrc and predict.rfsrc can handle missing data internally, this provides a dedicated, efficient solution for imputation tasks.

9. partial.rfsrc, partial

   Computes partial dependence functions to assess the marginal effect of one or more variables on the forest ensemble.

Home page, Vignettes, Discussions, Bug Reporting, Source Code, Beta Builds

1. The package home page, with vignettes, manuals, GitHub links, and additional documentation, is available at: https://www.randomforestsrc.org/index.html

2. Questions, comments, and general usage discussions (non-bug-related) can be posted at: https://github.com/kogalur/randomForestSRC/discussions/

3. Bug reports should be submitted at: https://github.com/kogalur/randomForestSRC/issues/

   Please use this only for bugs, and include the following with your report:

   • Output from sessionInfo().

   • A minimal reproducible example including:

     • A minimal dataset required to reproduce the error.

     • The smallest runnable code needed to reproduce the issue.

     • Version details of R and all relevant packages.

     • A random seed (via set.seed()) if randomness is involved.

4. The latest stable release of the package is available on CRAN: https://cran.r-project.org/package=randomForestSRC/

5. Development builds (unstable) with bug fixes and new features are hosted on GitHub: https://github.com/kogalur/randomForestSRC/

OpenMP Parallel Processing – Installation

This package supports OpenMP shared-memory parallel programming on systems where the architecture and operating system permit it. OpenMP is enabled by default.

Detailed instructions for configuring OpenMP parallel processing can be found at: https://www.randomforestsrc.org/articles/installation.html

Note that running the package with OpenMP (or Open MPI) may increase memory (RAM) usage. Users are advised to understand their system's hardware limits and to monitor resource consumption to avoid overtaxing CPU and memory capacity.

Reproducibility

Model reproducibility is determined by three components: the random seed, the forest topology (i.e., the structure of trees), and terminal node membership for the training data. These elements together allow the model and its terminal node statistics to be faithfully restored.

Other outputs, such as variable importance (VIMP) and performance metrics, rely on additional internal randomization and are not considered part of the model definition. As a result, such statistics are subject to Monte Carlo variability and may differ across runs, even with the same seed.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

Geurts, P., Ernst, D. and Wehenkel, L., (2006). Extremely randomized trees. Machine learning, 63(1):3-42.

Greenwald M. and Khanna S. (2001). Space-efficient online computation of quantile summaries. Proceedings of ACM SIGMOD, 30(2):58-66.

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

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

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

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

Ishwaran H., Kogalur U.B., Chen X. and Minn A.J. (2011). Random survival forests for high-dimensional data. Stat. Anal. Data Mining, 4:115-132

Ishwaran H., Gerds T.A., Kogalur U.B., Moore R.D., Gange S.J. and Lau B.M. (2014). Random survival forests for competing risks. Biostatistics, 15(4):757-773.

Ishwaran H. and Malley J.D. (2014). Synthetic learning machines. BioData Mining, 7:28.

Ishwaran H. (2015). The effect of splitting on random forests. Machine Learning, 99:75-118.

Ishwaran H. and Lu M. (2019). Standard errors and confidence intervals for variable importance in random forest regression, classification, and survival. Statistics in Medicine, 38, 558-582.

Lu M., Sadiq S., Feaster D.J. and Ishwaran H. (2018). Estimating individual treatment effect in observational data using random forest methods. J. Comp. Graph. Statist, 27(1), 209-219

Mantero A. and Ishwaran H. (2021). Unsupervised random forests. Statistical Analysis and Data Mining, 14(2):144-167.

Meinshausen N. (2006) Quantile regression forests, Journal of Machine Learning Research, 7:983-999.

O'Brien R. and Ishwaran H. (2019). A random forests quantile classifier for class imbalanced data. Pattern Recognition, 90, 232-249

Segal M.R. and Xiao Y. Multivariate random forests. (2011). Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery. 1(1):80-87.

Tang F. and Ishwaran H. (2017). Random forest missing data algorithms. Statistical Analysis and Data Mining, 10:363-377.

Zhang H., Zimmerman J., Nettleton D. and Nordman D.J. (2019). Random forest prediction intervals. The American Statistician. 4:1-5.

See Also

get.tree.rfsrc,

holdout.vimp.rfsrc,

imbalanced.rfsrc, impute.rfsrc,

max.subtree.rfsrc,

partial.rfsrc, plot.competing.risk.rfsrc, plot.rfsrc, plot.survival.rfsrc, plot.variable.rfsrc, predict.rfsrc, print.rfsrc,

quantreg.rfsrc,

rfsrc, rfsrc.cart, rfsrc.fast,

sidClustering.rfsrc,

subsample.rfsrc,

tune.rfsrc,

vimp.rfsrc

---

Wisconsin Prognostic Breast Cancer Data

Description

Recurrence of breast cancer from 198 breast cancer patients, all of which exhibited no evidence of distant metastases at the time of diagnosis. The first 30 features of the data describe characteristics of the cell nuclei present in the digitized image of a fine needle aspirate (FNA) of the breast mass.

Source

The data were obtained from the UCI machine learning repository, see http://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Prognostic).

Examples

data(breast, package = "randomForestSRC")
breast <- na.omit(breast)
o <- rfsrc(status ~ ., data = breast, nsplit = 10)
print(o)

---

Classification Performance Metrics

Description

Evaluate predicted classes and probabilities using confusion matrices, class-specific misclassification rates, ROC area under the curve (AUC), Brier score, log loss, and binary precision-recall summaries. get.bayes.rule() assigns classes from predicted probabilities.

Usage

get.confusion(y, class.or.prob)

get.misclass.error(y, yhat)

get.auc(y, prob)

get.brier.error(y, prob, normalized = TRUE, vector = FALSE)

get.logloss(y, prob, robust = TRUE)

get.pr.auc(truth, yhat)

get.pr.curve(truth, yhat)

get.bayes.rule(prob, class.relfrq = NULL)

Arguments

y	Observed class labels, usually a factor with levels in probability-column order. get.auc() also accepts nonfactor labels, with columns in sort(unique(y)) order. get.misclass.error() also accepts vectors of class labels.
prob	Numeric matrix with one row per observation and one column per class, in response-level order. Name columns with the class labels. get.auc() uses column positions; get.logloss() and get.bayes.rule() require column names. get.confusion() and get.brier.error() supply missing names from levels(y).
class.or.prob	A factor of predicted labels with the same levels as y, or a probability matrix as described for prob.
yhat	For get.misclass.error(), predicted class labels aligned with y. For precision-recall summaries, a numeric score vector or a matrix or data frame with one or two columns. A vector or single column gives positive-class scores; larger values favor that class. Scores need not lie in $[0,1]$. See Details for two-column inputs.
truth	Binary responses coded 0/1, or a two-level factor. For 0/1 labels, 1 is positive. For other factors, the less frequent class is positive, with ties resolved by the first factor level. The positive class is selected before excluding missing scores.
normalized	If TRUE, scale the Brier score so that equal class probabilities give score one. If FALSE, average the squared probability errors over classes. See Details.
vector	If TRUE, return scaled Brier contributions for each class instead of their sum.
robust	For get.logloss(), exclude infinite losses before averaging. Probabilities are unchanged.
class.relfrq	Optional two-class relative frequencies, in probability-column order. NULL selects the largest-probability rule; supplying frequencies selects the minority-prevalence rule described below.

Details

Supplying predictions

Supply responses and predictions in matching row order. For OOB evaluation of a grow object o, use o$yvar with o$predicted.oob for probabilities or o$class.oob for predicted classes. For new data, use the prediction object's yvar with predicted or class.

To compare metrics on the same observations, retain observed labels and finite probabilities. get.logloss() and get.misclass.error() require nonmissing observed labels.

For a combined binary summary with sensitivity, specificity, F1, G-mean, and random-reference comparisons, use get.imbalanced.performance.

Class assignments and misclassification

With class.relfrq = NULL, get.bayes.rule() selects the class with the largest probability, breaking ties at random. Rows with all probabilities missing receive NA. With two class frequencies supplied, it assigns the minority class when its probability is at least its relative frequency, and the majority class otherwise. If frequencies tie, the class in the first probability column is treated as the minority. Supply finite probabilities for this rule.

get.confusion() tabulates observed classes in rows and predicted classes in columns. The final class.error column gives each observed class's misclassification rate, rounded to four decimal places. Missing response/prediction pairs are excluded. Probability matrices use the largest-probability rule. For RFQ or another thresholded rule, supply the predicted factor instead, such as o$class.oob.

get.misclass.error() returns one error rate per observed class, in sort(unique(y)) order. A missing predicted label makes the corresponding class rate unavailable.

ROC area under the curve

get.auc() computes the multiclass AUC of Hand and Till (2001). For each class pair, it averages two rank-based AUCs, one using each class's probability as the score, then averages across available pairs. For complementary binary probabilities, the two directions agree. Larger values indicate better discrimination.

Missing labels and nonfinite scores are excluded from each rank calculation. Each direction requires at least two finite scores per class. Unavailable pairs are omitted; the result is NA when none remain.

Brier score

The Brier score measures squared probability error (Brier, 1950). For $J$ classes, let

$b_j=\mathrm{mean}_i\left[(I(y_i=j)-p_{ij})^2\right].$

get.brier.error() returns

$\frac{J}{J-1}\sum_{j=1}^{J} b_j$

for normalized = TRUE, or

$\frac{1}{J}\sum_{j=1}^{J} b_j$

for normalized = FALSE. Smaller values are better. Equal probabilities $p_{ij}=1/J$ give normalized score one. For complementary binary probabilities, the unnormalized score is the mean squared error of either class probability; the normalized score is four times that value.

With vector = TRUE, each element is $b_j$ multiplied by the selected scaling constant. Missing losses are omitted within each class. The scalar sums available contributions, returning NA when none remain.

Log loss

With every response level represented, get.logloss() averages the negative log probability of the observed class, $-\log(p_{i,y_i})$, using the natural logarithm (see Gneiting and Raftery, 2007). Smaller values are better. Missing losses are omitted. With robust = TRUE, infinite losses from zero probabilities are also omitted; robust = FALSE retains them. Probabilities are not clipped. An unused factor level contributes one zero term to the average.

Precision-recall summaries

For two-column scores, named columns match the original response levels; unnamed columns follow factor-level order, or 0, 1 order for nonfactor responses. The positive-class column is selected. Rows with a missing response or nonfinite selected score are excluded; both classes must remain.

get.pr.auc() returns the precision-recall area, calculated using analytic precision-recall interpolation, and a random reference area equal to the positive-class proportion among scored observations. get.pr.curve() returns recall, precision, and thresholds in decreasing recall order. Larger precision-recall areas are better. See Davis and Goadrich (2006) for background on precision-recall curves.

Value

get.confusion

A numeric matrix of class counts with an additional class.error column.

get.misclass.error

An unnamed numeric vector of class-specific error rates, in sort(unique(y)) order.

get.auc

A scalar AUC, or NA when unavailable.

get.brier.error

A scalar Brier score by default, or an unnamed length-$J$ vector of scaled class contributions in probability-column order when vector = TRUE.

get.logloss

A scalar mean log loss. It can be infinite with robust = FALSE, or NaN when no losses remain for averaging.

get.pr.auc

An unnamed numeric vector of length two: the model area followed by the random reference area. Both entries are NA when the calculation is unavailable.

get.pr.curve

A numeric matrix with columns recall, precision, and threshold, or NULL when the calculation is unavailable.

get.bayes.rule

A factor of predicted classes with levels given by the probability-column names.

References

Brier, G.W. (1950). Verification of forecasts expressed in terms of probability. Monthly Weather Review, 78, 1-3. doi:10.1175/1520-0493(1950)078<0001:VOFEIT>2.0.CO;2

Davis, J. and Goadrich, M. (2006). The relationship between precision-recall and ROC curves. Proceedings of the 23rd International Conference on Machine Learning, 233-240. doi:10.1145/1143844.1143874

Gneiting, T. and Raftery, A.E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102, 359-378. doi:10.1198/016214506000001437

Hand, D.J. and Till, R.J. (2001). A simple generalisation of the area under the ROC curve for multiple class classification problems. Machine Learning, 45, 171-186. doi:10.1023/A:1010920819831

See Also

rfsrc, predict.rfsrc, imbalanced.rfsrc, get.imbalanced.performance, get.brier.survival, get.auct.survival

Examples

## ------------------------------------------------------------
## A basic calculation from observed labels and probabilities
## ------------------------------------------------------------
y <- factor(c("no", "no", "no", "no", "yes", "yes"),
            levels = c("no", "yes"))
p <- c(.10, .20, .65, .35, .40, .85)
prob <- cbind(no = 1 - p, yes = p)
print(get.confusion(y, prob))
print(get.auc(y, prob))
print(get.brier.error(y, prob))


## ------------------------------------------------------------
## Class-specific errors and probability losses
## ------------------------------------------------------------
yhat <- get.bayes.rule(prob)
print(setNames(get.misclass.error(y, yhat), levels(y)))
print(get.brier.error(y, prob, normalized = FALSE))
print(setNames(get.brier.error(y, prob, vector = TRUE), colnames(prob)))
print(get.logloss(y, prob))

## Use supplied class frequencies for the binary RFQ decision rule.
class.frq <- as.numeric(prop.table(table(y)))
yhat.rfq <- get.bayes.rule(prob, class.relfrq = class.frq)
print(get.confusion(y, yhat.rfq))

## ------------------------------------------------------------
## Precision-recall from a score vector
## ------------------------------------------------------------
truth <- as.integer(y == "yes")
pr.auc <- get.pr.auc(truth, p)
print(setNames(pr.auc, c("model", "random")))
pr <- get.pr.curve(truth, p)
plot(pr[, "recall"], pr[, "precision"], type = "l",
     xlim = c(0, 1), ylim = c(0, 1),
     xlab = "Recall", ylab = "Precision")
abline(h = pr.auc[2], lty = 2)

## A single score column is also accepted.
print(get.pr.auc(truth, matrix(p, ncol = 1)))

## ------------------------------------------------------------
## Multiclass forest: OOB and test-data performance
## ------------------------------------------------------------
set.seed(17)
train <- c(1:35, 51:85, 101:135)
o <- rfsrc(Species ~ ., data = iris[train, ], ntree = 100)

## Select one common set of observed responses and finite OOB scores.
p.oob <- o$predicted.oob
keep <- !is.na(o$yvar) & rowSums(!is.finite(p.oob)) == 0
print(get.confusion(o$yvar[keep], o$class.oob[keep]))
print(get.auc(o$yvar[keep], p.oob[keep, , drop = FALSE]))

## Use the current responses and probabilities for test-data scoring.
p.test <- predict(o, newdata = iris[-train, ])
print(c(
  auc = get.auc(p.test$yvar, p.test$predicted),
  brier = get.brier.error(p.test$yvar, p.test$predicted),
  logloss = get.logloss(p.test$yvar, p.test$predicted)
))

---

Fast Saving and Loading of Random Forests

Description

Save and load a minimal forest object for prediction. These functions are designed to use substantially less disk space and load large forests faster than saving the full fitted object with saveRDS and loading it with readRDS. Functions are available for a single forest or a list of forests.

Usage

fast.save(o, path = NULL, testing = TRUE, units = "Mb", compact = FALSE,
          compress = 50)

fast.load(directory, path = NULL, testing = FALSE, units = "Mb")

fast.save.list(o, path = NULL, testing = FALSE, units = "Mb", compact = FALSE,
               compress = 50)

fast.load.list(directory, path = NULL, testing = FALSE, units = "Mb")

Arguments

o	For fast.save(), a fitted object with class beginning c("rfsrc", "grow") and a forest component. Supply the fitted object, not o$forest. For fast.save.list(), a nonempty list of these objects.
path	For saving, the destination directory, including its name. The default is file.path(getwd(), "forest"). An existing destination and all its contents are deleted. For loading, the parent directory containing directory; NULL uses the current working directory.
directory	Name of the saved directory relative to path. For fast.load.list(), this directory contains the individual forest subdirectories.
testing	Print component memory sizes using lsos(). Defaults to TRUE for fast.save() and FALSE for the other functions.
units	Units for the memory-size display, passed to format(object.size(...), units = units). The default is "Mb".
compact	If TRUE, omit terminal-node information from the saved forest and recompute it during prediction. Requires a non-anonymous forest containing the training predictors and responses. If FALSE (default), save the forest with its available terminal-node information. See Compact saving below.
compress	A single finite number between 0 and 100 specifying the FST compression level. The default is 50 for fast saving. Use 100 for maximum compression when smaller files are more important than saving speed. This setting is independent of compact.

Details

Saving and loading a forest

fast.save() addresses the large files and slow loading that can result from saving a full fitted forest with saveRDS(). It saves only o$forest, the minimal object used for prediction, including its tree information and any training data it contains. Predictions and performance summaries outside this component are excluded. fast.load() returns the forest for use with predict.rfsrc.

Large tabular components are saved efficiently with fst::write_fst(), and the remaining forest information with saveRDS(..., compress = FALSE). Saving requires the fst and data.table packages; loading requires fst. All files have a ‘.rda’ extension despite the different formats. Keep the directory's files together and load the forest with fast.load().

Compression

The compress setting is passed to every fst::write_fst() call, including each forest saved by fast.save.list(). Higher compression can reduce file sizes at the cost of slower saving; the reduction depends on the saved components. The remaining forest information in ‘forest.rda’ continues to use saveRDS(..., compress = FALSE).

Compression changes the encoding of the saved components, not the fitted forest or the choice of components. It can be combined with compact = TRUE. Use the same fast.load() or fast.load.list() call for any compression level; no compression argument is needed when loading. Invalid compression values are rejected before the destination is deleted.

Compact saving

With compact = TRUE, the saved forest includes the tree structure, factor-split information, training data, seeds, and other forest information, but omits stored terminal-node information. Prediction recomputes this information from the saved forest and training data without growing new trees. The fitted object in memory is unchanged.

The saved copy has terminal.qualts = FALSE and terminal.quants = FALSE, with nativeArrayTNDS removed. Use the same fast.load() or fast.load.list() call for either saving mode. Compact saving can further reduce disk space, but prediction may take longer and peak memory use may not decrease.

Use compact = FALSE for anonymous forests. Compact saving includes the training data. Eligibility checks run before the destination is deleted; for lists, every forest is checked first.

Shared imputer banks

Standalone fast.save() always writes a self-contained forest. Imputers fitted with share.xvar = TRUE use the same component serializers but store their repeated training predictors once per bank. Load those complete banks with load.impute.learn. Direct fast.load() calls on shared learner subdirectories stop because their predictors depend on the bank's common table and reconstruction mappings.

Directories

Both save functions delete an existing destination directory and all its contents. Use a dedicated directory for each forest or list. The destination and any missing parent directories are then created.

When saving, path is the complete destination. When loading, directory is its name and path is its parent. For a destination save.path, use fast.load(basename(save.path), path = dirname(save.path)). For the default destination, use fast.load("forest").

Lists of forests

fast.save.list() saves forests in subdirectories ‘forest1’, ‘forest2’, and so on, inside path. List names are not saved. fast.load.list() loads every entry listed by list.files() in that directory, so it should contain only saved forest subdirectories.

The loaded list is unnamed and follows directory-listing order, which can differ from the original order: ‘forest10’ can precede ‘forest2’. To load one forest, call fast.load() with its subdirectory name.

Value

fast.save

Saves the forest to disk and returns the memory-usage matrix produced by gc(FALSE).

fast.load

The saved forest component with its original class, usually c("rfsrc", "forest", family).

fast.save.list

An unnamed list of the memory-usage matrices returned by the individual fast.save() calls.

fast.load.list

An unnamed list of forest objects in directory-listing order.

Note

The sampling-size function is saved as sampfrac = o$forest$sampsize(1) and loaded as function(x) x * sampfrac. Use saveRDS() to save a custom function that is not proportional to sample size.

Compact saving checks that round(sampsize(n)) equals round(n * sampsize(1)) at the training sample size n. A mismatch stops saving before the destination is deleted. This check applies only at the training sample size.

See Also

rfsrc, predict.rfsrc, saveRDS, readRDS

Examples

## ------------------------------------------------------------
## Regression: save, load, and compare predictions
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  o <- rfsrc(mpg ~ ., data = mtcars)
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  p <- predict(o)
  pp <- predict(oo)
  print(summary(p$predicted - pp$predicted))
  print(summary(p$predicted.oob - pp$predicted.oob))
  unlink(save.path, recursive = TRUE)
}


## ------------------------------------------------------------
## Compact saving: the same forest, fewer saved components
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  o <- rfsrc(mpg ~ ., data = mtcars, ntree = 100)
  set.seed(19)
  reference <- predict(o, seed = -19)
  save.path <- tempfile("rfsrc-compact-")
  fast.save(o, path = save.path, compact = TRUE, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  print(oo[c("terminal.qualts", "terminal.quants")])
  print(list.files(save.path, pattern = "^nativeArrayTDNS_"))
  set.seed(19)
  restored <- predict(oo, seed = -19)
  print(all.equal(reference$predicted, restored$predicted))
  print(all.equal(reference$predicted.oob, restored$predicted.oob))
  unlink(save.path, recursive = TRUE)
}

## ------------------------------------------------------------
## Regression: a list of forests with different node sizes
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  o1 <- rfsrc(mpg ~ ., data = mtcars, nodesize = 1)
  o2 <- rfsrc(mpg ~ ., data = mtcars, nodesize = 10)
  print(o1)
  print(o2)
  models <- list(o1, o2)
  save.path <- tempfile("rfsrc-forest-list-")
  invisible(fast.save.list(models, path = save.path))
  oo <- fast.load.list(basename(save.path), path = dirname(save.path))
  print(predict(oo[[1]]))
  print(predict(oo[[2]]))
  unlink(save.path, recursive = TRUE)
}

## ------------------------------------------------------------
## RFQ for imbalanced classification
## ------------------------------------------------------------
## Use matching prediction seeds when comparing class labels.
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(breast, package = "randomForestSRC")
  dta <- na.omit(breast)
  o <- imbalanced(status ~ ., data = dta, ntree = 100)
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, seed = -19)
  set.seed(19)
  pp <- predict(oo, seed = -19)
  print(summary(p$predicted - pp$predicted))
  print(summary(p$predicted.oob - pp$predicted.oob))
  print(all.equal(as.character(p$class), as.character(pp$class)))
  print(all.equal(as.character(p$class.oob),
                  as.character(pp$class.oob)))
  unlink(save.path, recursive = TRUE)
}

## ------------------------------------------------------------
## Binary classification with rfq = TRUE
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(breast, package = "randomForestSRC")
  dta <- na.omit(breast)
  o <- rfsrc(status ~ ., data = dta, rfq = TRUE, ntree = 100,
             perf.type = "gmean", splitrule = "auc")
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, seed = -19)
  set.seed(19)
  pp <- predict(oo, seed = -19)
  print(summary(p$predicted - pp$predicted))
  print(summary(p$predicted.oob - pp$predicted.oob))
  print(all.equal(as.character(p$class), as.character(pp$class)))
  print(all.equal(as.character(p$class.oob),
                  as.character(pp$class.oob)))
  unlink(save.path, recursive = TRUE)
}

## ------------------------------------------------------------
## Anonymous RFQ: supply the same prediction data to both forests
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(breast, package = "randomForestSRC")
  dta <- na.omit(breast)
  o <- rfsrc.anonymous(status ~ ., data = dta, rfq = TRUE,
                       ntree = 100, perf.type = "gmean", splitrule = "auc")
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, newdata = dta, seed = -19)
  set.seed(19)
  pp <- predict(oo, newdata = dta, seed = -19)
  print(summary(p$predicted - pp$predicted))
  print(all.equal(as.character(p$class), as.character(pp$class)))
  unlink(save.path, recursive = TRUE)
}

## ------------------------------------------------------------
## Survival
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(pbc, package = "randomForestSRC")
  o <- rfsrc(Surv(days, status) ~ ., data = pbc, ntree = 100)
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, seed = -19)
  set.seed(19)
  pp <- predict(oo, seed = -19)
  print(summary(p$predicted - pp$predicted))
  print(summary(p$predicted.oob - pp$predicted.oob))
  unlink(save.path, recursive = TRUE)
}

## ------------------------------------------------------------
## Survival with save.memory = TRUE
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(pbc, package = "randomForestSRC")
  o <- rfsrc(Surv(days, status) ~ ., data = pbc,
             ntree = 100, save.memory = TRUE)
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, seed = -19)
  set.seed(19)
  pp <- predict(oo, seed = -19)
  print(summary(p$predicted - pp$predicted))
  print(summary(p$predicted.oob - pp$predicted.oob))
  unlink(save.path, recursive = TRUE)
}

## ------------------------------------------------------------
## Competing risks
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(wihs, package = "randomForestSRC")
  o <- rfsrc(Surv(time, status) ~ ., data = wihs, nsplit = 3, ntree = 100)
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, seed = -19)
  set.seed(19)
  pp <- predict(oo, seed = -19)
  print(summary(p$predicted - pp$predicted))
  print(summary(p$predicted.oob - pp$predicted.oob))
  unlink(save.path, recursive = TRUE)
}

## ------------------------------------------------------------
## Multivariate regression and classification
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(nutrigenomic, package = "randomForestSRC")
  ydta <- data.frame(diet = nutrigenomic$diet,
                     genotype = nutrigenomic$genotype,
                     nutrigenomic$lipids)
  o <- rfsrc(get.mv.formula(colnames(ydta)),
             data = data.frame(ydta, nutrigenomic$genes),
             ntree = 100, importance = TRUE, nsplit = 10)
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, seed = -19)
  set.seed(19)
  pp <- predict(oo, seed = -19)
  print(summary(get.mv.predicted(p, oob = FALSE) -
                get.mv.predicted(pp, oob = FALSE)))
  print(summary(get.mv.predicted(p) - get.mv.predicted(pp)))
  for (yn in names(p$classOutput)) {
    print(yn)
    print(all.equal(as.character(p$classOutput[[yn]]$class),
                    as.character(pp$classOutput[[yn]]$class)))
    print(all.equal(as.character(p$classOutput[[yn]]$class.oob),
                    as.character(pp$classOutput[[yn]]$class.oob)))
  }
  unlink(save.path, recursive = TRUE)
}


## Not run: 

## ------------------------------------------------------------
## Classification: optional alzheimers data from varPro
## ------------------------------------------------------------
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(alzheimers, package = "varPro")
  o <- rfsrc(Diagnosis ~ ., data = alzheimers)
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, seed = -19)
  set.seed(19)
  pp <- predict(oo, seed = -19)
  print(summary(p$predicted - pp$predicted))
  print(summary(p$predicted.oob - pp$predicted.oob))
  print(all.equal(as.character(p$class), as.character(pp$class)))
  print(all.equal(as.character(p$class.oob),
                  as.character(pp$class.oob)))
  unlink(save.path, recursive = TRUE)
}


## ------------------------------------------------------------
## Optional memory-intensive anonymous survival test
## ------------------------------------------------------------
## This test repeats each PBC row 250 times and can require substantial memory.
if (requireNamespace("fst", quietly = TRUE) &&
    requireNamespace("data.table", quietly = TRUE)) {
  data(pbc, package = "randomForestSRC")
  dta <- pbc[rep(seq_len(nrow(pbc)), each = 250), ]
  o <- rfsrc.anonymous(Surv(days, status) ~ ., data = dta)
  print(o)
  save.path <- tempfile("rfsrc-forest-")
  fast.save(o, path = save.path, testing = FALSE)
  oo <- fast.load(basename(save.path), path = dirname(save.path))
  set.seed(19)
  p <- predict(o, newdata = dta, seed = -19)
  set.seed(19)
  pp <- predict(oo, newdata = dta, seed = -19)
  print(summary(p$predicted - pp$predicted))
  unlink(save.path, recursive = TRUE)
}

## End(Not run)

---

Follicular Cell Lymphoma

Description

Competing risk data set involving follicular cell lymphoma.

Format

A data frame containing:

age	age
hgb	hemoglobin (g/l)
clinstg	clinical stage: 1=stage I, 2=stage II
ch	chemotherapy
rt	radiotherapy
time	first failure time
status	censoring status: 0=censored, 1=relapse, 2=death

Source

Table 1.4b, Competing Risks: A Practical Perspective.

References

Pintilie M., (2006) Competing Risks: A Practical Perspective. West Sussex: John Wiley and Sons.

Examples

data(follic, package = "randomForestSRC")
follic.obj <- rfsrc(Surv(time, status) ~ ., follic, nsplit = 3, ntree = 100)

---

Extract a Single Tree from a Forest and plot it on your browser

Description

Extracts a single tree from a forest which can then be plotted on the users browser. Works for all families. Missing data not permitted.

Usage

## S3 method for class 'rfsrc'
get.tree(object, tree.id, target, m.target = NULL,
   time, surv.type = c("mort", "rel.freq", "surv", "years.lost", "cif", "chf"),
   class.type = c("bayes", "rfq", "prob"),
   ensemble = FALSE, oob = TRUE, show.plots = TRUE, do.trace = FALSE)

Arguments

object	An object of class (rfsrc, grow).
tree.id	Integer specifying the tree to extract.
target	For classification: integer or character indicating the class of interest (defaults to the first class). For competing risks: integer between 1 and J (number of event types) specifying the event of interest (default is the first event type).
m.target	Character string specifying the target outcome for multivariate families. If unspecified, a default is selected.
time	For survival: time point at which the predicted value is evaluated (depends on surv.type).
surv.type	For survival: specifies the type of predicted value returned. See Details.
class.type	For classification: specifies the type of predicted value. See Details.
ensemble	Logical. If TRUE, prediction is based on the ensemble of all trees. If FALSE (default), prediction is based on the specified tree.
oob	Logical. Use OOB predicted values (TRUE) or in-bag values (FALSE). Only applies when ensemble=TRUE.
show.plots	Logical. Should plots be displayed?
do.trace	Number of seconds between progress updates.

Details

Extracts a specified tree from a forest and converts it into a hierarchical structure compatible with the data.tree package. Plotting the resulting object renders an interactive tree visualization in the user's web browser.

Left-hand splits are shown. For continuous variables, the left split is displayed as an inequality (e.g., x < value); the right split is the reverse. For factor variables, the left daughter node is defined by a set of levels assigned to it; the right daughter is its complement.

Terminal nodes are highlighted with color and display both sample size and predicted value. By default, the predicted value corresponds to the prediction from the selected tree, and the sample size refers to the in-bag cases reaching the terminal node. If ensemble = TRUE, the predicted value equals the forest ensemble prediction, allowing visualization of the full forest predictor over the selected tree's partition. In this case, sample sizes refer to all observations (not just in-bag cases).

Predicted values displayed in terminal nodes are defined as follows:

1. For regression: the mean of the response.

2. For classification: depends on the class.type argument and target class:

   • If class.type = "bayes", the predicted class with the most votes, or the RFQ classifier threshold in two-class problems.

   • If class.type = "prob", the class probability for the target class.

3. For multivariate families: the predicted value for the outcome specified by m.target, using the logic above depending on whether the outcome is continuous or categorical.

4. For survival:

   • mort: estimated mortality (Ishwaran et al., 2008).

   • rel.freq: relative frequency of mortality.

   • surv: predicted survival probability at the specified time (time).

5. For competing risks:

   • years.lost: expected number of life years lost.

   • cif: cumulative incidence function.

   • chf: cause-specific cumulative hazard function.

   For cif and chf, predictions are evaluated at the time point given by time, and all metrics are specific to the event type indicated by target.

Value

Invisibly, returns an object with hierarchical structure formatted for use with the data.tree package.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

Many thanks to @dbarg1 on GitHub for the initial prototype of this function

Examples

## ------------------------------------------------------------
## survival/competing risk
## ------------------------------------------------------------

## survival - veteran data set but with factors
## note that diagtime has many levels
data(veteran, package = "randomForestSRC")
vd <- veteran
vd$celltype=factor(vd$celltype)
vd$diagtime=factor(vd$diagtime)
vd.obj <- rfsrc(Surv(time,status)~., vd, ntree = 100, nodesize = 5)
plot(get.tree(vd.obj, 3))

## competing risks
data(follic, package = "randomForestSRC")
follic.obj <- rfsrc(Surv(time, status) ~ ., follic, nsplit = 3, ntree = 100)
plot(get.tree(follic.obj, 2))

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

airq.obj <- rfsrc(Ozone ~ ., data = airquality)
plot(get.tree(airq.obj, 10))

## ------------------------------------------------------------
## two-class imbalanced data (see imbalanced function)
## ------------------------------------------------------------

data(breast, package = "randomForestSRC")
breast <- na.omit(breast)
f <- as.formula(status ~ .)
breast.obj <- imbalanced(f, breast)

## compare RFQ to Bayes Rule
plot(get.tree(breast.obj, 1, class.type = "rfq", ensemble = TRUE))
plot(get.tree(breast.obj, 1, class.type = "bayes", ensemble = TRUE))

## ------------------------------------------------------------
## classification
## ------------------------------------------------------------

iris.obj <- rfsrc(Species ~., data = iris, nodesize = 10)

## equivalent
plot(get.tree(iris.obj, 25))
plot(get.tree(iris.obj, 25, class.type = "bayes"))

## predicted probability displayed for terminal nodes
plot(get.tree(iris.obj, 25, class.type = "prob", target = "setosa"))
plot(get.tree(iris.obj, 25, class.type = "prob", target = "versicolor"))
plot(get.tree(iris.obj, 25, class.type = "prob", target = "virginica"))


## ------------------------------------------------------------
## multivariate regression
## ------------------------------------------------------------

mtcars.mreg <- rfsrc(Multivar(mpg, cyl) ~., data = mtcars)
plot(get.tree(mtcars.mreg, 10, m.target = "mpg"))
plot(get.tree(mtcars.mreg, 10, m.target = "cyl"))


## ------------------------------------------------------------
## multivariate mixed outcomes
## ------------------------------------------------------------

mtcars2 <- mtcars
mtcars2$carb <- factor(mtcars2$carb)
mtcars2$cyl <- factor(mtcars2$cyl)
mtcars.mix <- rfsrc(Multivar(carb, mpg, cyl) ~ ., data = mtcars2)
plot(get.tree(mtcars.mix, 5, m.target = "cyl"))
plot(get.tree(mtcars.mix, 5, m.target = "carb"))

## ------------------------------------------------------------
## unsupervised analysis
## ------------------------------------------------------------

mtcars.unspv <- rfsrc(data = mtcars)
plot(get.tree(mtcars.unspv, 5))

---

Hodgkin's Disease

Description

Competing risk data set involving Hodgkin's disease.

Format

A data frame containing:

age	age
sex	gender
trtgiven	treatment: RT=radition, CMT=Chemotherapy and radiation
medwidsi	mediastinum involvement: N=no, S=small, L=Large
extranod	extranodal disease: Y=extranodal disease, N=nodal disease
clinstg	clinical stage: 1=stage I, 2=stage II
time	first failure time
status	censoring status: 0=censored, 1=relapse, 2=death

Source

Table 1.6b, Competing Risks: A Practical Perspective.

References

Pintilie M., (2006) Competing Risks: A Practical Perspective. West Sussex: John Wiley and Sons.

Examples

data(hd, package = "randomForestSRC")

---

Hold out variable importance (VIMP)

Description

Hold out VIMP is calculated from the error rate of mini ensembles of trees (blocks of trees) grown with and without a variable. Applies to all families.

Usage

## S3 method for class 'rfsrc'
holdout.vimp(formula, data,
  ntree = function(p, vtry){1000 * p / vtry},
  nsplit = 10,
  ntime = 50,
  sampsize = function(x){x * .632},
  samptype = "swor",
  block.size = 10,
  vtry = 1,
  ...)

Arguments

formula	A symbolic description of the model to be fit.
data	Data frame containing the y-outcome and x-variables.
ntree	Specifies the number of trees used to grow the forest. Can be a function of data dimension and number of holdout variables, or a fixed numeric value.
nsplit	Non-negative integer specifying the number of random split points used to split a node. A value of zero corresponds to deterministic splitting, which is significantly slower.
ntime	Integer value used for survival settings to constrain ensemble calculations to a grid of ntime time points.
sampsize	Specifies the size of the subsampled data. Can be either a function or a numeric value.
samptype	Type of bootstrap used when subsampling.
vtry	Number of variables randomly selected to be held out when growing a tree. Can also be a list for targeted holdout variable importance analysis. See details for more information.
block.size	Specifies the number of trees in a block when calculating holdout variable importance.
...	Further arguments passed to rfsrc.

Details

Holdout variable importance (holdout VIMP) measures the importance of a variable by comparing prediction error between two forests (blocks of trees): one in which selected variables are held out during tree growing (the holdout forest) and one in which no variables are held out (the baseline forest).

For each variable-block combination, the bootstrap samples used to grow the trees are the same in both forests. The difference in out-of-bag (OOB) prediction error between the holdout and baseline forests gives the holdout VIMP for that variable-block pair. The final holdout VIMP for a variable is the average of these differences over all blocks in which the variable was held out.

The option vtry controls how many variables are held out per tree. The default is one, meaning a single variable is held out per tree. Larger values of vtry increase the number of times each variable is held out, reducing the required total number of trees. However, interpretation of holdout VIMP changes when vtry exceeds one, and this option should be used cautiously.

High accuracy requires a sufficiently large number of trees. As a general guideline, we recommend using ntree = 1000 * p / vtry, where p is the number of features. Accuracy also depends on block.size, which determines how many trees comprise a block. Smaller values yield better accuracy but are computationally more demanding. The most accurate setting is block.size = 1. Ensure that block.size does not exceed ntree / p, otherwise insufficient trees may be available for certain variables.

Targeted holdout VIMP analysis can be requested by specifying vtry as a list with two components: a vector of variable indices (xvar) and a logical flag joint indicating whether to compute joint VIMP. For example, to compute holdout VIMP only for variables 1, 4, and 5 individually:

vtry = list(xvar = c(1, 4, 5), joint = FALSE)

To compute the joint effect of removing these three variables together:

vtry = list(xvar = c(1, 4, 5), joint = TRUE)

Targeted analysis is useful when the user has prior knowledge of variables of interest and can significantly reduce computation. Joint VIMP quantifies the combined importance of specific groups of variables. See the Iris example below for illustration.

Value

Invisibly a list with the following components (which themselves can be lists):

importance	Holdout VIMP.
baseline	Prediction error for the baseline forest.
holdout	Prediction error for the holdout forest.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Lu M. and Ishwaran H. (2018). Expert Opinion: A prediction-based alternative to p-values in regression models. J. Thoracic and Cardiovascular Surgery, 155(3), 1130–1136.

See Also

vimp.rfsrc

Examples

## ------------------------------------------------------------
## regression analysis
## ------------------------------------------------------------

## new York air quality measurements
airq.obj <- holdout.vimp(Ozone ~ ., data = airquality, na.action = "na.impute")
print(airq.obj$importance)

## ------------------------------------------------------------
## classification analysis
## ------------------------------------------------------------

## iris data
iris.obj <- holdout.vimp(Species ~., data = iris)
print(iris.obj$importance)

## iris data using brier prediction error
iris.obj <- holdout.vimp(Species ~., data = iris, perf.type = "brier")
print(iris.obj$importance)

## ------------------------------------------------------------
## illustration of targeted holdout vimp analysis
## ------------------------------------------------------------

## iris data - only interested in variables 3 and 4
vtry <- list(xvar = c(3, 4), joint = FALSE)
print(holdout.vimp(Species ~., data = iris, vtry = vtry)$impor)

## iris data - joint importance of variables 3 and 4
vtry <- list(xvar = c(3, 4), joint = TRUE)
print(holdout.vimp(Species ~., data = iris, vtry = vtry)$impor)

## iris data - joint importance of variables 1 and 2
vtry <- list(xvar = c(1, 2), joint = TRUE)
print(holdout.vimp(Species ~., data = iris, vtry = vtry)$impor)


## ------------------------------------------------------------
## imbalanced classification (using RFQ)
## ------------------------------------------------------------

if (library("caret", logical.return = TRUE)) {

  ## experimental settings
  n <- 400
  q <- 20
  ir <- 6
  f <- as.formula(Class ~ .)
 
  ## simulate the data, create minority class data
  d <- twoClassSim(n, linearVars = 15, noiseVars = q)
  d$Class <- factor(as.numeric(d$Class) - 1)
  idx.0 <- which(d$Class == 0)
  idx.1 <- sample(which(d$Class == 1), sum(d$Class == 1) / ir , replace = FALSE)
  d <- d[c(idx.0,idx.1),, drop = FALSE]

  ## VIMP for RFQ with and without blocking
  vmp1 <- imbalanced(f, d, importance = TRUE, block.size = 1)$importance[, 1]
  vmp10 <- imbalanced(f, d, importance = TRUE, block.size = 10)$importance[, 1]

  ## holdout VIMP for RFQ with and without blocking
  hvmp1 <- holdout.vimp(f, d, rfq =  TRUE,
               perf.type = "g.mean", block.size = 1)$importance[, 1]
  hvmp10 <- holdout.vimp(f, d, rfq =  TRUE,
               perf.type = "g.mean", block.size = 10)$importance[, 1]
  
  ## compare VIMP values
  imp <- 100 * cbind(vmp1, vmp10, hvmp1, hvmp10)
  legn <- c("vimp-1", "vimp-10","hvimp-1", "hvimp-10")
  colr <- rep(4,20+q)
  colr[1:20] <- 2
  ylim <- range(c(imp))
  nms <- 1:(20+q)
  par(mfrow=c(2,2))
  barplot(imp[,1],col=colr,las=2,main=legn[1],ylim=ylim,names.arg=nms)
  barplot(imp[,2],col=colr,las=2,main=legn[2],ylim=ylim,names.arg=nms)
  barplot(imp[,3],col=colr,las=2,main=legn[3],ylim=ylim,names.arg=nms)
  barplot(imp[,4],col=colr,las=2,main=legn[4],ylim=ylim,names.arg=nms)

}

## ------------------------------------------------------------
## multivariate regression analysis
## ------------------------------------------------------------
mtcars.mreg <- holdout.vimp(Multivar(mpg, cyl) ~., data = mtcars,
                                    vtry = 3,
                                    block.size = 1,
                                    samptype = "swr",
                                    sampsize = dim(mtcars)[1])
print(mtcars.mreg$importance)

## ------------------------------------------------------------
## mixed outcomes analysis
## ------------------------------------------------------------

mtcars.new <- mtcars
mtcars.new$cyl <- factor(mtcars.new$cyl)
mtcars.new$carb <- factor(mtcars.new$carb, ordered = TRUE)
mtcars.mix <- holdout.vimp(cbind(carb, mpg, cyl) ~., data = mtcars.new,
                                   ntree = 100,
                                   block.size = 2,
                                   vtry = 1)
print(mtcars.mix$importance)

##------------------------------------------------------------
## survival analysis
##------------------------------------------------------------

## Primary biliary cirrhosis (PBC) of the liver
data(pbc, package = "randomForestSRC")
pbc.obj <- holdout.vimp(Surv(days, status) ~ ., pbc,
                                nsplit = 10,
                                ntree = 1000,
                                na.action = "na.impute")
print(pbc.obj$importance)

##------------------------------------------------------------
## competing risks
##------------------------------------------------------------

## WIHS analysis
## cumulative incidence function (CIF) for HAART and AIDS stratified by IDU

data(wihs, package = "randomForestSRC")
wihs.obj <- holdout.vimp(Surv(time, status) ~ ., wihs,
                                 nsplit = 3,
                                 ntree = 100)
print(wihs.obj$importance)

---

Ames Iowa Housing Data

Description

Data from the Ames Assessor's Office used in assessing values of individual residential properties sold in Ames, Iowa from 2006 to 2010. This is a regression problem and the goal is to predict "SalePrice" which records the price of a home in thousands of dollars.

References

De Cock, D., (2011). Ames, Iowa: Alternative to the Boston housing data as an end of semester regression project. Journal of Statistics Education, 19(3), 1–14.

Examples

## load the data
data(housing, package = "randomForestSRC")

## the original data contains lots of missing data, so impute it
## use missForest, can be slow so grow trees with small training sizes
housing2 <- impute(data = housing, mf.q = 1, sampsize = function(x){x * .1})

## same idea ... but directly use rfsrc.fast and multivariate missForest 
housing3 <- impute(data = housing, mf.q = .5, fast = TRUE)

## even faster, but potentially less acurate
housing4 <- impute(SalePrice~., housing, splitrule = "random", nimpute = 1)

---

Random Forests for Imbalanced Two-Class Classification

Description

Fits a random forest for a two-class outcome using the random forest quantile classifier (RFQ), balanced random forest (BRF), or standard classification. RFQ changes the probability threshold used to assign class labels; BRF changes the sampling used to grow the trees. The splitting rule, performance measure, and variable importance can be specified separately. RFQ and BRF use normalized Brier performance by default for prediction error and variable importance.

Companion helpers summarize performance, select a probability threshold, or extract the class-frequency threshold from response labels.

Usage

## S3 method for class 'rfsrc'
imbalanced(formula, data, ntree = 3000, 
  method = c("rfq", "brf", "standard"), splitrule = "auc",
  perf.type = NULL, block.size = NULL, fast = FALSE,
  ratio = NULL, ...)

get.imbalanced.performance(obj, prob = NULL, threshold = NULL,
                           confusion = FALSE, robust = FALSE)
get.imbalanced.optimize(obj, prob = NULL, newdata = NULL,
                        measure = c("gmean", "F1", "F1mod", "F1modgmean"),
                        ngrid = 1000, plot.it = TRUE)
get.rfq.threshold(y)

## S3 method for class 'imbalanced.performance'
print(x, digits = 4,
                                      show.confusion = TRUE, ...)

Arguments

formula	A model formula with one factor response having two levels and one or more predictors.
data	A data frame containing the response and predictors. Data-frame subclasses are converted to a plain data frame.
ntree	Number of trees to grow. The default is 3000.
method	Classification method. "rfq", the default, uses the RFQ decision rule of O'Brien and Ishwaran (2019). "brf" uses class-weighted sampling with replacement and a reduced tree sample size. "standard" uses random forest classification without enabling the RFQ rule. See Details.
splitrule	Splitting rule passed to the forest grow function. The default is "auc". Other choices include "gini" and "entropy"; see rfsrc for supported classification rules. This setting is separate from perf.type.
perf.type	Performance measure used for the forest's prediction error (err.rate), variable importance (VIMP), and downstream inference based on these quantities. With NULL, the default is "brier" (normalized Brier score) for RFQ and BRF, and "default" (misclassification error) for standard classification. Set perf.type = "gmean" to use the previous RFQ/BRF default. Choices accepted by this wrapper are "none", "default", "standard", "misclass", "brier", "gmean", and the legacy spelling "g.mean". The latter is converted to "gmean". Use "none" to omit forest performance calculations.
block.size	Block size passed to the selected forest grow function for performance and VIMP calculations; see rfsrc. This wrapper passes NULL explicitly when the argument is not set. Specify the same block size when comparing VIMP across methods.
fast	Logical. Use rfsrc.fast for the RFQ or standard method when TRUE. The BRF branch always calls rfsrc and does not use this argument.
ratio	Optional, experimental majority-class undersampling for the RFQ and standard methods. Supply a number between zero and one. The majority class is sampled without replacement, while the minority class is bootstrapped with replacement. Small requested fractions may be increased by the sampling helper. This is a majority-class sampling fraction, rather than a requested final minority prevalence. Ignored for BRF. See Details for OOB support.
obj	For the performance and optimization helpers, a grow or prediction object with a single two-class outcome and observed responses, or a two-level response factor when prob is supplied.
prob	Optional probabilities aligned with the responses in obj. A numeric vector or one-column matrix or data frame contains minority-class probabilities. A two-column matrix or data frame contains both class probabilities; named columns are matched to factor levels, and unnamed columns follow factor-level order. NULL uses the object's predictions as described in Details.
threshold	Cutoff for minority-class probabilities in get.imbalanced.performance. NULL, the default, uses the minority proportion among the observed evaluation responses. "auto" uses that prevalence for an RFQ object and $0.5$ otherwise. "bayes" always uses $0.5$. A single finite numeric value specifies the cutoff directly; values outside $[0,1]$ are clipped to that interval. With direct response and probability inputs, "auto" uses $0.5$ because no fitted classifier is supplied.
confusion	Return a list including the confusion matrix instead of the default named vector of performance measures?
robust	Add one to the diagonal counts when calculating the four classification rates and their derived summaries? See Details.
newdata	Optional threshold-selection data for get.imbalanced.optimize, supplied with a grow object and prob = NULL. Must contain the predictors and observed response. NULL uses predictions already in obj. Ignored when prob is supplied.
measure	Performance measure to maximize during threshold selection: "gmean" (the default), "F1", "F1mod", or "F1modgmean". See Details for definitions.
ngrid	Number of equally spaced thresholds from zero to one, including both endpoints. Must be an integer of at least two.
plot.it	Draw the four threshold-search diagnostic panels?
y	Class-label vector for get.rfq.threshold, typically a two-level factor. Missing responses are omitted when counting classes.
...	Additional arguments passed to rfsrc or, when selected, rfsrc.fast. Examples include mtry, nodesize, nsplit, importance, and na.action. BRF sets its own case.wt, sampsize, and samptype. When ratio is supplied, the wrapper constructs samp and sets bootstrap = "by.user". Additional arguments to the print method are currently unused.
x	For the print method, the result of get.imbalanced.performance.
digits	Number of decimal places used in the performance display.
show.confusion	Print model and expected random-reference confusion tables when a confusion matrix is present in x? Obtain that matrix with confusion = TRUE in get.imbalanced.performance. This option does not recompute it.

Details

Choice of classifier

Let $\pi$ be the training proportion of the minority class and $\widehat p(x)$ its predicted probability. RFQ assigns an observation to the minority class when $\widehat p(x) \geq \pi$. The class-frequency threshold can identify minority observations whose predicted probability is below one half. The standard method instead uses class assignment based on the largest predicted probability. RFQ here is a classification rule, distinct from the continuous-response quantile estimation performed by quantreg.

BRF calls the grow function with replacement sampling, a tree sample size of twice the smaller observed class count, and case weights that give the two classes equal total sampling weight. The wrapper does not construct separate fixed-size samples from each class. Both fast and ratio are ignored in this branch.

Performance and variable importance

The default performance measure for RFQ and BRF is the normalized Brier score, selected by perf.type = "brier". For binary outcomes $Y_i \in \{0,1\}$, with $Y_i = 1$ denoting the minority class, and predicted minority-class probabilities $\widehat p_i$, the overall score is

$\mathrm{BS}^{*} = \frac{4}{n}\sum_{i=1}^{n} (Y_i - \widehat p_i)^2.$

The mean is taken over the observations being scored. Smaller values are better. This is four times the binary Brier score; the normalization does not reweight the classes or change VIMP rankings.

Brier performance responds to changes in predicted probabilities even when the assigned classes stay the same. This can improve the separation of predictive variables from noise in difficult imbalanced problems. Request permutation VIMP with importance = "permute"; it measures the increase in the selected prediction error after a variable is permuted.

G-mean remains available with perf.type = "gmean". It is the geometric mean of minority sensitivity and majority specificity, $\sqrt{\mathrm{sensitivity}\, \mathrm{specificity}}$. Larger G-mean values are better; the forest uses one minus G-mean as overall prediction error. Unlike Brier performance, it depends on the class assignments. Changing perf.type changes the performance measure, while method selects the classifier and splitrule selects the splitting criterion.

The selected metric is used for err.rate and VIMP and is inherited by prediction and VIMP calculations unless overridden. subsample also uses the fitted forest's metric for inference on VIMP and prediction error. The helper get.imbalanced.performance continues to report all its performance measures, regardless of perf.type. Selecting a classification threshold with get.imbalanced.optimize remains a separate operation; changing that threshold does not change a Brier score calculated from the same probabilities.

For comparisons, use a common perf.type, splitting rule, and block.size, unless the comparison is intended to study those settings. To compare performance metrics, vimp can score the same fitted forest under each metric. Compare recovery or ranking of predictive variables, rather than raw VIMP magnitudes on different scales. Increase ntree as needed to assess the stability of OOB performance and VIMP. The simulation example compares Brier and G-mean VIMP using known signal and noise variables.

Sampling and missing values

With ratio = NULL, the RFQ and standard branches pass missing-data handling to the selected grow function. BRF and calls using ratio first apply complete-case omission to the supplied data frame. This omission currently includes columns that are not selected by the formula and occurs before the grow function handles na.action.

With ratio supplied, each tree uses a bootstrap sample of the minority observations and a without-replacement sample of majority observations. At ratio = 1, every majority observation is in-bag for every tree, so majority-class OOB performance is unavailable. A single minority observation is also always in-bag in this sampling scheme. Use a separate evaluation sample when OOB support is absent.

Detailed performance summaries

get.imbalanced.performance calculates performance from observed responses and predicted probabilities.

Probability inputs.

Supply a grow or prediction object, or a two-level response factor with prob. With prob = NULL, the helper uses predicted.oob when present and not entirely missing; otherwise it uses predicted. When only some OOB predictions are missing, those entries remain missing.

A numeric vector or one-column matrix or data frame supplies minority-class probabilities; no column name is required. Two-column inputs supply both class probabilities; named columns are matched to factor levels, and unnamed columns follow factor-level order. Rows must align with the responses.

Evaluation classes.

All observed evaluation responses determine class counts, the minority class, and its proportion, before excluding unavailable predictions. The first factor level is the minority when counts tie.

Probability checks.

Finite probabilities below zero are set to zero, and those above one are set to one. Finite two-column rows must sum to one within numerical tolerance. Scoring excludes rows with missing responses or nonfinite probabilities; both columns must be finite for two-column inputs. Unavailable scores are NA.

Classification threshold.

An observation is assigned to the minority class when its predicted minority-class probability is at least the threshold, and to the majority class otherwise. The threshold choices are:

NULL

The evaluation minority proportion for every method (default).

"auto"

The evaluation minority proportion when forest$rfq is TRUE; $0.5$ otherwise, including direct response/probability inputs.

"bayes"

$0.5$ for every method.

Numeric

A single finite cutoff. Values below zero are set to zero, and values above one are set to one.

For RFQ test predictions, NULL and "auto" use evaluation prevalence, not training prevalence. The RFQ flag, rather than splitrule or perf.type, determines the "auto" choice. Changing the threshold preserves the evaluation minority definition and the fitted forest.

Performance measures.

The summary includes sensitivity (sens), specificity (spec), precision (prec), negative predictive value (npv), G-mean, misclass, brier, brier.norm, auc, logloss, and pr.auc. For complementary binary probabilities, brier is the mean squared error of the minority probability; brier.norm is four times this value. Log loss excludes infinite losses without clipping zero probabilities.

F-measures.

F1 is the harmonic mean of precision and sensitivity; F1mod is the harmonic mean of all four classification rates. F1gmean averages F1 and G-mean; F1modgmean averages F1mod and G-mean. A harmonic mean is zero if any component rate is zero and all are defined; it is NA if any required rate is undefined.

Robust summaries.

robust = TRUE adds one to each diagonal confusion count when calculating the four rates and their derived summaries. Returned confusion counts, misclassification error, and probability-based scores are unchanged.

Confusion matrix.

confusion = TRUE includes the confusion matrix: rows are observed classes and columns are predicted classes, with 0 denoting the majority and 1 the minority.

Random-score reference

The performance summary has a "rand" attribute containing reference metrics and expected confusion counts. This reference uses a score $U$ uniformly distributed on $[0,1]$, independent of class, and predicts the minority class when $U$ reaches the chosen threshold. For threshold $t$, its sensitivity is $1-t$, specificity is $t$, and ROC AUC is $1/2$. The separately returned PR-AUC reference is the proportion of positive observations among the rows scored by the PR helper. Its unnormalized Brier, normalized Brier, and log-loss references are $1/3$, $4/3$, and $1$, respectively. These differ from a predictor that always returns probability $1/2$.

Nonlinear rate summaries are computed from the reference rates (or smoothed expected counts under robust = TRUE); they are not finite-sample expectations of those nonlinear statistics. The performance print method compares the model and reference with Delta and percentage Gain, oriented so that positive values indicate improvement. The fitted-forest printer has its own probability reference; it need not use the same baseline.

Threshold selection

get.imbalanced.optimize searches ngrid equally spaced thresholds from zero to one, using G-mean, F1, F1mod, or F1modgmean as selected by measure. It uses the first grid point attaining the largest available value. It returns one row of performance metrics at the selected threshold and does not change the forest or its stored class predictions. plot.it = FALSE suppresses its four diagnostic panels. ngrid must be an integer of at least two. When the chosen measure is unavailable at every grid point, the function stops with an explanatory error. A valid optimum of zero is retained. The caller's graphics settings are restored after plotting.

A grow object ordinarily supplies OOB probabilities for this search. Supplying newdata searches using that data's observed responses and predictions, making it a threshold-selection sample. For an independent assessment of the chosen threshold, apply it unchanged to a separate test sample. The examples select a threshold from training OOB predictions and carry it forward to the test data.

Printing performance

Printing a performance result displays two metric tables. The first contains classification rates and their F1 and G-mean combinations; the second contains overall performance, including G-mean, AUC, PR-AUC, misclassification, Brier scores, and log loss. The same G-mean row appears in both tables; F1-related combinations appear only in the first.

Delta and percentage Gain compare the model with its random-score reference, with positive values indicating improvement. Larger rates, G-mean, and AUC values are better; smaller errors and losses are better. digits controls the display precision, and show.confusion = FALSE omits the optional confusion tables. These settings do not change the stored performance values.

Value

imbalanced returns an object from the selected forest grow function. Usual components include the class-probability matrices predicted and predicted.oob, class assignments class and class.oob, prediction error err.rate, and requested VIMP importance. Prediction error and VIMP use the selected perf.type; for RFQ and BRF the default is normalized Brier. The first column of err.rate and importance contains the overall error and VIMP, respectively. Availability depends on the selected grow interface and sampling options; see rfsrc and rfsrc.fast.

Performance and threshold helpers

get.imbalanced.performance

A named numeric vector of class counts, imbalance ratio, resolved numeric threshold, and performance measures. With confusion = TRUE, a named list also includes confusion, the confusion matrix with a class.error column. Both forms have class "imbalanced.performance" and a "rand" attribute containing the random-reference metrics and expected confusion counts.

get.imbalanced.optimize

A one-row data frame of performance measures at the selected threshold. Extract the cutoff from the threshold column; the row name is its grid index.

get.rfq.threshold

The smaller class count divided by the total observed class count, as a numeric scalar. Returns NULL unless the class-frequency table has exactly two entries. No forest or probability predictions are required.

The performance and optimization helpers return NULL when the evaluation response is not a two-level factor.

print.imbalanced.performance displays the summaries and returns x invisibly.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Chen, C., Liaw, A. and Breiman, L. (2004). Using random forest to learn imbalanced data. University of California, Berkeley, Technical Report 110.

Kubat, M., Holte, R. and Matwin, S. (1997). Learning when negative examples abound. Machine Learning, ECML-97: 146-153.

O'Brien, R. and Ishwaran, H. (2019). A random forests quantile classifier for class imbalanced data. Pattern Recognition, 90, 232-249.

See Also

rfsrc, rfsrc.fast, get.imbalanced.performance, vimp, subsample

Examples

## ------------------------------------------------------------
## RFQ, BRF, and standard classification
## ------------------------------------------------------------
data(breast, package = "randomForestSRC")
dta <- na.omit(breast)
f <- status ~ .

## Use 300 trees for illustration; imbalanced() defaults to 3000.
set.seed(19)
rfq <- imbalanced(f, dta, ntree = 300)
print(rfq)
print(get.imbalanced.performance(rfq, confusion = TRUE))

## Equivalent direct RFQ call using rfsrc().
## RFQ uses equal case weights and sampling without replacement.
set.seed(19)
rfq.direct <- rfsrc(f, dta, ntree = 300,
                    rfq = TRUE, splitrule = "auc", perf.type = "brier",
                    block.size = NULL, bootstrap = "by.root",
                    samptype = "swor", sampsize = function(n) n * .632,
                    case.wt = NULL)
print(all.equal(rfq$predicted.oob, rfq.direct$predicted.oob))

## Extract the class-frequency threshold from the observed responses.
print(get.rfq.threshold(rfq$yvar))

set.seed(20)
brf <- imbalanced(f, dta, ntree = 300, method = "brf")

## Equivalent direct BRF call using rfsrc().
## Each class has equal total sampling weight; sampling is with replacement.
class.frequency <- table(dta$status)
case.wt <- min(class.frequency) /
  as.numeric(class.frequency[as.character(dta$status)])
set.seed(20)
brf.direct <- rfsrc(f, dta, ntree = 300,
                    rfq = FALSE, splitrule = "auc", perf.type = "brier",
                    block.size = NULL, bootstrap = "by.root",
                    samptype = "swr", sampsize = 2L * min(class.frequency),
                    case.wt = case.wt)
print(all.equal(brf$predicted.oob, brf.direct$predicted.oob))

standard <- imbalanced(f, dta, ntree = 300, method = "standard")

## Compare all three using the helper's default evaluation rule.
metrics <- c("sens", "spec", "gmean", "misclass", "brier.norm", "auc")
comparison <- rbind(
  RFQ = as.numeric(get.imbalanced.performance(rfq)[metrics]),
  BRF = as.numeric(get.imbalanced.performance(brf)[metrics]),
  Standard = as.numeric(get.imbalanced.performance(standard)[metrics])
)
colnames(comparison) <- metrics
print(comparison)

## Supply a cutoff explicitly to inspect another probability rule.
print(get.imbalanced.performance(standard, threshold = .5))

## Let the saved classifier setting select prevalence or 0.5.
print(get.imbalanced.performance(rfq, threshold = "auto"))
print(get.imbalanced.performance(standard, threshold = "auto"))
print(get.imbalanced.performance(brf, threshold = "auto"))

## Request the 0.5 Bayes cutoff explicitly, including for RFQ.
print(get.imbalanced.performance(rfq, threshold = "bayes"))

## Vectors, one-column objects, and named two-column matrices.
class.frequency <- table(rfq$yvar)
minority <- names(class.frequency)[which.min(class.frequency)]
p <- rfq$predicted.oob[, minority]
print(get.imbalanced.performance(rfq$yvar, p))
print(get.imbalanced.performance(rfq$yvar, data.frame(p = p)))
print(get.imbalanced.performance(
  rfq$yvar, rfq$predicted.oob[, rev(levels(rfq$yvar)), drop = FALSE]
))

## ------------------------------------------------------------
## Train/test assessment with a training-selected threshold
## ------------------------------------------------------------
set.seed(23)
trn.id <- sort(unlist(lapply(
  split(seq_len(nrow(dta)), dta$status),
  function(idx) idx[sample.int(length(idx),
    size = floor(2 * length(idx) / 3))]
)))
trn <- dta[trn.id, , drop = FALSE]
tst <- dta[-trn.id, , drop = FALSE]
print(rbind(training = table(trn$status), test = table(tst$status)))

fit <- imbalanced(f, trn, ntree = 300)
selected <- get.imbalanced.optimize(fit, measure = "gmean",
                                     ngrid = 501, plot.it = FALSE)
print(selected)
threshold.selected <- selected$threshold

pred <- predict(fit, newdata = tst)
print(get.imbalanced.performance(pred))
print(get.imbalanced.performance(pred, threshold = threshold.selected))

## Predicting labels for new rows does not require their outcomes.
class.frequency <- table(trn$status)
minority <- names(class.frequency)[which.min(class.frequency)]
majority <- setdiff(levels(trn$status), minority)
pred.x <- predict(fit, newdata = tst[, fit$xvar.names, drop = FALSE])
selected.class <- factor(
  ifelse(pred.x$predicted[, minority] >= threshold.selected,
         minority, majority),
  levels = levels(trn$status)
)
print(head(selected.class))

## ------------------------------------------------------------
## Explicit permutation VIMP and blocking
## ------------------------------------------------------------
## RFQ uses normalized Brier performance by default.
rfq.vimp <- imbalanced(f, dta, ntree = 300,
                       importance = "permute", block.size = 1)
rfq.block <- imbalanced(f, dta, ntree = 300,
                        importance = "permute", block.size = 10)
print(cbind(unblocked = rfq.vimp$importance[, 1],
            blocked = rfq.block$importance[, 1]))

## Small B is for illustration; increase it for final inference.
smp <- subsample(rfq.block, B = 25, verbose = FALSE)
print(extract.subsample(smp)$var.sel.Z)
plot.subsample(smp, cex.axis = .7, alpha = .5)

## ------------------------------------------------------------
## Strong imbalance: probability signals and permutation VIMP
## ------------------------------------------------------------
## x1 identifies a small subgroup; x2 changes risk within it.
## Both affect the outcome probability.  The other variables are noise.
## Expected minority prevalence is .04 * (.15 + .5 / 2) = .016.
## Within the subgroup, all true probabilities exceed this prevalence,
## so the population RFQ class rule uses x1 alone.  Brier performance
## can also detect the probability information supplied by x2.
local({
  ## One thread makes the blocked VIMP comparison reproducible.
  old.options <- options(rf.cores = 1)
  on.exit(options(old.options))

  set.seed(47)
  n <- 4000
  p <- 20
  x <- matrix(runif(n * p), ncol = p)
  x[, 1] <- as.numeric(x[, 1] < .04)
  colnames(x) <- paste0("x", seq_len(p))
  prob <- x[, 1] * (.15 + .5 * x[, 2])
  y <- factor(rbinom(n, 1, prob), levels = c(0, 1),
              labels = c("majority", "minority"))
  sim.data <- data.frame(y = y, x)

  ## Grow one RFQ forest, then use the same trees for both metrics.
  ## Brier is explicit here; it is also the RFQ default.
  sim.fit <- imbalanced(y ~ ., sim.data, ntree = 500, mtry = p,
                         splitrule = "auc", perf.type = "brier",
                         block.size = 50, seed = -103)
  brier <- vimp(sim.fit, importance = "permute", perf.type = "brier",
                 block.size = 50, seed = -107)
  gmean <- vimp(sim.fit, importance = "permute", perf.type = "gmean",
                 block.size = 50, seed = -107)

  counts <- table(y)
  cat("\nClass counts and majority/minority ratio\n")
  print(c(n.majority = unname(counts["majority"]),
          n.minority = unname(counts["minority"]),
          iratio = unname(counts["majority"] / counts["minority"])))

  ## Compare signals with noise within each metric's own scale.
  vimp.values <- cbind(Brier = brier$importance[, 1],
                       Gmean = gmean$importance[, 1])
  cat("\nPermutation VIMP: true signals and largest noise value\n")
  print(rbind(vimp.values[1:2, , drop = FALSE],
               largest.noise = apply(vimp.values[-(1:2), , drop = FALSE],
                                     2, max)), digits = 5)

  ## A smaller rank is better; the ideal ranks are 1 and 2.
  vimp.ranks <- apply(vimp.values, 2, function(z) {
    rank(-z, ties.method = "average")
  })
  cat("\nRanks of the two true signals among all 20 variables\n")
  print(vimp.ranks[1:2, , drop = FALSE])
})

---

Learn a predictive imputer for test-time imputation and OOD scoring

Description

Learns a predictive imputer from training data for later use on new data.

If the training data contain missing values, the function first imputes them using impute. It then fits one saved full-sweep learner per selected target on the completed training data and reuses those learners later to update missing values in new data without refitting on the test set.

The same saved learner bank can also be used to score new data for out-of-distribution (OOD) behavior. Note that OOD scores are available even when new data have missing values. Each selected target is reconstructed from its saved conditional learner and compared with the observed value. Target-wise discrepancies are calibrated against a training reference calculated from out-of-bag predictions computed during training.

If the training data are complete and target.mode = "all", the initial training-data imputation step is skipped and the full-sweep learners are fit directly from the complete training data.

If supervised.formula is supplied, the function also fits an internal supervised forest from the training data. The supervised forest is fit after training imputation and provides auxiliary predictors for test-time imputation and OOD scoring. The auxiliary predictors use supervised information learned from the training outcomes. Supervised outcomes supplied with new data are dropped and are not used at deployment time. Leave supervised.formula unspecified to use the learned imputer without these auxiliary predictors.

Usage

impute.learn.rfsrc(formula, data,
  ntree = 100, nodesize = 1, nsplit = 10,
  nimpute = 2, fast = FALSE, blocks,
  mf.q, max.iter = 10, eps = 0.01,
  ytry = NULL, always.use = NULL, verbose = TRUE,
  ...,
  supervised.formula = NULL,
  supervised.args = list(),
  full.sweep.options = list(ntree = 100, nsplit = 10),
  target.mode = c("missing.only", "all"),
  deployment.xvars = NULL,
  anonymous = TRUE,
  learner.prefix = "impute.learner.",
  learner.root = "learners",
  out.dir = NULL,
  wipe = TRUE,
  keep.models = is.null(out.dir),
  keep.ximp = FALSE,
  save.on.fit = !is.null(out.dir),
  save.ood = TRUE,
  weight = NULL,
  compress = 50,
  share.xvar = FALSE)

save.impute.learn.rfsrc(object, path, wipe = TRUE, verbose = TRUE,
  compress = 50, share.xvar = NULL)

load.impute.learn.rfsrc(path, targets = NULL, lazy = TRUE, verbose = TRUE)

## S3 method for class 'impute.learn.rfsrc'
predict(object, newdata,
  max.predict.iter = 3L,
  eps = 1e-3,
  targets = NULL,
  restore.integer = TRUE,
  cache.learners = c("session", "none", "all"),
  verbose = TRUE,
  ...)

impute.ood.rfsrc(object, newdata,
  targets = NULL,
  max.predict.iter = 3L,
  eps = 1e-3,
  cache.learners = c("all", "session", "none"),
  weight = NULL,
  aggregate = c("bounded.product",  "weighted.mean",
    "weighted.lp", "weighted.lp.log", "top.k"),
  aggregate.args = list(),
  return.details = FALSE,
  return.reconstruction = FALSE,
  verbose = TRUE,
  ...)

Arguments

formula	An optional symbolic model description passed to impute for the initial training-data imputation when supervised.formula is not supplied. It follows the same guidelines as in impute. This formula does not specify the full-sweep learner bank, the OOD targets, or the supervised auxiliary forest. When supervised.formula is supplied, the raw predictor block is taken from the right-hand side of supervised.formula, and formula is not used.
data	Training data, converted to a plain data frame before processing. Matrices, tibbles, and data.table objects can be supplied. Column names must be unique and nonempty, each column must be a vector, and numeric values must be finite or missing. Bare vectors and nested matrix or list columns are not supported. Variables that are not real-valued are coerced to factors before fitting when possible; otherwise fitting stops with an error. Rows and columns that are entirely missing are dropped before training begins. If supervised.formula is supplied, data should contain both the supervised response columns and the raw predictor columns. The learned imputer is then built on the raw predictor block defined by the right-hand side of supervised.formula.
ntree, nodesize, nsplit, nimpute, fast, blocks, max.iter, ytry, always.use, verbose	Arguments passed to impute for the initial training-data imputation. The argument full.sweep is controlled internally and should not be supplied. fast is also forwarded to the saved target learners and, unless overridden by supervised.args, to the supervised forest. verbose controls progress throughout training, prediction, and scoring. The same verbose, max.predict.iter, eps, and cache.learners controls are also used by predict.impute.learn and impute.ood. Iteration and forest-size controls must be scalar integers in their permitted ranges; max.iter must be at least one.
mf.q	Controls the imputation engine used by impute. With mf.q = 1 and always.use = NULL, targets are updated one at a time using missForest. Other positive settings use the multivariate missForest generalization; a fraction below one controls the proportion of missing-data variables grouped as responses, and a value above one specifies their requested group size. A non-NULL always.use selects the multivariate branch, including an empty vector or a vector with no matching column names. If mf.q is omitted, training uses on-the-fly imputation. The selected method is recorded in manifest$train.imputation.
eps	Finite nonnegative convergence threshold. In impute.learn this controls the initial missForest imputation. In predict.impute.learn and impute.ood it controls early stopping for the prediction-time sweep.
...	For impute.learn, additional arguments passed to impute. For predict.impute.learn and impute.ood, additional arguments are currently ignored.
supervised.formula	Optional supervised learning formula used to augment the learned imputer for improved OOD detection in supervised settings. The left-hand side defines the supervised response and the right-hand side defines the raw predictor block to be learned by impute.learn. The supervised forest is fit internally after the raw predictor block has been completed. Its training out-of-bag predicted values and test-time predicted values are appended internally as auxiliary predictors. Training auxiliary values use out-of-bag predictions when available, with ordinary forest predictions and training means used as fallbacks. At deployment, auxiliary values are computed once from the initialized raw predictors and held fixed during the imputation passes. Predictors must be raw column names; non-syntactic names can be enclosed in backticks. Interactions and transformations are not supported. These auxiliary predictors are internal and are not expected in newdata.
supervised.args	Optional named list of arguments passed to the internal supervised rfsrc fit. Entries for formula, data, and forest are controlled internally and are ignored if supplied. All other entries override the defaults used for this supervised fit, including defaults inherited from top-level arguments such as fast and internally supplied defaults such as perf.type = "none".
full.sweep.options	A named list of options used when fitting saved target learners after initial training-data imputation. Defaults are ntree = 100, nodesize = NULL, and nsplit = 10, independently of the corresponding settings used for initial imputation. Recognized entries include ntree, nodesize, nsplit, mtry, splitrule, bootstrap, sampsize, samptype, perf.type, rfq, save.memory, importance, proximity, terminal.qualts, and terminal.quants. To request on-the-fly prediction for the saved target learners, supply terminal.qualts = FALSE and terminal.quants = FALSE in this list. These settings are passed to the fitting function selected by anonymous; they do not change the options for initial imputation or for the supervised auxiliary forest, which uses supervised.args. Unknown entries are ignored with a warning. Duplicate names produce a warning, and the last entry for each name is used.
target.mode	Determines which raw variables receive a saved full-sweep learner. The default "missing.only" saves learners only for variables that were missing in the training data. The option "all" saves a learner for every retained raw variable. If the training data are complete, target.mode = "all" must be used. When supervised.formula is supplied, the default is promoted internally to "all" so that every retained raw predictor can later be updated and scored. For the broadest OOD coverage, target.mode = "all" is recommended so every deployment-time variable can be reconstructed.
deployment.xvars	Controls which raw predictors are assumed to be available later when the saved imputer is used on new data. If NULL, all raw non-target columns are used. If a character vector, the same raw predictor set is used for all targets. If a named list, names must be target variables and each entry gives the raw predictor set for that target; targets omitted from the list fall back to all raw non-target columns. Unnamed or duplicated list entries are not allowed. Entries named for non-target variables, and predictor names not found in the training data, are ignored with a warning. When supervised.formula is supplied, the internal auxiliary predicted.* columns are added automatically on top of these raw predictor sets. In other words, deployment.xvars restricts the raw predictors only.
anonymous	If TRUE, uses rfsrc.anonymous when fitting the saved target-wise full-sweep learner bank. The internal supervised forest, when requested through supervised.formula, is fit with rfsrc because it must be retained for later prediction-time auxiliary variables. Thus anonymous usually reduces the size of the target-wise learner bank, but not necessarily the size of the supervised forest.
learner.prefix, learner.root	Names used when writing saved full-sweep learners to disk. If supervised.formula is supplied, the internal supervised forest is saved under the same learner root. learner.root must be a relative directory path within the imputer directory. learner.prefix must be a single file-name component. Absolute paths and parent-directory components are not allowed. The top-level name ‘manifest.rds’ is reserved.
out.dir	Optional output directory. If supplied and save.on.fit = TRUE, the manifest and the saved full-sweep learners are written to this directory during fitting. This requires the fst package because learner components are serialized in FST format. If supervised.formula is supplied, the internal supervised forest is also saved there.
wipe	If TRUE, replaces an existing saved-imputer directory after the new imputer has been written and verified. If FALSE, other files already inside that directory are retained. Saving uses a temporary bundle and requires additional disk space.
keep.models	If TRUE, keeps the fitted full-sweep learners in memory in the returned object. At least one storage mode must be enabled: either keep.models = TRUE or out.dir with save.on.fit = TRUE. If supervised.formula is supplied, the internal supervised forest is kept in memory under the same rule.
keep.ximp	If TRUE, keeps the completed raw training predictor table in the returned object. This is not required for later prediction.
save.on.fit	If TRUE and out.dir is supplied, writes the imputer to disk during fitting.
compress	A single number from 0 to 100 controlling lossless FST compression of saved learners. The default, 50, balances file size and saving speed. Use 100 for stronger compression at the cost of slower saving. Applies to the target learners and the optional supervised forest. In impute.learn.rfsrc, this controls saving during fitting; in save.impute.learn.rfsrc, it controls the current save call and defaults to 50 independently of any earlier setting. It does not change forest fitting or prediction.
save.ood	If TRUE, computes and stores an OOD reference. The reference is built during training from target-wise out-of-bag reconstruction discrepancies, their target-wise calibrated training scores, and a default row-level weighted mean using the saved OOD target weights. If no weight is supplied at fit time, equal target weights are used. The saved target-wise training-score matrix allows impute.ood to rebuild a calibrated row-level percentile for arbitrary target subsets, test-time weight overrides, and alternate row aggregates besides the weighted mean.
share.xvar	Share the completed training predictor table across non-anonymous target learners instead of writing a separate ‘xvar.rda’ for each learner. The fitting default is FALSE. Set share.xvar = TRUE and anonymous = FALSE when fitting to capture the common table and exact learner row/column mappings. The option applies to saving during fitting and to later explicit saving. For save.impute.learn, NULL uses the object's shared-storage setting; FALSE writes self-contained learners. An older bank fitted without shared mappings cannot be converted by setting this argument alone. See Details.
object	An object returned by impute.learn or load.impute.learn.
path	For save.impute.learn, the destination directory. For load.impute.learn, the directory containing the saved imputer. Relative paths are resolved from the current working directory. Saving restrictions are described in Details; they do not apply to loading. Both operations require the fst package. Use load.impute.learn to load the complete imputer, including its shared predictors when present.
targets	Optional subset of target variables to load, update, or score. Unknown names are ignored with a warning. For impute.ood, row-level percentile calibration is rebuilt for the requested target subset from the saved target-wise training OOD scores whenever those scores are available in the manifest. A load-time subset limits which learners are available for predictor completion. It can therefore differ from scoring a subset of a fully loaded bank when other predictors are missing.
lazy	If TRUE, saved learners are loaded only when they are needed. If FALSE, all saved learners are loaded at once. If supervised.formula was used at fit time, the internal supervised forest follows the same lazy versus eager loading rule.
newdata	New data to be imputed or scored, converted to a plain data frame before processing. The column-name and vector-column requirements for data also apply. Missing columns are added and extra columns are dropped to match the training schema. A zero-row table returns correctly typed empty output without calling forests. Retained numeric values must be finite or missing. Nonmissing values that cannot be converted to a required numeric type become missing with a warning; their row indices are recorded in conversion.issues. Unseen factor levels are converted to NA for harmonization, but they are also tracked row-wise. In supervised mode, newdata should contain only the raw predictor columns learned by the imputer. The internal auxiliary columns are created automatically and should not be supplied by the user. If supervised response columns are supplied in newdata and are not part of the learned raw predictor block, they are treated as extra columns and dropped; they are not used to compute auxiliary predictors, imputations, or OOD scores. In impute.ood, observed raw target values are used to compute target-wise discrepancies. Missing raw target values are completed for predictor-side reconstruction but do not themselves contribute a finite target-wise OOD score. Any row containing an unseen factor level is flagged and its row-level OOD score is set to the maximum value.
max.predict.iter	Maximum number of full-sweep passes applied to newdata before the saved learner bank is used for OOD reconstruction or returned prediction-time imputations. Must be a nonnegative integer. Zero performs initialization without iterative target updates; OOD reconstruction is still performed.
restore.integer	If TRUE, generated imputations in integer-supported variables are rounded to the integer grid. This includes columns stored as R integers and numeric columns whose observed finite training values were all integer-valued. Observed numeric values in newdata are preserved, including fractional values, both as predictors during imputation and in the returned data. Restoration applies to initialization values and model-based updates. Numeric training columns retain numeric storage. Columns stored as R integers during training are returned as integer vectors only when every nonmissing completed value is exactly integral and within R's integer range; otherwise numeric storage preserves observed fractions and out-of-range values. Set to FALSE to retain unrounded generated numeric values. Factor columns are always conformed back to the training schema. The package operates on real-valued and factor variables; inputs that are not real-valued are coerced to factors during preprocessing when possible, otherwise an error is raised.
cache.learners	How saved learners are reused during prediction or OOD scoring. For predict.impute.learn, the default "session" loads each needed learner once per call. The option "none" reloads a learner every time it is needed. The option "all" loads all requested learners before work starts. For impute.ood, "all" is the default because the saved learner bank is typically reused once for predictor-side completion and again for target reconstruction.
weight	Optional nonnegative target weights used for row-level OOD aggregation. In impute.learn, these weights define the default row-level OOD weighting scheme stored in the manifest. In impute.ood, they define the active row-level weighting scheme for the current scoring call. If supplied as a named vector, entries are matched to targets by name; omitted targets are set to zero, and extra names are ignored. If omitted in impute.ood, the saved training-time OOD weights are used automatically. Because the fit stores target-wise training OOD scores, score.percentile can be recalibrated automatically for test-time weight overrides rather than being limited to the original training-time weights.
aggregate	Row-level aggregation metric used by impute.ood to combine calibrated target-wise OOD scores. The default "bounded.product" applies a weighted product of the form $1-\prod_j \max(1-u_j,\varepsilon)^{\tilde w_j}$, where $u_j$ is the calibrated score for target $j$ and $\tilde w_j$ denotes weights normalized over the positively weighted, scoreable targets in the row. "weighted.mean" is the weighted average. "weighted.lp" applies a weighted Minkowski $L_p$ aggregation to the calibrated target scores. "weighted.lp.log" first applies the tail-stretching transform $-\log\{\max(1-u_j,\varepsilon)\}$ to each calibrated target score $u_j$ and then applies the weighted $L_p$ aggregation. "top.k" averages only the $k$ largest scoreable target scores among the positively weighted targets. When the fit stores target-wise training OOD scores, score.percentile is rebuilt for the requested aggregate as well as the requested targets and weights.
aggregate.args	Optional named list of tuning arguments for aggregate. Recognized entries are p for "weighted.lp" and "weighted.lp.log", k (or top.k) for "top.k", and eps for "weighted.lp.log" and "bounded.product". The default values are p = 2, k = 1, and eps = 1e-12. The power must be finite and at least one; the probability floor must satisfy $0 < \varepsilon < 1$. The top-$k$ value must be finite, at least one, and no larger than .Machine$integer.max; it is rounded to an integer. Unknown entries are ignored with a warning.
return.details	If TRUE, impute.ood returns the per-target discrepancy and calibrated target-score matrices and a richer info list.
return.reconstruction	If TRUE, impute.ood returns the target-wise reconstructed values used during OOD scoring. The reconstructed values are returned on the raw target scale, with one column per scored target. Typically these values are computed internally and discarded.

Details

Training

Training begins by converting variables that are not real-valued to factors when possible; otherwise fitting stops with an error. Rows and columns that are entirely missing are removed before the training schema is stored. The imputer is then fitted in two stages:

1. Complete the training data. Missing values are imputed using impute, with the same options as that function. With mf.q = 1 and always.use = NULL, targets are updated one at a time. Other positive settings use the multivariate missForest generalization.

   If mf.q is omitted, on-the-fly imputation is used. It is supervised when an initial-imputation formula is used, and unsupervised otherwise. Complete training data with target.mode = "all" skip this stage.

2. Fit the saved learners. A full forest sweep is fitted on the completed training data. For each target selected by target.mode, a forest is fitted using rows where that target was originally observed and predictors selected by deployment.xvars. All saved learners use the same completed training table; fitting them does not further update that table.

Training stops if every requested target learner fails. If some succeed, a partial bank is returned with one summary warning. manifest$learners records each target's status and error, and printing the imputer reports successful and unavailable learners.

Choosing predictors

By default, deployment.xvars = NULL uses every non-target column as a predictor. Restrict deployment.xvars when the training data include outcomes, future-only variables, identifiers, or other fields that could introduce leakage or will be unavailable in new data.

Supervised auxiliary predictors

Supplying supervised.formula fits an internal supervised forest. Its out-of-bag predicted values for training data and predicted values for new data provide auxiliary predictors for the saved learners. These variables can influence both imputation and OOD scoring, but are not imputation targets. Leave supervised.formula unspecified to use the basic unsupervised learned imputer.

The auxiliary variables are created automatically and added to each saved target's predictor set. deployment.xvars restricts only the raw predictors; users do not need to supply auxiliary columns in new data.

Ordinary forest predictions fill unavailable auxiliary predictions, including when the supervised forest has no OOB predictions, as with supervised.args = list(bootstrap = "none"). This fallback applies only to auxiliary predictors. OOD reference discrepancies always use the final target learners' OOB predictions.

Saving and loading

A saved imputer consists of a manifest and a directory of learners, with an optional shared predictor table. Learner components use the same serializers as fast.save and fast.load, so saving and loading require the fst package. The supervised forest, when used, is saved alongside the target learners. The explicit save method can save learners from memory or from the directory where the imputer is already saved.

Each imputer has its own directory, which can be inside a project directory; for example, file.path(getwd(), "imputer"). Saving accepts a new directory, an empty directory, or a directory recognized as a saved imputer by its manifest. Other nonempty directories are rejected with either setting of wipe.

Saving cannot replace the current working directory, home directory, filesystem root, or their ancestors. A complete imputer can be saved back to the same directory. A different destination must neither contain nor be inside the existing imputer directory. An object loaded with a target subset must be saved to a different directory. Loading has none of these location restrictions.

Set compress = 100 for stronger FST compression when saving a deployment bundle. Supply this directly to impute.learn or save.impute.learn, rather than in full.sweep.options. The manifest and remaining RDS metadata use their existing compression settings. Loading requires no compression argument.

Both training-time and explicit saving write and verify the new learners before replacing an existing destination. A failed staging operation leaves the previous saved imputer unchanged. If replacement fails, restoration of the previous directory is attempted; an unrecovered backup path is reported.

Shared training predictors

With share.xvar = TRUE, the completed training table used by the final target learners is saved once in ‘_shared_xvar/xvar.fst’. In supervised mode this table also includes the auxiliary training predictors actually used in those fits. Each target stores its exact training-row indices, predictor names and column encodings. These reconstruct the learner's training predictors before prediction. Reconstruction is checked against the fitted forest before its separate ‘xvar.rda’ is omitted. Trees, responses, terminal information and OOD references keep their existing representations. The single supervised forest remains self-contained.

Shared storage requires anonymous = FALSE. It is independent of compress, terminal.qualts and terminal.quants; their existing settings still apply. keep.ximp continues to control the returned raw completed table. An in-memory shared-capable fit separately holds the common predictor table needed for deferred saving. Saved shared banks and eagerly loaded shared objects can be saved again without reconstructing row identities from predictor values.

Load the whole imputer with load.impute.learn. Shared learner subdirectories are not independently portable, and direct fast.load calls on them stop with an explanatory message. Keep the manifest, shared table and learners together when moving or archiving the bank. A compatible shared-aware package version is required. The updated loader also accepts the existing independent-learner format. Saving a shared object with share.xvar = FALSE writes that independent format.

Lazy loading leaves the common table on disk until first use. Each prediction or OOD operation loads it at most once and reuses it across learner accesses, including when cache.learners = "none". Reconstructed learners remain normal forest objects: caching or eagerly loading many learners can still use memory for their individual predictor tables. Shared storage primarily reduces repeated serialization and on-disk duplication. Missing shared files, checksum mismatches and invalid mappings cause errors rather than silently substituting initialization fills.

Imputing new data

Prediction applies the saved learners in three steps:

1. Match the training schema. Columns and types in newdata are matched to the training schema. Supervised response columns outside the learned raw predictor block are dropped and are not used for prediction or OOD scoring.

2. Initialize missing values. Missing raw values are filled with training means or modes. In supervised mode, the saved supervised forest then computes auxiliary predictors from this initialized table. The auxiliary predictors stay fixed throughout the subsequent passes.

3. Update with saved learners. Full-sweep passes update missing values in the selected targets. Each target update requires one valid prediction per requested row. Failed or unavailable predictions leave the previous imputed values unchanged and are recorded in target.issues. A pass with no valid model updates is reported separately from convergence.

With target.mode = "missing.only", a variable that was complete in training but is missing in new data receives only an initialization value. Use target.mode = "all" when missing values may appear later in any raw variable. Complete training data also require this setting because there are no missing variables from which to select the saved targets.

Integer restoration applies only to generated values, including initialization values in columns without a target update. Observed numeric entries are unchanged.

OOD training references

With save.ood = TRUE, each saved learner's out-of-bag predictions are compared with observed training values to form target-wise reconstruction discrepancies:

• Continuous and integer targets use absolute reconstruction error.

• Factor targets use the negative log probability of the observed class when OOB class probabilities are available. If only OOB class labels are available, unordered factors use misclassification error and ordered factors use scaled distance between ordered levels.

Each target's discrepancy rule is saved with its reference and used unchanged for new data. For probability-based discrepancies, the saved probability floor is 1e-12; finite zero probabilities use this floor. Unavailable or invalid predictions give missing discrepancies, without switching to another rule. Each learner entry records n.oob.finite and n.oob.nonfinite.

OOB predictions come from the final target forests after training-data imputation and auxiliary predictor construction. Those earlier steps are not repeated for each OOB sample. Missing OOB predictions are never replaced by ordinary training predictions when constructing these references.

The training-time row reference combines calibrated target scores using a weighted mean. If weight is omitted, all saved OOD targets receive weight 1. Named weights are matched to targets, and omitted targets receive weight 0. These weights are also the defaults for later OOD scoring.

OOD scoring

impute.ood first completes the predictors in newdata using the same schema matching, initialization, and full-sweep passes as predict.impute.learn. It then reconstructs each requested raw target from its saved learner and compares the reconstruction with the observed value in newdata. A target that is missing in a row does not contribute to that row's OOD score.

Discrepancies are converted to target-wise OOD scores using the saved target-specific training references. Values strictly above a nonempty reference's maximum, including positive infinity, receive its largest stored probability. Ties, including equality at the largest quantile, follow the reference's quantile-grid convention. Missing discrepancies and empty references are unscored.

In supervised mode, the auxiliary variables are predictors for reconstruction; the supervised response is not used in new data.

Older saved imputers that do not record the factor discrepancy rule must be refitted before OOD scoring those targets. They remain usable for imputation and numeric-only OOD scoring.

Row-level OOD scores

impute.ood returns two row-level summaries:

score

Combines calibrated target scores over the targets that are observed and scoreable in each row. The default is a bounded product rule. Alternatives are a weighted mean, weighted $L_p$, log-tail weighted $L_p$, and top-$k$ rules. These options allow greater sensitivity to sparse but severe coordinate shifts. Unless overridden, scoring uses the weights saved by impute.learn.

score.percentile

Calibrates the row score against a training reference rebuilt from the saved target-wise training OOD scores. The reference uses the requested target subset, weights, and row aggregate, so percentile calibration remains available when any of these settings change. It is not conditioned on each new row's pattern of observed targets.

Unseen factor levels

Unseen factor levels are tracked by row when matching new data to the training schema. impute.ood flags these rows and assigns them the maximum row-level score. If an unseen level occurs in a scored target, its target-level discrepancy is also maximal.

Value

impute.learn returns an object of class c("impute.learn.rfsrc", "impute.learn"). The object contains a manifest, optionally the fitted full-sweep learners, optionally the internal supervised forest when supervised.formula is used, optionally the completed raw training predictor table, and optionally a path to the saved imputer on disk. If save.ood = TRUE, the manifest also contains an ood component storing compact target-wise OOD references, the saved row-by-target training OOD score matrix used for later percentile recalibration, and the default OOD aggregation weights. When supervised mode is active, the manifest also records the supervised family, response names, and the internally created auxiliary predicted.* variable names.

Shared-capable fits also contain manifest$shared.xvar, the common-table specification, and per-learner xvar.ref reconstruction mappings. shared.xvar.data holds the common table for in-memory fits with keep.models = TRUE and eagerly loaded shared objects. It is NULL for disk-backed lazy objects. This component is separate from ximp.

load.impute.learn returns an object of the same class.

predict.impute.learn returns a data frame with imputed values overlaid on the raw predictor table, retaining its row names. An attribute named "impute.learn.info" contains prediction-time diagnostics such as the number of sweep passes, pass-difference history, caching mode, disk-load counts, schema harmonization details, dropped supervised response columns when present, row-wise unseen-factor flags, supervised-auxiliary diagnostics when present, and any targets skipped because a learner was unavailable or a prediction failed. pass.updated.cells and pass.failed.cells count accepted updates and failed updates in each pass; converged and stopping.reason distinguish convergence from initialization-only, empty-input, iteration-limit, and failed-update stopping. conversion.issues records the row indices of nonmissing values that became missing during numeric conversion. In both prediction and OOD diagnostics, n.disk.loads counts successful target-learner load operations, including repeated loads with cache.learners = "none". disk.load.targets lists the distinct targets loaded. Supervised forest loading is reported separately in info$supervised and is not included in the target-learner count.

impute.ood returns an object of class c("impute.ood.rfsrc", "impute.ood"). It is a list with the following components:

• score: the row-level aggregate of calibrated target-wise OOD scores under the requested aggregate and weight. Larger values indicate greater out-of-distribution behavior. For aggregate = "weighted.lp.log", the probability floor bounds the mathematical aggregate by $-\log(\varepsilon)$ for $0 < \varepsilon < 1$. With the default eps = 1e-12, the upper bound is approximately 27.63.

• score.percentile: the percentile of score relative to a row-level training reference rebuilt from the saved target-wise training OOD scores for the requested targets, weights, and row aggregate. For legacy fitted objects that do not contain those saved training scores, the original saved row-level reference is used when possible; otherwise NA.

• targets.used: the number of weighted targets that contributed to each row-level score.

• target.score: optional matrix of target-wise calibrated OOD scores, returned when return.details = TRUE.

• target.delta: optional matrix of raw target-wise reconstruction discrepancies, returned when return.details = TRUE.

• target.reconstruction: when return.reconstruction = TRUE, a data frame containing the saved learners' predictions for the scored targets.

• reconstructed.data: when return.reconstruction = TRUE, the harmonized raw table with scored targets replaced by their reconstructions. Other columns retain their harmonized values. Integer restoration applies to generated reconstructions, preserving observed numeric values in the other columns.

• completed.data: when both return.details and return.reconstruction are TRUE, the raw predictor table after prediction-time imputation. This is distinct from the target reconstruction table.

• info: a list of diagnostics including harmonization details, dropped supervised response columns when present, row-wise unseen-factor flags, learner-loading information, supervised-auxiliary diagnostics when present, the active row aggregate and its arguments, whether the saved row-level calibration was used, and any target-specific issues.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Stekhoven D.J. and Buhlmann P. (2012). MissForest–non-parametric missing value imputation for mixed-type data. Bioinformatics, 28(1):112–118.

Tang F. and Ishwaran H. (2017). Random forest missing data algorithms. Statistical Analysis and Data Mining, 10:363–377.

See Also

impute.rfsrc, rfsrc, predict.rfsrc.

Examples

## ------------------------------------------------------------
## small data example: uses missForest for impute engine
## ------------------------------------------------------------

set.seed(101)
aq <- airquality[, c("Ozone", "Solar.R", "Wind", "Temp", "Month")]
aq$Month <- factor(aq$Month)

id <- sample(1:nrow(aq), 100)
train <- aq[id, ]
test <- aq[-id, ]

## training the imputer
fit <- impute.learn(
  data = train,
  ntree = 25,
  mf.q = 1,
  max.iter = 5,
  full.sweep.options = list(ntree = 25, nsplit = 5)
)

## test time imputation
test.imp <- predict(fit, test, max.predict.iter = 2, verbose = FALSE)
print(head(test.imp))


## ------------------------------------------------------------
## Shared training predictors for a deployment bank
## ------------------------------------------------------------
## Keep these objects separate from the airquality examples below.
local({
  if (requireNamespace("fst", quietly = TRUE) &&
      requireNamespace("data.table", quietly = TRUE)) {
    set.seed(17)
    train <- data.frame(x1 = rnorm(100), x2 = runif(100),
                        group = factor(rep(c("a", "b"), 50)))
    bank.dir <- tempfile("shared-imputer-")
    fit <- impute.learn(
      data = train, target.mode = "all", anonymous = FALSE,
      share.xvar = TRUE, compress = 100, out.dir = bank.dir,
      full.sweep.options = list(ntree = 50, nsplit = 5,
                                terminal.qualts = FALSE, terminal.quants = FALSE),
      verbose = FALSE
    )
    imp <- load.impute.learn(bank.dir, verbose = FALSE)
    test <- train[1:5, ]
    test$x1[c(1, 3)] <- NA_real_
    print(predict(imp, test, verbose = FALSE))
    unlink(bank.dir, recursive = TRUE)
  }
})



## OOD scoring is most informative when every deployment-time
## variable can be reconstructed, so target.mode = "all" is recommended.
## Optional named OOD weights can also be supplied here. Any omitted
## targets receive weight 0, and the saved weights are reused
## automatically later by impute.ood().

ood.fit <- impute.learn(
  data = train,
  ntree = 25,
  mf.q = 1,
  max.iter = 5,
  target.mode = "all",
  save.ood = TRUE,
  full.sweep.options = list(ntree = 25, nsplit = 5),
  verbose = FALSE
)

ood <- impute.ood(ood.fit, test, return.details = TRUE, verbose = FALSE)
print(head(ood$score))
print(head(ood$score.percentile))

## try a more spike-sensitive row aggregate
ood.lp <- impute.ood(ood.fit, test,
                     aggregate = "weighted.lp",
                     aggregate.args = list(p = 4),
                     verbose = FALSE)
print(head(ood.lp$score.percentile))

## ------------------------------------------------------------
## supervised OOD example: regression benchmark
## the user supplies only the raw x-columns at test time
## ------------------------------------------------------------

friedman1_sim <- function(n = 150, p = 10, sigma = 1) {
  X <- matrix(runif(n * p), nrow = n)
  y <- 10 * sin(pi * X[, 1] * X[, 2]) +
       20 * (X[, 3] - 0.5)^2 +
       10 * X[, 4] + 5 * X[, 5] +
       rnorm(n, sd = sigma)
  list(X = X, y = y)
}


trn <- data.frame(friedman1_sim())
tst <- data.frame(friedman1_sim())
xvars <- setdiff(names(trn), "y")


## impute data using missForest, construct a supervised forest
## - supervised forests are used to create auxiliary variables
## - improves test time OOD in supervised problems

sup.fit <- impute.learn(
  data = trn,
  mf.q = 1,
  supervised.formula = y ~ .,
  supervised.args = list(ntree = 50, nsplit = 5),
  full.sweep.options = list(ntree = 25, nsplit = 5),
  save.ood = TRUE,
  verbose = FALSE
)

## add some missing values to the test data
xnew <- tst[, xvars, drop = FALSE]
xnew[sample(seq_len(nrow(xnew)), 5), xvars[1]] <- NA
xnew[sample(seq_len(nrow(xnew)), 5), xvars[2]] <- NA

## imputation
xnew.imp <- predict(sup.fit, xnew, max.predict.iter = 2, verbose = FALSE)
print(head(xnew.imp))

## OOD score
ood.sup <- impute.ood(sup.fit, xnew, verbose = FALSE)
print(head(ood.sup$score.percentile))

## ------------------------------------------------------------
## Save the learned imputer to disk and load it later.
## This explicit save example writes learners kept in memory.
## Uses missForest for the impute engine.
## ------------------------------------------------------------

bundle.dir <- file.path(tempdir(), "aq.imputer")

fit <- impute.learn(
  data = train,
  ntree = 25,
  mf.q = 1,
  max.iter = 5,
  full.sweep.options = list(ntree = 25, nsplit = 5),
  keep.models = TRUE,
  verbose = FALSE
)

## Stronger compression for a deployment bundle.
save.impute.learn(fit, bundle.dir, compress = 100, verbose = FALSE)
imp <- load.impute.learn(bundle.dir, lazy = TRUE, verbose = FALSE)
test.imp <- predict(imp, test, max.predict.iter = 2, verbose = FALSE)

unlink(bundle.dir, recursive = TRUE)



## ------------------------------------------------------------
## Challenging example with factors, uses save/reload
## ------------------------------------------------------------

## load pbc, convert everything to factors
data(pbc, package = "randomForestSRC")
dta <- data.frame(lapply(pbc, factor))
dta$days <- pbc$days
dta$status <- dta$status

## split the data into unbalanced train/test data (25/75)
## the train/test data have the same levels, but different labels
idx <- sample(1:nrow(dta), round(nrow(dta) * .25))
train <- dta[idx,]
test <- dta[-idx,]

## even harder ... factor level not previously encountered in training
levels(test$stage) <- c(levels(test$stage), "fake")
test$stage[sample(seq_len(nrow(test)), 10)] <- "fake"

## train forest
fit <- suppressWarnings(
  impute.learn(Surv(days, status) ~ ., train,
               target.mode = "all",
               save.ood = TRUE,
               keep.models = TRUE)
)

## save/reload
bundle.dir <- file.path(tempdir(), "pbc.imputer")
save.impute.learn(fit, bundle.dir, verbose = FALSE)
imp <- load.impute.learn(bundle.dir, lazy = TRUE, verbose = FALSE)
test.imp <- predict(imp, test, max.predict.iter = 2, verbose = FALSE)
ood <- impute.ood(imp, test, return.details = TRUE, verbose = FALSE)
print(which(ood$info$unseen.rows))
print(summary(test.imp))
unlink(bundle.dir, recursive = TRUE)

---

Random Forest Imputation

Description

Imputes missing values using random forests. The default is on-the-fly imputation (OTFI), with supervised or unsupervised splitting, including an optional autoencoder. Alternatively, missForest predicts missing values one variable or group of variables at a time. Returns an imputed data frame.

Usage

## S3 method for class 'rfsrc'
impute(formula, data,
  ntree = 100, nodesize = 1, nsplit = 10,
  nimpute = 2, fast = FALSE, blocks,
  mf.q, max.iter = 10, eps = 0.01,
  ytry = NULL, always.use = NULL, verbose = TRUE,
  full.sweep = FALSE, restore.integer = TRUE,
  ...)

Arguments

formula	Model formula specifying the responses and predictors for OTFI. Omit it for unsupervised splitting, when no outcome is designated. Ignored for missForest.
data	Data to be imputed. Matrices, tibbles, and data.table objects are converted to a plain data frame. Column names must be unique and nonempty, each column must be a vector, and numeric values must be finite or missing. Bare vectors and nested matrix or list columns are not supported.
ntree	Number of trees grown for each imputation forest.
nodesize	Minimum terminal node size in each tree.
nsplit	Non-negative integer for specifying random splitting.
nimpute	Number of OTFI iterations. Ignored for missForest, whose iterations are controlled by max.iter and eps.
fast	If TRUE, uses fast forest settings to reduce computation. May reduce imputation accuracy.
blocks	Number of row blocks for the main imputation fits. May improve speed for large data, but can reduce imputation accuracy. If omitted, fits use all rows. See Details.
mf.q	Selects missForest when supplied. With mf.q = 1 and always.use = NULL, variables are imputed one at a time. Other positive values select multivariate missForest and specify an approximate fraction (below 1) or number (above 1) of variables with missing values to include in each response group. Takes precedence over splitrule = "auto". See Details.
max.iter	Maximum number of missForest iterations. Must be a positive integer.
eps	Minimum decrease in the change between successive missForest imputations required to continue iterating. Must be finite and nonnegative. See Details.
ytry	Number of response variables used for unsupervised OTFI. With the default splitting method, pseudo-responses are selected for each candidate splitting variable; the default is min(p - 1, max(25, ceiling(sqrt(p)))). With splitrule = "auto", response copies are selected at each node; the default is p, and a smaller positive integer can be supplied. Here p is the number of variables remaining after all-missing columns are removed. Used when formula and mf.q are omitted.
always.use	Names of variables to include in every response group for multivariate missForest. With mf.q supplied, any non-NULL value selects the multivariate version, even if mf.q = 1. Empty vectors and names not present in the data add no responses but still select this version. Ignored for OTFI.
verbose	If TRUE, prints progress during missForest imputation and the optional final sweep.
full.sweep	If TRUE, performs one additional pass after either imputation method. Each variable with missing values in the original data is predicted from the current imputed values of the other variables. Only originally missing entries are replaced. See Details.
restore.integer	If TRUE, rounds imputed values to integers for columns stored as R integers or whose observed numeric values are all integer-valued. Observed values are preserved. Numeric columns retain numeric storage. Original integer columns are returned as integer vectors only when all nonmissing completed values are integral and within R's integer range; otherwise they retain numeric storage. Set to FALSE to keep unrounded imputations. See Details for the integer-detection tolerance.
...	Additional forest controls for the main imputation: mtry, splitrule, bootstrap, sampsize, and samptype. Use splitrule = "auto" for autoencoder splitting when formula and mf.q are omitted. Supplying mf.q selects missForest and ignores the "auto" setting. Controls for the optional final sweep are supplied separately in the named list full.sweep.options. This list accepts the same controls, plus ntree (default 500), nodesize (default NULL), and nsplit (default 10).

Details

Two approaches are available. The default, on-the-fly imputation (OTFI), handles missing values during tree construction and obtains final imputations from terminal nodes. Supplying mf.q selects missForest, which repeatedly uses the other variables to impute each variable or group of variables with missing values. Both approaches support numeric and categorical data. The function returns the imputed data frame; fitted forests and performance summaries are not returned.

On-the-fly imputation (OTFI)

OTFI is the approach underlying na.action = "na.impute" in rfsrc and predict.rfsrc. Tang and Ishwaran (2017, Sections 2.3–2.4) describe the method for both the supervised and unsupervised case.

Split statistics are calculated from the observed values of the variables involved. If an observation is missing the value of the chosen splitting variable, a value is drawn from that variable's observed inbag values in the node to determine its daughter node. This draw is temporary: it is discarded after assignment, and the entry is treated as missing again at the next node. Final replacements are obtained separately from observed values in the terminal nodes and combined across trees.

With a response formula, OTFI uses supervised splitting guided by the specified outcome. For example, Ozone ~ . uses Ozone to guide the splits while imputing missing responses and predictors.

With no formula, OTFI uses unsupervised splitting to learn relationships among the variables. At each node, each of mtry candidate splitting variables is evaluated against ytry randomly selected pseudo-responses using a multivariate composite splitting rule. The best split is selected.

Additionally, for unsupervised data, setting splitrule = "auto" implements a autoencoder type method which regresses variables against themselves. Temporary copies serve as multivariate responses, and the original variables serve as predictors. Missing entries remain missing in both copies for OTFI to process. All response variables are used by default (ytry = p). The temporary responses are discarded after imputation; only the completed original variables are returned or used in the optional final sweep. Other splitrule settings are ignored when the formula is omitted.

nimpute controls the number of OTFI iterations. With a formula, splitrule = "random" can be used for completely random splitting.

missForest

missForest is an unsupervised imputation procedure in the sense that no outcome needs to be designated: every variable can be imputed using the others. Its individual forests use supervised regression or classification. The formula argument is ignored; an outcome column included in data is treated like any other variable.

The implementation first fills missing values using a small forest with completely random splitting. These initial imputations supply values for the predictors in subsequent fits.

With mf.q = 1 and always.use = NULL, the method follows the variable-by-variable approach of Stekhoven and Buehlmann (2012). For each variable with missing values, a forest is fit using the rows where that variable was observed. Its missing entries are replaced by predictions from this forest. Later fits use the imputations already updated for earlier variables, and the process repeats.

The multivariate extension groups variables with missing values and imputes each group using the remaining variables (Tang and Ishwaran, 2017, Section 2.5). Group sizes are approximate and controlled by mf.q. Within each group, response entries that were originally missing are set back to missing, while the current imputations remain available in the predictors. A forest with multivariate composite splitting imputes the responses. Groups are formed anew at each iteration. Variables named in always.use are added to every response group.

For both versions, iteration is controlled by max.iter and eps. After each pass, changes at originally missing entries are summarized by factor disagreement and scaled numeric changes, averaged over the variables. Iteration stops when this measure fails to decrease by at least eps, or when max.iter is reached. The last completed pass is retained.

Data handling and optional refinements

Rows and columns with no observed values are removed before imputation. If no missing values remain, the data are returned without fitting a forest. If the entire data frame is missing, it is returned without imputation.

For large data sets, blocks restricts the main imputation fits to separate groups of rows. This can reduce computation but limits the information available to each fit. Autoencoder fits use only response variables with observed values in that block and limit ytry accordingly. The missForest initialization and optional final sweep use all retained rows.

Setting full.sweep = TRUE adds a final pass after either method. Each variable that originally had missing values is fit on its observed rows using the current imputed predictors. Only its originally missing entries are replaced, and later fits use earlier updates from the sweep. Sweep controls are supplied in full.sweep.options; see ... above.

With restore.integer = TRUE, integer support is determined from the original data before imputation. Numeric values are compared with the nearest integer using a fixed absolute floating-point tolerance. Rounding applies only to imputed entries: after the missForest initialization, before and after the optional final sweep, and before the data are returned. Values are not rounded between the individual missForest updates.

Value

A data frame with missing entries replaced by imputed values and observed values preserved. All-missing rows and columns are removed as described in Details.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Stekhoven D.J. and Buehlmann P. (2012). MissForest–non-parametric missing value imputation for mixed-type data. Bioinformatics, 28(1):112–118. doi:10.1093/bioinformatics/btr597.

Tang F. and Ishwaran H. (2017). Random forest missing data algorithms. Statistical Analysis and Data Mining, 10(6):363–377. doi:10.1002/sam.11348.

See Also

impute.learn.rfsrc, rfsrc, rfsrc.fast, predict.rfsrc

Examples

## ------------------------------------------------------------
## example of survival imputation
## ------------------------------------------------------------

## default OTFI: unsupervised splitting
data(pbc, package = "randomForestSRC")
pbc1.d <- impute(data = pbc)

## autoencoder OTFI: regress the variables against themselves
pbc.auto <- impute(data = pbc, splitrule = "auto")

## optionally use fewer response variables at each node
pbc.auto.sub <- impute(data = pbc, splitrule = "auto", ytry = 3)

## OTFI with splits guided by the survival outcome
f <- as.formula(Surv(days, status) ~ .)
pbc2.d <- impute(f, data = pbc, nsplit = 3)

## random splitting can be reasonably good
pbc3.d <- impute(f, data = pbc, splitrule = "random", nimpute = 5)

## optional final sweep after OTFI
pbc3.fs <- impute(f, data = pbc, splitrule = "random", nimpute = 5,
                  full.sweep = TRUE)

## ------------------------------------------------------------
## example of regression imputation
## ------------------------------------------------------------

air1.d <- impute(data = airquality, nimpute = 5)
air2.d <- impute(Ozone ~ ., data = airquality, nimpute = 5)
air3.d <- impute(Ozone ~ ., data = airquality, fast = TRUE)

## final sweep with custom options (e.g., larger forest)
air3.fs <- impute(Ozone ~ ., data = airquality, nimpute = 5,
                  full.sweep = TRUE,
                  full.sweep.options = list(ntree = 1000, nodesize = 5, nsplit = 0,
                                            mtry = 3, splitrule = "random"))

## ------------------------------------------------------------
## multivariate missForest imputation
## ------------------------------------------------------------

data(pbc, package = "randomForestSRC")

## missForest algorithm - uses 1 variable at a time for the response
pbc.d <- impute(data = pbc, mf.q = 1)

## multivariate missForest - use about 10 percent of variables
## with missing values as responses in each group
pbc.mv <- impute(data = pbc, mf.q = .10)

## missForest but faster by using random splitting
pbc.fast <- impute(data = pbc, mf.q = 1, splitrule = "random")

## missForest + final sweep
pbc.fast.fs <- impute(data = pbc, mf.q = 1, splitrule = "random",
                      full.sweep = TRUE)

---

Acquire Maximal Subtree Information

Description

Extract maximal subtree information from a RF-SRC object. Used for variable selection and identifying interactions between variables.

Usage

## S3 method for class 'rfsrc'
max.subtree(object,
  max.order = 2, sub.order = FALSE, conservative = FALSE, ...)

Arguments

object	An object of class (rfsrc, grow) or (rfsrc, forest).
max.order	Non-negative integer specifying the maximum interaction order for which minimal depth is calculated. Defaults to 2. Set max.order=0 to return first-order depths only. When max.order=0, conservative is automatically set to FALSE.
sub.order	Logical. If TRUE, returns the minimal depth of each variable conditional on every other variable. Useful for investigating variable interdependence. See Details.
conservative	Logical. If TRUE, uses a conservative threshold for selecting variables based on the marginal minimal depth distribution (Ishwaran et al., 2010). If FALSE, uses the tree-averaged distribution, which is less conservative and typically identifies more variables in high-dimensional settings.
...	Additional arguments passed to or from other methods.

Details

The maximal subtree for a variable x is the largest subtree in which the root node splits on x. The largest possible maximal subtree is the full tree (root node), though multiple maximal subtrees may exist for a variable. A variable may also have no maximal subtree if it is never used for splitting. See Ishwaran et al. (2010, 2011) for further discussion.

The minimal depth of a maximal subtree-called the first-order depth-quantifies the predictive strength of a variable. It is defined as the distance from the root node to the parent of the closest maximal subtree for x. Smaller values indicate stronger predictive impact. A variable is flagged as strong if its minimal depth is below the mean of the minimal depth distribution.

The second-order depth is the distance from the root to the second-closest maximal subtree of x. To request depths beyond first order, use the max.order option (e.g., max.order = 2 returns both first and second-order depths). Set max.order = 0 to retrieve first-order depths for each variable in each tree.

Set sub.order = TRUE to obtain the relative minimal depth of each variable j within the maximal subtree of another variable i. This returns a p x p matrix (with p the number of variables) whose entry (i,j) is the normalized relative depth of j in i's subtree. Entry (i,i) gives the depth of i relative to the root. Read the matrix across rows to assess inter-variable relationships: small (i,j) entries suggest interactions between variables i and j.

For competing risks, all analyses are unconditional (non-event specific).

Value

Invisibly returns a list with the following components:

order	Matrix of order depths for each variable up to max.order, averaged over trees. The matrix has p rows and max.order columns, where p is the number of variables. If max.order = 0, returns a matrix of dimension p x ntree containing first-order depths for each variable by tree.
count	Average number of maximal subtrees per variable, normalized by tree size.
nodes.at.depth	List of vectors recording the number of non-terminal nodes at each depth level for each tree.
sub.order	Matrix of average minimal depths of each variable relative to others (i.e., conditional minimal depth matrix). NULL if sub.order = FALSE.
threshold	Threshold value for selecting strong variables based on the mean of the minimal depth distribution.
threshold.1se	Conservative threshold equal to the mean minimal depth plus one standard error.
topvars	Character vector of selected variable names using the threshold criterion.
topvars.1se	Character vector of selected variable names using the threshold.1se criterion.
percentile	Percentile value of minimal depth for each variable.
density	Estimated density of the minimal depth distribution.
second.order.threshold	Threshold used for selecting strong second-order depth variables.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

Ishwaran H., Kogalur U.B., Chen X. and Minn A.J. (2011). Random survival forests for high-dimensional data. Statist. Anal. Data Mining, 4:115-132.

See Also

holdout.vimp.rfsrc, vimp.rfsrc

Examples

## ------------------------------------------------------------
## survival analysis
## first and second order depths for all variables
## ------------------------------------------------------------

data(veteran, package = "randomForestSRC")
v.obj <- rfsrc(Surv(time, status) ~ . , data = veteran)
v.max <- max.subtree(v.obj)

# first and second order depths
print(round(v.max$order, 3))

# the minimal depth is the first order depth
print(round(v.max$order[, 1], 3))

# strong variables have minimal depth less than or equal
# to the following threshold
print(v.max$threshold)

# this corresponds to the set of variables
print(v.max$topvars)

## ------------------------------------------------------------
## regression analysis
## try different levels of conservativeness
## ------------------------------------------------------------

mtcars.obj <- rfsrc(mpg ~ ., data = mtcars)
max.subtree(mtcars.obj)$topvars
max.subtree(mtcars.obj, conservative = TRUE)$topvars

---