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

## ------------------------------------------------------------
## Standard analysis
## ------------------------------------------------------------

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

Classification Performance Metrics

Description

Evaluate classification predictions using confusion matrices, class-specific misclassification rates, ROC area under the curve (AUC), Brier score, log loss, and binary precision-recall summaries. The helpers take observed responses and predicted classes or probabilities directly. A class-assignment helper is also provided.

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

Details

Supplying predictions

These helpers score the predictions supplied to them; they do not fit a forest or automatically choose OOB predictions. For a grow object, use o$yvar with o$predicted.oob for OOB probability scores, or o$class.oob for its stored class assignments. For new-data evaluation, use the prediction object's responses and predicted or class components.

Keep responses and predictions aligned by row. For a common scoring sample, restrict to observed responses and finite probabilities before calling the helpers. In particular, get.logloss() and get.misclass.error() require nonmissing observed labels. The examples illustrate this preparation for OOB predictions.

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() assigns each row to a class with the largest probability, sampling among tied maximizers. A row whose probabilities are all missing receives NA. With two class relative frequencies supplied, it assigns the minority class when that class's probability is at least its supplied relative frequency. Otherwise it assigns the majority class. If the frequencies tie, the first probability column is treated as the minority. Supply finite probabilities for this frequency-based rule.

get.confusion() has observed classes in rows and predicted classes in columns. Its final class.error column is one minus the fraction correctly classified within each observed class, rounded to four decimal places. Missing response/prediction pairs do not enter the counts. When supplied a probability matrix, it uses the largest-probability rule. To evaluate RFQ or another thresholded rule, supply that rule's predicted factor instead, for example o$class.oob.

get.misclass.error() returns one misclassification rate per observed class, in sort(unique(y)) order. It returns class-specific rates rather than an overall error rate. Missing predicted labels make the corresponding class rate unavailable.

ROC area under the curve

get.auc() uses the pairwise multiclass construction of Hand and Till. For each pair of classes, it calculates a rank-based AUC using each class's probability column as its score and averages the two directions. It then averages over available class pairs. For complementary two-class probabilities, the two directions agree. Larger values indicate better discrimination.

Missing observed labels and nonfinite scores are excluded from the relevant rank calculations. Each direction requires at least two finite scores in each of the two classes. Unavailable pairs are omitted; the result is NA when no pair is available.

Brier score

For J classes, let

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

The scalar returned by get.brier.error() is

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

with normalized = TRUE, and

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

with 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 usual mean squared error of either class probability, and the normalized score is four times that value.

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

Log loss

With every response level represented, get.logloss() averages -\log(p_{i,y_i}), using the natural logarithm. Smaller values are better. With robust = TRUE, infinite losses from zero probabilities are excluded, rather than replaced using a probability floor. Use robust = FALSE to retain these infinite losses. Missing losses are excluded in either case. An unused factor level contributes one zero term to the helper's average.

Precision-recall summaries

These helpers apply to two classes. For truth coded 0/1, including a factor with those two labels, 1 is positive. For other two-level factors, the less frequent observed class is positive; ties select the first factor level. This choice is made before excluding missing scores.

A score vector or single column always refers to the positive class. With two columns, named columns are matched to the original response levels; unnamed columns follow factor-level order, or 0, 1 order for a nonfactor response. The positive-class column is then selected. Rows with a missing response or nonfinite selected score are excluded.

get.pr.auc() returns the precision-recall area and a random reference area equal to the positive-class proportion among the scored observations. The area uses the analytic precision-recall interpolation implemented by the package. get.pr.curve() supplies recall, precision, and threshold values for plotting; its rows run from high to low recall. Both classes must have a usable score. Larger precision-recall areas are better.

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.

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 Objects

Description

Save and load the forest component of a fitted random forest using a directory of files. The loaded forest can be passed to predict(). Corresponding helpers save and load lists of forests.

Usage

fast.save(o, path = NULL, testing = TRUE, units = "Mb")

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

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

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

Arguments

Details

Saving and loading a forest

fast.save() saves o$forest, including its retained training data and tree information. fast.load() returns the reconstructed forest object, which can be used directly with predict.rfsrc.

Components outside o$forest, such as the grow object's prediction and performance summaries, are not saved. To preserve the complete R object instead, use saveRDS and readRDS.

Saving requires the fst and data.table packages; loading requires fst. Large tabular components are written with fst::write_fst(), and the remaining forest metadata are stored with saveRDS(..., compress = FALSE). The files use a ‘.rda’ suffix despite these different formats. Keep the entire directory together and read it using fast.load().

Directories

Both save helpers remove the destination directory and all its contents before writing. Use a dedicated directory for the forest or forest list, rather than a directory containing other work. The destination is then created, including parent directories as needed.

Saving takes the full destination in path. Loading takes its name in directory and its parent in path. For a destination stored in save.path, the matching load call is fast.load(basename(save.path), path = dirname(save.path)). With the default save destination, use fast.load("forest").

Lists of forests

fast.save.list() saves the elements in subdirectories ‘forest1’, ‘forest2’, and so on, inside path. Original list names are not stored.

fast.load.list() loads every entry returned by list.files() in the list directory and returns an unnamed list in that order. Keep only saved forest subdirectories there. Directory-name ordering can differ from the original numeric index order: for example, ‘forest10’ can precede ‘forest2’. Use fast.load() on a named subdirectory when selecting a particular saved element.

Sampling-size function

The saved sampling-size information is sampfrac = o$forest$sampsize(1). Loading reconstructs the function as function(x) x * sampfrac. This preserves a proportional sampling-size rule; it does not preserve an arbitrary nonlinear custom function. Use saveRDS() when that original function must be retained.

Value

fast.save

Used for its file-writing side effect. Its return value is the memory-usage matrix from the final gc(FALSE) call, not the saved forest.

fast.load

The reconstructed forest component, retaining its saved class, ordinarily c("rfsrc", "forest", family). It is a forest object rather than the original grow object.

fast.save.list

An unnamed list of the return values from the individual fast.save() calls. The forests are saved in the destination subdirectories.

fast.load.list

An unnamed list of reconstructed forest objects in the directory-listing order described above.

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)
}


## ------------------------------------------------------------
## 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)) {
    cat("check coherence for outcome:", yn, "\n")
    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:

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

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:

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

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):

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.

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

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 ordinary 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

G-mean is the geometric mean of minority sensitivity and majority specificity, \sqrt{\mathrm{sensitivity}\, \mathrm{specificity}}. Larger values indicate better balance between the two class-specific rates. Misclassification and Brier scores are error measures for which smaller values are better. Changing perf.type selects the forest performance measure; it does not select the RFQ decision rule.

Request permutation VIMP explicitly with importance = "permute". For comparisons, use a common perf.type, splitting rule, and block.size, unless the comparison is intended to study those settings. Increase ntree as needed to assess the stability of OOB performance and VIMP.

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 accepts a fitted or predicted forest object, or a two-level factor supplied as obj together with prob. For a forest object, it prefers predicted.oob when that component exists and is not entirely missing; otherwise it uses predicted. A partially missing OOB component is not filled with full-ensemble predictions.

A numeric vector or one-column matrix or data frame represents minority-class probabilities, with one value per response. No column name is needed for these inputs. A two-column matrix or data frame contains both class probabilities. Named columns are matched to the response factor levels; unnamed columns follow factor-level order. The response and probability rows must already be aligned.

Finite probabilities are restricted to [0,1] using the existing clipping rule. Finite two-column rows must sum to one within numerical tolerance. Rows with a missing response or a nonfinite probability are excluded from scoring; both columns must be finite for a two-column input. Unavailable scores are returned as NA. Class counts, minority identity, and the default threshold are determined from all observed evaluation responses before excluding unavailable predictions.

The helper determines the minority class and its proportion from the response supplied for evaluation. If the class counts tie, the first factor level is selected. The default threshold = NULL uses this evaluation proportion for all forest methods.

With threshold = "auto", the saved forest$rfq flag selects the cutoff: TRUE uses the same evaluation prevalence as NULL, and otherwise the cutoff is 0.5. This applies to both grow and prediction objects. For RFQ test predictions the prevalence is therefore calculated from the evaluation responses, not substituted from the training responses. Splitting rules and perf.type do not determine whether the classifier is RFQ. Direct response/probability inputs have no fitted-classifier flag, so "auto" uses 0.5. Use threshold = "bayes" to request 0.5 explicitly for any classifier.

All settings threshold the probabilities rather than reading the stored class component. Supply a single finite numeric threshold to examine a different cutoff; finite values outside [0,1] are clipped to that interval. Threshold ties are assigned to the minority class. A supplied threshold changes the cutoff, not how the minority class is identified or how the forest was fitted.

Besides the four rates sens, spec, prec, and npv, the summary includes misclass, brier, brier.norm, auc, logloss, and pr.auc. For complementary binary probabilities, brier is the mean squared error of the minority probability and brier.norm is four times this value. The current log-loss helper excludes infinite losses rather than clipping zero probabilities.

F1 is the harmonic mean of precision and sensitivity; F1mod is the harmonic mean of sensitivity, specificity, precision, and negative predictive value. F1gmean averages F1 and G-mean; F1modgmean averages F1mod and G-mean. A harmonic mean is zero when all its component rates are defined and at least one is zero. It is NA when any required rate is undefined.

robust = TRUE adds one to the two diagonal counts when calculating the four rates and their derived summaries. It does not change the returned confusion counts, misclassification error, or probability-based scores. Use confusion = TRUE to return the confusion matrix in addition to the metrics. Its 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, and requested performance and importance values. 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, subsample

Examples

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

set.seed(19)
rfq <- imbalanced(f, dta, ntree = 300)
print(rfq)
print(get.imbalanced.performance(rfq, confusion = TRUE))

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

brf <- imbalanced(f, dta, ntree = 300, method = "brf")
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", "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.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)

## ------------------------------------------------------------
## Optional sampling and grow-interface choices
## ------------------------------------------------------------
rfq.sample <- imbalanced(f, dta, ntree = 300, ratio = .25)
print(get.imbalanced.performance(rfq.sample))
rfq.fast <- imbalanced(f, dta, ntree = 300, fast = TRUE)
print(get.imbalanced.performance(rfq.fast))

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)

save.impute.learn.rfsrc(object, path, wipe = TRUE, verbose = TRUE)

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

Details

A predictive imputer is calculated in two stages.

The training data are first processed. Variables that are not real-valued are coerced 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.

If the resulting training data contain missing values, the first stage uses impute as the imputation engine to complete the training data. Options are specified exactly as in impute. In particular, mf.q = 1 with always.use = NULL updates one target at a time, other positive settings give the multivariate missForest generalization, and if mf.q is omitted, on the fly imputation is used if formula is specified, otherwise default unsupervised imputation is used. If the training data are already complete and target.mode = "all", this initial imputation step is skipped.

In the second stage, a full forest sweep is fit on the completed training data. For each target selected by target.mode, rows where that target was observed are used to fit a forest with that target on the left-hand side and the predictors selected by deployment.xvars on the right-hand side. These saved forests all use the same completed training table; fitting the learner bank does not further update that table.

By default, deployment.xvars = NULL allows every non-target column to be used as a predictor. This is convenient, but it can also introduce leakage if the training data include outcomes, future-only variables, identifiers, or any fields that will not be available when the learned imputer is applied to new data. Restrict deployment.xvars when that is a concern.

If supervised.formula is supplied, the function also fits an internal supervised forest from the training data. Out-of-bag predicted values for the training sample, and predicted values for new data, are appended internally as auxiliary predictors. These auxiliary variables are not used as imputation targets, however, they are available as predictors for the saved forests and can therefore influence both test-time imputation and OOD scoring. Users who want only the basic unsupervised learned imputer should leave supervised.formula unspecified.

When supervised.formula is used, deployment.xvars continues to restrict only the raw predictors. The internally created auxiliary variables are added automatically to the predictor sets for the saved raw targets. This lets the learned imputer benefit from supervised signal without asking the user to prepare and pass those auxiliary columns manually.

When the imputer is saved to disk, each full-sweep learner is written separately using fast.save. Loading uses fast.load. In practice this gives a small manifest plus a directory of saved learners. The fst package is therefore required for save and load operations. The explicit save method can write learners either from memory or by reloading them from an attached saved path. If supervised.formula is used, the internal supervised forest is saved alongside the target-wise learner bank. Both training-time and explicit saving write and read back 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.

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

Prediction starts by matching newdata to the training schema, filling missing values with training means or modes, and then applying one or more full-sweep passes. Only the targets selected by target.mode are updated by saved learners. In supervised mode, this same test-time preparation first initializes the raw predictor table, then computes the internal auxiliary variables from the saved supervised forest, and then runs the forest test-time sweep. The auxiliary predictors remain fixed throughout these passes. Supervised response columns, if present in newdata but not in the learned raw predictor block, are dropped before this step and play no role at prediction or scoring time. A target update requires one valid prediction per requested row. Failed or unavailable predictions retain the preceding imputed values and are recorded in target.issues. A pass with no valid model updates is reported separately from convergence.

If target.mode = "missing.only", a variable that was complete in training but missing in new data is initialized from the training fit but does not receive a model-based update. Use target.mode = "all" if missing values may appear later in any raw variable. Complete training data also require target.mode = "all", because otherwise there are no missing variables from which to determine the saved targets.

If save.ood = TRUE, the fit also stores an OOD reference in the manifest. For each saved target, the out-of-bag prediction from the fitted learner is compared with the observed training value to form a target-wise reconstruction discrepancy. Continuous and integer targets use absolute reconstruction error. Factor targets use negative log predictive probabilities. Unavailable or invalid probabilities remain missing discrepancies; a genuine zero probability is scored using the probability floor. Each learner entry records n.oob.finite and n.oob.nonfinite when OOD reference construction is enabled.

The row-level OOD calibration stored at fit time is built by aggregating the target-wise training scores with a weighted mean using weight. If weight is omitted at fit time, all saved OOD targets receive weight 1. If a named vector is supplied, entries are matched by target name, omitted saved OOD targets receive weight 0, and the resulting weighting scheme is carried forward for later deployment-time scoring.

impute.ood first completes the predictor side of newdata using the same harmonization, initialization, and iterative sweep logic used by predict.impute.learn. It then reconstructs each requested raw target directly from its saved learner and compares the reconstruction with the observed raw target value supplied in newdata. If a scored raw target is missing in a row, that target does not contribute to that row's OOD score. Raw target discrepancies are converted to target-wise OOD scores using the saved target-specific training references. Discrepancies strictly above a nonempty reference's maximum, including positive infinity, receive its largest stored probability. The existing convention for ties is retained, including equality at the largest quantile. Missing discrepancies and empty references remain unscored. In supervised mode, auxiliary variables participate as predictors in the reconstruction, but the supervised response itself is not used at deployment time.

The row-level OOD score combines those calibrated target-wise scores over the targets that are both observed and scoreable for that row. By default, impute.ood uses a bounded product rule, but the row aggregate can be changed to weighted mean, a weighted L_p rule, a log-tail weighted L_p rule, or a top-k rule. This makes it possible to explore row scores that are more sensitive to sparse but severe coordinate shifts. By default, impute.ood reuses the same OOD weights saved during impute.learn, so a pipeline can fix its weighting scheme once upstream and carry it forward automatically.

A second component, score.percentile, is obtained by rebuilding the row-level training reference from the saved target-wise training OOD scores using the requested target subset, the active weight vector, and the active row aggregate. This means percentile calibration remains available when the user leaves the saved weights in place, overrides them at test time, scores only a subset of the saved OOD targets, or experiments with alternate row aggregates.

Unseen factor levels are tracked row-wise during harmonization. Because such values are immediate anomalies relative to the training schema, impute.ood flags those rows and assigns them the maximum row-level score. If the unseen level occurs in a scored target itself, the corresponding target-level discrepancy is also treated as 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.

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.

• 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 before schema restoration.

• 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))



## 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
)

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)

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)

Impute Only Mode

Description

Fast imputation mode. A random forest is grown and used to impute missing data. No ensemble estimates or error rates are calculated. Optionally, a final sweep can be performed to re-fit each variable that had original missingness on the final covariates and overwrite only its originally-missing entries.

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

Details

1. For a table containing both observed and missing values, observations and variables with all values missing are removed before imputation. A complete table or a table whose entries are all missing is returned without imputation after container normalization.

2. A forest is grown and used solely for imputation. No ensemble statistics (e.g., error rates) are computed. Use this function when imputation is the only goal.

3. For standard imputation (not missForest), splits are based only on non-missing data. If a split variable has missing values, they are temporarily imputed by randomly drawing from in-bag, non-missing values to allow node assignment.

4. If mf.q is specified, multivariate missForest imputation is applied (Stekhoven and B\"uhlmann, 2012). A fraction (or integer count) of variables are selected as multivariate responses, predicted using the remaining variables with multivariate composite splitting. Each round imputes a disjoint set of variables, and the full cycle is repeated until convergence, controlled by max.iter and eps. Setting mf.q = 1 reverts to response-wise missForest when always.use = NULL. A non-NULL always.use, including an empty or entirely unmatched vector, selects the multivariate branch. The change statistic averages factor disagreement and scaled numeric changes over originally missing entries. Iteration stops when this statistic no longer decreases by at least eps, or when max.iter is reached. The last completed pass is retained.

5. If no formula is provided and mf.q is omitted, unsupervised splitting is used. The default ytry is sqrt(p), where p is the number of retained variables. For each of mtry candidate variables, a random subset of ytry variables is selected as pseudo-responses. A multivariate composite splitting rule is applied, and the split is made on the variable yielding the best result (Tang and Ishwaran, 2017).

6. If no missing values remain after preprocessing, the function returns the processed data without further action.

7. Additional forest controls forwarded through ... are mtry, splitrule, bootstrap, sampsize, and samptype. The optional final sweep has its own full.sweep.options list.

8. Optional final sweep: if full.sweep = TRUE, a post-imputation sweep is performed for every variable with original missingness. Each such variable is re-fit on its observed rows using the final imputed covariates, and predictions overwrite only the originally missing entries. Defaults for the sweep are ntree = 500, nodesize = NULL, nsplit = 10, and can be customized via full.sweep.options passed through .... This applies to both standard and missForest modes.

9. If restore.integer = TRUE, integer support is detected from the original retained data before imputation begins. Numeric values are compared with the nearest integer using an absolute floating-point tolerance that does not increase with their magnitude. Restoration is applied after the missForest rough-fix initialization, before and after the optional final sweep, and once more before returning the completed data. The missForest iterative updates themselves are not rounded between variable or block updates.

Value

Invisibly, the data frame containing the original data with imputed data overlaid.

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.learn.rfsrc, rfsrc, rfsrc.fast

Examples

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

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

## imputation using outcome splitting
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 (standard imputation)
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 10 percent of variables as responses
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

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:

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

Extracting Multivariate Values

Description

Extract predictions, performance errors, variable importance (VIMP), and case-specific values from fitted random forests. The helpers collect response-specific results from multivariate regression, multivariate classification, and mixed-outcome forests. They also provide a common interface for univariate results. A helper for constructing a multivariate formula is included.

Usage

get.mv.predicted(obj, oob = TRUE)

get.mv.error(obj, standardize = FALSE, pretty = TRUE, block = FALSE)

get.mv.error.block(obj, standardize = FALSE)

get.mv.vimp(obj, standardize = FALSE, pretty = TRUE)

get.mv.cserror(obj, standardize = FALSE)

get.mv.csvimp(obj, standardize = FALSE)

get.mv.formula(ynames)

Arguments

Details

Predictions

get.mv.predicted() combines the stored predictions into a matrix with observations in rows. Each regression response supplies one column. Each classification response supplies one probability column per class, named response.class. The responses follow obj$yvar.names; class columns retain their stored order.

The default uses OOB predictions when available. An existing OOB component is retained even if some or all of its entries are missing. The fallback to predicted occurs only when that entire component is NULL. Use oob = FALSE to explicitly select full-ensemble or new-data predictions.

For right-censored survival, the helper extracts the stored mortality prediction. For competing risks, it extracts the event-specific predicted values, with names of the form response.event. It does not extract the time-indexed survival, cumulative hazard, or cumulative incidence arrays.

Performance errors and variable importance

get.mv.error() extracts the final stored performance error for each response. The error measure is the one used when the object was fitted or predicted; the helper does not recalculate it or change perf.type. With pretty = FALSE, classification output includes all and the class-specific errors. Survival output retains its stored error columns, including event-specific columns for competing risks.

get.mv.error.block() returns the full stored block-error sequence for each response. It is equivalent to get.mv.error(obj, standardize = standardize, block = TRUE). It does not choose a block size or form new blocks.

get.mv.vimp() collects the stored importance components. With pretty = TRUE, predictors index rows and results are combined across responses in columns. With pretty = FALSE, each response has its own matrix, preserving classification all and class-specific columns and competing-risk event columns. Request importance during fitting or use vimp before extracting it.

To select a response or predictor, subset the returned vector, matrix, or response-named list. These extraction calls do not fit trees, run prediction, or average errors or VIMP across responses.

Case-specific values

get.mv.cserror() obtains the case-specific error as cse.num / cse.den from the stored response component. get.mv.csvimp() similarly obtains case-specific VIMP as csv.num / csv.den. These quantities must have been retained during fitting or prediction; the helpers do not generate them from the ensemble predictions. In particular, case-specific error is not obtained by applying a loss to predicted.oob and the observed response.

With one response, the case-specific values are returned directly; with multiple responses, they are returned in a response-named list. Case-specific VIMP has observations in rows and the stored VIMP variables in columns. Variable names are attached when available from the importance component. Both helpers return NULL for right-censored survival and competing risks.

Standardization

For a regression response Y, standardize = TRUE divides each extracted error or importance value by var(Y, na.rm = TRUE). This uses the response values on the supplied object: training responses for a grow result and evaluation responses for a test-prediction result. It leaves predictions and all classification and survival values unchanged.

get.mv.error() and get.mv.vimp() divide by this variance directly; a zero or unavailable variance can therefore produce nonfinite values. The two case-specific helpers instead use divisor one when that variance is zero or NA. No standardization is applied by default.

Availability of stored results

Supply an object containing the requested output for its responses. Missing predictions and numerical NA values remain as stored. get.mv.error() returns NULL when all response errors are absent. Otherwise its compact output uses NA for an absent response error, while list output retains a NULL entry.

get.mv.vimp(), get.mv.cserror(), and get.mv.csvimp() return NULL when the first response has no corresponding component, even if a later response has one. When a response-named list is returned, absent later components remain NULL.

Constructing a multivariate formula

get.mv.formula(ynames) returns a formula of the form Multivar(y1, y2, ...) ~ .. The named responses may be continuous, factors, or a mixture; the forest fit determines their types from the supplied data. The dot denotes the remaining data columns. This helper constructs the formula only.

Value

get.mv.predicted

A numeric matrix with one row per observation and columns for the response predictions, class probabilities, or event-specific predictions. A single prediction column remains a matrix.

get.mv.error

A named numeric vector by default, or a response-named list when pretty = FALSE. Without blocking, each entry is the final stored error or final error row. With block = TRUE, list entries contain the stored block-error sequences. Returns NULL if all error components are absent.

get.mv.error.block

A response-named list of block-error vectors or matrices, or NULL when all are absent.

get.mv.vimp

A numeric matrix by default, or a response-named list of matrices when pretty = FALSE. Returns NULL under the availability rule described above.

get.mv.cserror

Case-specific error values for a single response, retaining their stored vector or array dimensions, or a response-named list for multiple responses. Returns NULL for survival families or under the availability rule above.

get.mv.csvimp

A case-by-variable matrix for a single response, or a response-named list of these matrices for multiple responses. Returns NULL for survival families or under the availability rule above.

get.mv.formula

An R formula with the supplied response names on the left of ~ and a dot on the right.

See Also

rfsrc, predict.rfsrc, vimp, subsample, classification.performance

Examples

## ------------------------------------------------------------
## A basic multivariate analysis
## ------------------------------------------------------------
o <- rfsrc(cbind(Ozone, Temp) ~ ., data = na.omit(airquality))
print(head(get.mv.predicted(o)))
print(get.mv.error(o))


## ------------------------------------------------------------
## Select a response from the stored results
## ------------------------------------------------------------
pred.oob <- get.mv.predicted(o)
print(head(pred.oob[, "Temp", drop = FALSE]))
print(get.mv.error(o)["Temp"])
print(get.mv.error(o, standardize = TRUE))
print(head(get.mv.predicted(o, oob = FALSE)))

## ------------------------------------------------------------
## Formula construction, VIMP, and block errors
## ------------------------------------------------------------
f <- get.mv.formula(c("Ozone", "Temp"))
print(f)

o.vimp <- rfsrc(f, data = na.omit(airquality), ntree = 100,
                importance = "permute", block.size = 10)
print(get.mv.vimp(o.vimp))
print(get.mv.vimp(o.vimp, standardize = TRUE))
print(get.mv.vimp(o.vimp, pretty = FALSE)[["Temp"]])
print(head(get.mv.error.block(o.vimp)[["Temp"]]))

## Optional case-specific values are NULL when not saved in the object.
print(get.mv.cserror(o.vimp))
print(get.mv.csvimp(o.vimp))

## ------------------------------------------------------------
## Mixed outcomes: retain class-specific entries
## ------------------------------------------------------------
f.mix <- get.mv.formula(c("Sepal.Length", "Species"))
mix <- rfsrc(f.mix, data = iris, ntree = 100,
             importance = "permute", block.size = 10)
print(colnames(get.mv.predicted(mix)))
print(get.mv.error(mix))
print(get.mv.error(mix, pretty = FALSE)[["Species"]])
print(get.mv.vimp(mix, pretty = FALSE)[["Species"]])

## ------------------------------------------------------------
## Extract predictions and errors for held-out observations
## ------------------------------------------------------------
dta <- na.omit(airquality)
set.seed(17)
train <- sample(seq_len(nrow(dta)), floor(.7 * nrow(dta)))
fit <- rfsrc(f, data = dta[train, ], ntree = 100)
p.test <- predict(fit, newdata = dta[-train, ])
print(head(get.mv.predicted(p.test, oob = FALSE)))
print(get.mv.error(p.test))

Nutrigenomic Study

Description

Investigates the effects of five dietary treatments on 21 liver lipids and 120 hepatic gene expressions in wild-type and PPAR-alpha deficient mice. The examples fit a multivariate regression forest with lipid expression profiles as the responses and gene expression, diet, and genotype as predictors.

Details

The examples use get.mv.formula to construct the response formula, and get.mv.error and get.mv.vimp to extract response-specific performance errors and variable importance (VIMP). For these regression responses, standardize = TRUE divides each error or importance value by the variance of the corresponding response. See Extracting Multivariate Values (multivariate.values) for the helper interfaces and returned formats.

References

Martin P.G. et al. (2007). Novel aspects of PPAR-alpha-mediated regulation of lipid and xenobiotic metabolism revealed through a nutrigenomic study. Hepatology, 45(3), 767–777.

See Also

rfsrc, predict.rfsrc, Extracting Multivariate Values (multivariate.values)

Examples

## ------------------------------------------------------------
## multivariate regression forests using Mahalanobis splitting
## lipids (all real values) used as the multivariate y
## ------------------------------------------------------------

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

## parse into y and x data
ydta <- nutrigenomic$lipids
xdta <- data.frame(nutrigenomic$genes,
                   diet = nutrigenomic$diet,
                   genotype = nutrigenomic$genotype)

## multivariate mixed forest call
obj <- rfsrc(get.mv.formula(colnames(ydta)),
             data.frame(ydta, xdta),
             importance=TRUE, nsplit = 10,
             splitrule = "mahalanobis")
print(obj)

## ------------------------------------------------------------
## plot the standarized performance and VIMP values
## ------------------------------------------------------------

## acquire the error rate for each of the 21-coordinates 
## standardize to allow for comparison across coordinates
serr <- get.mv.error(obj, standardize = TRUE)

## acquire standardized VIMP 
svimp <- get.mv.vimp(obj, standardize = TRUE)

par(mfrow = c(1,2))
plot(serr, xlab = "Lipids", ylab = "Standardized Performance")
matplot(svimp, xlab = "Genes/Diet/Genotype", ylab = "Standardized VIMP")


## ------------------------------------------------------------
## plot some trees
## ------------------------------------------------------------

plot(get.tree(obj, 1))
plot(get.tree(obj, 2))
plot(get.tree(obj, 3))


## ------------------------------------------------------------
##
## Compare above to (1) user specified covariance matrix
##                  (2) default composite (independent) splitting
##
## ------------------------------------------------------------

## user specified sigma matrix
obj2 <- rfsrc(get.mv.formula(colnames(ydta)),
              data.frame(ydta, xdta),
              importance = TRUE, nsplit = 10,
              splitrule = "mahalanobis",
              sigma = cov(ydta))
print(obj2)

## default independence split rule
obj3 <- rfsrc(get.mv.formula(colnames(ydta)),
              data.frame(ydta, xdta),
              importance=TRUE, nsplit = 10)
print(obj3)

## compare vimp
imp <- data.frame(mahalanobis  = rowMeans(get.mv.vimp(obj,  standardize = TRUE)),
                  mahalanobis2 = rowMeans(get.mv.vimp(obj2, standardize = TRUE)),
                  default      = rowMeans(get.mv.vimp(obj3, standardize = TRUE)))

print(head(100 * imp[order(imp$mahalanobis, decreasing = TRUE), ], 15))

Compute Partial Dependence Values

Description

Compute observation-level predictions after setting a predictor to specified values. The saved forest is reused. Results can be averaged for partial dependence curves or retained for custom plots. Regression, classification, multivariate and mixed outcomes, right-censored survival, and competing risks have family-specific outputs.

Usage

partial.rfsrc(object, oob = TRUE,
  partial.type = NULL, partial.xvar = NULL, partial.values = NULL,
  partial.xvar2 = NULL, partial.values2 = NULL,
  partial.time = NULL, get.tree = NULL, seed = NULL, do.trace = FALSE, ...)

get.partial.plot.data(o, target, m.target = NULL, granule = FALSE)

Arguments

Details

Partial predictions and averaging

For a primary predictor X_j, the calculation replaces its value by z for every training observation while leaving the other predictor values at that observation's values. Write the resulting prediction as P_i(z) = \widehat f(z, X_{i,-j}). The default summary returned by get.partial.plot.data() is the average of the available P_i(z) values. With OOB prediction, the prediction function for row i uses that row's OOB trees.

partial() returns the observation-level arrays and does not average them. The extraction helper replaces infinite values by NA and omits missing values when averaging. A grid value with no available predictions has an unavailable mean.

Predictor grids and factors

The returned partial.values gives the actual sorted, unique grid. The input order and duplicated values are not retained. Reducing the grid reduces the number of partial predictions and the size of the result.

Factor codes refer to the level order saved during training. For example, if the original levels are c("low", "middle", "high"), the code for "high" is 3 even when a plotting subset contains only "middle" and "high". Use match(labels, levels(training.factor)) to obtain these codes. The extraction helper returns numeric grid values, so factor labels can be added in the custom plotting call.

Fixing additional predictors

partial.xvar2 and partial.values2 fix one or more additional coordinates for every observation. For example, vary wind over a grid while holding temperature at a specified value. To make a two-dimensional display, repeat the call at several fixed temperature values. A single call does not expand partial.values2 into a Cartesian product of grids.

Responses and survival summaries

Regression outputs contain predicted responses. Classification outputs contain an "all" column followed by columns named for the class levels. Multivariate and mixed fits return named response components; m.target in the extraction helper selects among these components. The class selector target is then resolved using that response's class labels, not response names. A single observation or grid point retains its observation-by-grid meaning during extraction.

For right-censored survival, "mort" returns mortality, "surv" returns survival probabilities, and "chf" returns cumulative hazards. At this low-level interface, "rel.freq" requests the same raw output as "mort"; neither partial() nor get.partial.plot.data() normalizes this output.

Competing-risk outputs are expected years lost, cumulative incidence, or cause-specific cumulative hazard, selected by partial.type. The time-dependent outputs are calculated on the forest's stored time grid and then selected in R for partial.time. The returned partial.time retains the requested times.

Custom plots

With granule = FALSE, the extraction helper returns one mean for each primary grid value. Time-dependent survival summaries have one column per requested time. With granule = TRUE, values are ordered by primary grid value and then observation; the result contains raw predictions, without the display shrinkage used for boxplots by plot.variable(partial = TRUE).

The automatic plotting wrapper also returns partial.summary, containing the partial mean, empirical row-averaging standard error, and available row count for each grid value. See plot.variable.rfsrc for that scale's interpretation.

Use plot.variable() for an automatically selected predictor grid and built-in displays. Use partial() when explicit grids, several evaluation times, additional fixed predictors, or raw observation-level results are needed.

Value

partial() returns an object of class (rfsrc, partial, family). It contains call, family, partial.values, partial.time, yvar.names, and event.info, together with the applicable prediction components below. Let n be the number of training rows, K the number of sorted unique primary values, T the number of requested times, and J the number of events.

get.partial.plot.data() returns a list with x, yhat, and partial.time. Normally x is the primary grid and yhat is a vector of means or a grid-by-time matrix. With granule = TRUE, x is repeated across observations and yhat contains the corresponding raw values, with one column per time for time-dependent summaries.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

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

See Also

plot.variable.rfsrc

Examples

## ------------------------------------------------------------
##
## Regression: explicit grid and raw observation-level output.
##
## ------------------------------------------------------------
set.seed(19)
dta <- na.omit(airquality)
fit <- rfsrc(Ozone ~ ., data = dta, ntree = 100)
wind <- unique(as.numeric(quantile(fit$xvar$Wind, seq(0, 1, length.out = 11))))
p <- partial(fit, partial.xvar = "Wind", partial.values = wind)
pd <- get.partial.plot.data(p)
print(head(data.frame(Wind = pd$x, partial = pd$yhat)))
plot(pd$x, pd$yhat, type = "b", xlab = "Wind", ylab = "Partial prediction")
raw <- get.partial.plot.data(p, granule = TRUE)
boxplot(raw$yhat ~ raw$x, xlab = "Wind", ylab = "Individual partial predictions")

## ------------------------------------------------------------
##
## Vary wind while fixing temperature at three values.
##
## ------------------------------------------------------------
temperature <- as.numeric(quantile(fit$xvar$Temp, c(.25, .5, .75)))
slices <- lapply(temperature, function(tt) {
  pp <- partial(fit, partial.xvar = "Wind", partial.values = wind,
                partial.xvar2 = "Temp", partial.values2 = tt)
  get.partial.plot.data(pp)$yhat
})
matplot(wind, do.call(cbind, slices), type = "l", lty = 1,
        xlab = "Wind", ylab = "Partial prediction")
legend("topright", legend = paste("Temp =", temperature),
       col = seq_along(temperature), lty = 1)

## ------------------------------------------------------------
##
## regression: partial effects for three variables simultaneously
## (can be slow, so modify accordingly)
##
## ------------------------------------------------------------

n <- 200
x <- matrix(rnorm(n * 3), ncol = 3)
y <- x[, 1] + x[, 1] * x[, 2] + x[, 1] * x[, 2] * x[, 3]
o <- rfsrc(y ~ ., data = data.frame(y = y, x))

## define target x values
x1 <- seq(-3, 3, length = 40)
x2 <- x3 <- seq(-3, 3, length = 10)

## extract second order partial effects
pdta <- do.call(rbind,
          lapply(x3, function(x3v) {
            cat("outer loop x3 = ", x3v, "\n")
            do.call(rbind,lapply(x2, function(x2v) {
              o <- partial(o,
                      partial.xvar = "X1",
                      partial.values = x1,
                      partial.xvar2 = c("X2", "X3"),
                      partial.values2 = c(x2v, x3v))
              cbind(x1, x2v, x3v, get.partial.plot.data(o)$yhat)
            }))
          }))
pdta <- data.frame(pdta)
colnames(pdta) <- c("x1", "x2", "x3", "effectSize")

## coplot of partial effects
coplot(effectSize ~ x1|x2*x3, pdta, pch = 16, overlap = 0)


## ------------------------------------------------------------
##
## Classification: extract a named class from a univariate response.
##
## ------------------------------------------------------------
iris.fit <- rfsrc(Species ~ ., data = iris, ntree = 100)
sp <- partial(iris.fit, partial.xvar = "Sepal.Length",
              partial.values = seq(4.5, 7.5, length.out = 11))
sp.data <- get.partial.plot.data(sp, target = "versicolor")
plot(sp.data$x, sp.data$yhat, type = "b", ylim = c(0, 1),
     xlab = "Sepal length", ylab = "Partial probability of versicolor")

## ------------------------------------------------------------
##
## Multivariate regression: one partial call, two extracted responses.
##
## ------------------------------------------------------------
mv <- rfsrc(cbind(Ozone, Temp) ~ ., data = dta, ntree = 100)
mp <- partial(mv, partial.xvar = "Wind", partial.values = wind)
print(names(mp$regrOutput))
mp.ozone <- get.partial.plot.data(mp, m.target = "Ozone")
mp.temp <- get.partial.plot.data(mp, m.target = "Temp")
print(head(data.frame(Wind = mp.temp$x, Temp = mp.temp$yhat)))

## ------------------------------------------------------------
##
## Mixed outcomes: extract a named class from the selected response.
##
## ------------------------------------------------------------
car.mix <- mtcars
car.mix$am <- factor(car.mix$am)
mix <- rfsrc(Multivar(mpg, am) ~ ., data = car.mix, ntree = 100)
mix.p <- partial(mix, partial.xvar = "wt", partial.values = c(2, 3, 4))
mix.am <- get.partial.plot.data(mix.p, m.target = "am", target = "1")
print(data.frame(wt = mix.am$x, probability = mix.am$yhat))

## ------------------------------------------------------------
##
## Factors: keep the training level codes when choosing labels.
##
## ------------------------------------------------------------
car.data <- mtcars
car.data$cyl <- factor(car.data$cyl)
car.fit <- rfsrc(mpg ~ ., data = car.data, ntree = 100)
lev <- levels(car.fit$xvar$cyl)
code <- match(lev, levels(car.fit$xvar$cyl))
cp <- partial(car.fit, partial.xvar = "cyl", partial.values = code)
cd <- get.partial.plot.data(cp, granule = TRUE)
boxplot(cd$yhat ~ factor(cd$x, levels = code, labels = lev),
        xlab = "Cylinders", ylab = "Individual partial predictions")

## ------------------------------------------------------------
##
## Right-censored survival: several stored evaluation times.
##
## ------------------------------------------------------------
data(veteran, package = "randomForestSRC")
vfit <- rfsrc(Surv(time, status) ~ ., data = veteran, ntree = 100)
grid <- vfit$time.interest
when <- grid[unique(round(seq(1, length(grid), length.out = 3)))]
karno <- unique(as.numeric(quantile(vfit$xvar$karno)))
vp <- partial(vfit, partial.type = "surv", partial.xvar = "karno",
              partial.values = karno, partial.time = when)
vd <- get.partial.plot.data(vp)
print(vd$partial.time)
matplot(vd$x, vd$yhat, type = "l", lty = 1,
        xlab = "Karnofsky score", ylab = "Partial survival probability")
legend("bottomright", legend = paste("Time =", vd$partial.time),
       col = seq_along(vd$partial.time), lty = 1)

## ------------------------------------------------------------
##
## Survival curves at selected quantiles of Karnofsky score.
## Reuse vfit from the preceding right-censored survival example.
##
## ------------------------------------------------------------
karno.probs <- c(.10, .50, .90)
karno.values <- quantile(vfit$xvar$karno, probs = karno.probs,
                         na.rm = TRUE)
print(karno.values)

## Set Karnofsky score to each requested value for every reference row.
## Keep the other predictors at their observed values and use OOB trees.
surv.partial <- partial(vfit,
                         oob = TRUE,
                         partial.type = "surv",
                         partial.xvar = "karno",
                         partial.values = as.numeric(karno.values),
                         partial.time = vfit$time.interest)
surv.data <- get.partial.plot.data(surv.partial)

## yhat has one row per distinct Karnofsky value and one column per time.
## Transpose it so that each plotted column is a survival curve over time.
## Each curve averages over the same reference rows.
curve.id <- seq_along(surv.data$x)
matplot(surv.data$partial.time, t(surv.data$yhat),
        type = "s", col = curve.id, lty = curve.id, lwd = 2,
        ylim = c(0, 1), xlab = "Time", ylab = "Survival probability",
        main = "Partial survival curves by Karnofsky score")

## Use the returned grid for labels: partial() removes tied quantile values.
legend("topright",
       legend = paste("Karnofsky =", format(surv.data$x, trim = TRUE)),
       col = curve.id, lty = curve.id, lwd = 2, bty = "n")

## ------------------------------------------------------------
##
## Competing risks: select an event code from the saved event map.
##
## ------------------------------------------------------------
data(follic, package = "randomForestSRC")
crfit <- rfsrc(Surv(time, status) ~ ., data = follic, ntree = 100, nsplit = 3)
age <- unique(as.numeric(quantile(crfit$xvar$age)))
crp <- partial(crfit, partial.type = "years.lost",
               partial.xvar = "age", partial.values = age)
event <- crp$event.info$event.type[1]
crd <- get.partial.plot.data(crp, target = event)
plot(crd$x, crd$yhat, type = "b", xlab = "Age", ylab = "Partial years lost")

Primary Biliary Cirrhosis (PBC) Data

Description

Data from the Mayo Clinic trial in primary biliary cirrhosis (PBC) of the liver conducted between 1974 and 1984. A total of 424 PBC patients, referred to Mayo Clinic during that ten-year interval, met eligibility criteria for the randomized placebo controlled trial of the drug D-penicillamine. The first 312 cases in the data set participated in the randomized trial and contain largely complete data.

Source

Flemming and Harrington, 1991, Appendix D.1.

References

Flemming T.R and Harrington D.P., (1991) Counting Processes and Survival Analysis. New York: Wiley.

Examples

data(pbc, package = "randomForestSRC")
pbc.obj <- rfsrc(Surv(days, status) ~ ., pbc, nsplit = 3)

Systolic Heart Failure Data

Description

The data involve 2231 patients with systolic heart failure who underwent cardiopulmonary stress testing at the Cleveland Clinic. The primary end point was all-cause death. In total, 39 variables were measured for each patient, including baseline clinical values and exercise stress test results. A key variable of interest is peak VO2 (mL/kg per min), the peak respiratory exchange ratio. More details regarding the data can be found in Hsich et al. (2011).

References

Hsich E., Gorodeski E.Z.,Blackstone E.H., Ishwaran H. and Lauer M.S. (2011). Identifying important risk factors for survival in systolic heart failure patients using random survival forests. Circulation: Cardio. Qual. Outcomes, 4(1), 39-45.

Examples

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

## random survival forest analysis
o <- rfsrc(Surv(ttodead, died)~., peakVO2)
print(o)

## partial effect of peak V02 on mortality
partial.o <- partial(o,
       partial.type = "mort",
       partial.xvar = "peak.vo2",
       partial.values = o$xvar$peak.vo2,
       partial.time = o$time.interest)
pdta.m <- get.partial.plot.data(partial.o)


## partial effect of peak V02 on survival
pvo2 <- quantile(o$xvar$peak.vo2)
partial.o <- partial(o,
       partial.type = "surv",
       partial.xvar = "peak.vo2",
       partial.values = pvo2,
       partial.time = o$time.interest)
pdta.s <- get.partial.plot.data(partial.o)
     

## compare the two plots
par(mfrow=c(1,2))    

plot(lowess(pdta.m$x, pdta.m$yhat, f = 2/3),
     type = "l", xlab = "peak VO2", ylab = "adjusted mortality")
rug(o$xvar$peak.vo2)

matplot(pdta.s$partial.time, t(pdta.s$yhat), type = "l", lty = 1,
          xlab = "years", ylab = "peak VO2 adjusted survival")
legend("bottomleft", legend = paste0("peak VO2 = ", pvo2),
       bty = "n", cex = .75, fill = 1:5)

Plots for Competing Risks

Description

Plot useful summary curves from a random survival forest competing risk analysis.

Usage

## S3 method for class 'rfsrc'
plot.competing.risk(x, plots.one.page = FALSE, ...)

Arguments

Details

Given a random survival forest object from a competing risk analysis (Ishwaran et al. 2014), plots from top to bottom, left to right: (1) cause-specific cumulative hazard function (CSCHF) for each event, (2) cumulative incidence function (CIF) for each event, and (3) continuous probability curves (CPC) for each event (Pepe and Mori, 1993).

Does not apply to right-censored data. Whenever possible, out-of-bag (OOB) values are displayed.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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.

Pepe, M.S. and Mori, M., (1993). Kaplan-Meier, marginal or conditional probability curves in summarizing competing risks failure time data? Statistics in Medicine, 12(8):737-751.

See Also

follic, hd, rfsrc, wihs

Examples

## ------------------------------------------------------------
## follicular cell lymphoma
## ------------------------------------------------------------

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

## ------------------------------------------------------------
## Hodgkin's Disease
## ------------------------------------------------------------

  data(hd, package = "randomForestSRC")
  hd.obj <- rfsrc(Surv(time, status) ~ ., hd, nsplit = 3, ntree = 100)
  print(hd.obj)
  plot.competing.risk(hd.obj)

## ------------------------------------------------------------
## competing risk analysis of pbc data from the survival package
## events are transplant (1) and death (2)
## ------------------------------------------------------------

if (library("survival", logical.return = TRUE)) {
   data(pbc, package = "survival")
   pbc$id <- NULL
   plot.competing.risk(rfsrc(Surv(time, status) ~ ., pbc))
}

Plot Conditional Quantiles and CRPS Diagnostics

Description

Plots observed responses against their conditional quantiles. An optional inset compares the forest's continuous ranked probability score (CRPS) curve with a predictor-free reference based on the training responses.

Usage

## S3 method for class 'rfsrc'
plot.quantreg(x, prbL = .25, prbU = .75,
  m.target = NULL, crps = TRUE, subset = NULL,
  xlab = NULL, ylab = NULL, ...,
  inset.args = list(), crps.null = TRUE, quantreg.tau = NULL)

Arguments

Details

The main panel places observed response values on the horizontal axis and the requested conditional quantiles on the vertical axis. A point marks the middle quantile, and a vertical segment with endpoint marks spans the lower and upper quantiles. The dashed diagonal is the identity line. Horizontal jitter reduces overplotting; it does not change the response values used for CRPS. Rows with unavailable responses or requested quantiles are omitted from the main panel with a warning.

An explicit quantreg.tau adds pinball losses to the main panel. Each loss is averaged over the requested subset, excluding unavailable responses and quantiles separately at that level. It uses the original response values, not their jittered display positions. These are the same values returned by get.pinball.error(x, tau = quantreg.tau, subset = subset, m.target = m.target) for the selected response. The reporting levels are independent of the lower and upper quantiles shown in the plot. The session option quantreg.tau sets defaults for printing and loss extraction; plotting keeps this annotation off unless levels are supplied explicitly.

The inset shows the average squared CDF error, integrated by the trapezoidal rule from the first response-grid value to each successive grid value and divided by the width of that interval. Its horizontal axis is the response threshold and its right-hand vertical axis is standardized CRPS. Lower values indicate less error over the corresponding integration interval. This is the finite-grid curve computed by get.quantile.crps(), with the same dependence on the reporting grid. The first value is unavailable because its integration interval has zero width. An entirely unavailable curve is omitted with a warning.

The null reference uses the empirical CDF of all finite training responses for the selected outcome,

F_0(t) = \frac{1}{n_0}\sum_{j=1}^{n_0} I(Y_j^{\mathrm{train}}\leq t),

where n_0 is the number of finite training responses. Repeated response values retain their frequencies. This same predictor-free CDF is used for every evaluated observation. Training responses are read from x$forest$yvar, or from x$yvar for a grow object without saved training responses in its forest. Test responses are used only for evaluation, not to construct the reference. Distinct reporting-grid values are not a replacement for the training sample.

Both curves use the same response grid, subset, integration rule, and standardization. At each threshold, both scores exclude the same rows with an unavailable observed response or forest CDF prediction. The reference distribution itself is not restricted by subset. In training plots, the full-training empirical reference includes the evaluated training response; it is an in-sample benchmark, not an OOB or leave-one-out null estimate. If training responses are unavailable, the null curve is omitted with a warning.

Graphical arguments in ... control the main panel. Customize the inset separately with inset.args, for example, inset.args = list(ylim = c(0, 0.3), lwd = 1.5). Default inset limits cover both curves. Limits supplied inside inset.args change the displayed inset window only; they do not recompute or truncate either score integral. Main-panel limits also leave both score calculations unchanged. Inset graphics settings are restored before returning so additional drawing refers to the main panel, and successive calls can be used in a multi-panel layout.

Value

Used for its graphical side effect. Invisibly returns NULL.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

See Also

quantreg, get.quantile(), get.quantile.crps(), plot.default

Examples

## Univariate quantiles, with forest and null CRPS curves.
library(randomForestSRC)
set.seed(19)
dta <- na.omit(airquality)
prob <- c(.25, .50, .75)
q <- quantreg(Temp ~ ., data = dta, prob = prob, ntree = 100)
plot.quantreg(q)

## Add pinball losses for this plot only; the interval and inset are unchanged.
plot.quantreg(q, quantreg.tau = c(.2, .5, .8))
print(get.pinball.error(q, tau = c(.2, .5, .8)))

## Customize the main panel directly, and the inset separately.
plot.quantreg(q, main = "Temperature",
              xlim = c(50, 100), ylim = c(45, 110), pch = 19,
              inset.args = list(xlim = c(60, 95), ylim = c(0, .35),
                                lwd = 1.5))

## Keep the inset, but suppress its null reference.
plot.quantreg(q, crps.null = FALSE)

## Multivariate forest: select each response for plotting only.
mv <- quantreg(cbind(Ozone, Temp) ~ ., data = dta,
               splitrule = "mahalanobis", prob = prob, ntree = 100)
print(names(get.quantile(mv, pretty = FALSE)))
op <- par(mfrow = c(1, 2))
plot.quantreg(mv, m.target = "Ozone")
plot.quantreg(mv, m.target = "Temp", main = "Temperature",
              xlim = c(50, 100), ylim = c(45, 110),
              inset.args = list(ylim = c(0, .35)))
par(op)

Plot Error Rate and Variable Importance from a RF-SRC analysis

Description

Plot out-of-bag (OOB) error rates and variable importance (VIMP) from a RF-SRC analysis. This is the default plot method for the package.

Usage

## S3 method for class 'rfsrc'
plot(x, m.target = NULL,
  plots.one.page = TRUE, sorted = TRUE, verbose = TRUE,  ...)

Arguments

Details

Plot cumulative OOB error rates as a function of number of trees and variable importance (VIMP) if available. Note that the default settings are now such that the error rate is no longer calculated on every tree and VIMP is only calculated if requested. To get OOB error rates for ever tree, use the option block.size = 1 when growing or restoring the forest. Likewise, to view VIMP, use the option importance when growing or restoring the forest.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

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

Examples

## ------------------------------------------------------------
## classification example
## ------------------------------------------------------------

iris.obj <- rfsrc(Species ~ ., data = iris,
     block.size = 1, importance = TRUE)
plot(iris.obj)

## ------------------------------------------------------------
## competing risk example
## ------------------------------------------------------------

## use the pbc data from the survival package
## events are transplant (1) and death (2)
if (library("survival", logical.return = TRUE)) {
  data(pbc, package = "survival")
  pbc$id <- NULL
  plot(rfsrc(Surv(time, status) ~ ., pbc, block.size = 1))
}

## ------------------------------------------------------------
## multivariate mixed forests
## ------------------------------------------------------------

mtcars.new <- mtcars
mtcars.new$cyl <- factor(mtcars.new$cyl)
mtcars.new$carb <- factor(mtcars.new$carb, ordered = TRUE)
mv.obj <- rfsrc(cbind(carb, mpg, cyl) ~., data = mtcars.new, block.size = 1)
plot(mv.obj, m.target = "carb")
plot(mv.obj, m.target = "mpg")
plot(mv.obj, m.target = "cyl")

Plot Subsampling Confidence Intervals for Variable Importance

Description

Display confidence intervals for variable importance (VIMP) from a subsample object. Select predictors, an outcome, and an interval method using the saved replicate estimates. Prediction-error and joint-VIMP rows can also be displayed when requested during subsampling.

Usage

## S3 method for class 'rfsrc'
plot.subsample(x, alpha = .01, xvar.names,
 standardize = TRUE, normal = TRUE, jknife = FALSE, target, m.target = NULL,
 pmax = 75, main = "", sorted = TRUE, show.plots = TRUE, ...)

Arguments

Details

Choosing an interval display

The default uses normal intervals and the subsampling standard error. Use normal = TRUE, jknife = TRUE for normal intervals with the jackknife standard error, or normal = FALSE for nonparametric intervals. All use the estimates already stored in x; changing alpha or the display method does not perform another round of subsampling.

The subsampling standard error measures replicate dispersion around the replicate mean. The jackknife calculation measures dispersion around the full-data estimate, incorporating the displacement between those two centers. The jackknife standard error is not necessarily larger for a finite replicate set. The nonparametric interval reverses empirical quantiles of the centered subsampling roots. See subsample.rfsrc for the calculations.

The display summarizes uncertainty in VIMP, not the distribution of the observed predictor values. Its intervals are calculated separately for each displayed statistic. A positive lower VIMP endpoint identifies a positive-importance result under that interval procedure; the display does not apply a multiple-testing correction. Prediction-error rows summarize the error itself.

Selecting predictors and outcomes

xvar.names selects rows and m.target selects a response. For a classification response or competing-risk outcome, target additionally selects its class-specific or event-specific statistic. These selections are made after resampling, so one saved object can be used for several displays.

Set alpha explicitly when comparing results with print or extract.subsample; their default is .05, whereas the plotting default is .01.

Graphical customization

The default display is horizontal. xlim controls its horizontal VIMP/error axis and ylim controls its vertical row positions. These arguments always refer to the displayed axes; the wrapper handles the internal limit reversal used by bxp for horizontal boxes. horizontal = FALSE gives a vertical display, with VIMP/error on the vertical axis.

Use boxfill (or col) for the box fill, border for its border, and whisklty and whisklwd for interval whiskers. By default, boxes are red when their lower endpoint is positive and blue otherwise. User styles are retained when the intervals are redrawn above the guides. outline = FALSE remains the default. The whiskers are confidence limits and the boxes show the inner interval summaries, rather than quartiles of the observed predictor values.

cex.axis, col.axis, and las control axis labels. Set xaxt = "n" or yaxt = "n" to suppress one axis, or axes = FALSE to suppress both. A scalar ylab is a vertical-axis title. names supplies one row label per selected statistic before sorting and trimming; labels follow the same ordering as the statistics. at supplies positions for the final displayed intervals.

show.plots = FALSE returns the plot data invisibly without opening a graphics device. extract.subsample(x, raw = TRUE) returns the interval matrices and replicate estimates directly.

Value

Invisibly returns a boxplot-summary list. Its stats component is the selected five-row confidence-interval matrix, in displayed order, and names contains the row labels. The remaining components are the boxplot scaffold, not additional subsampling confidence intervals. The same list is returned with show.plots = FALSE.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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.

Politis, D.N. and Romano, J.P. (1994). Large sample confidence regions based on subsamples under minimal assumptions. The Annals of Statistics, 22(4):2031-2050.

Shao, J. and Wu, C.J. (1989). A general theory for jackknife variance estimation. The Annals of Statistics, 17(3):1176-1197.

See Also

subsample.rfsrc, extract.subsample, bxp

Examples

## Small settings are for illustration; increase B for final inference.
set.seed(19)
dta <- na.omit(airquality)
o <- rfsrc(Ozone ~ ., data = dta, ntree = 100,
            importance = "permute", block.size = 1)
smp <- subsample(o, B = 25, verbose = FALSE)

## Three interval displays, all at the same confidence level.
plot.subsample(smp, alpha = .05, main = "Subsampling normal intervals")
plot.subsample(smp, alpha = .05, jknife = TRUE,
                 main = "Jackknife normal intervals")
plot.subsample(smp, alpha = .05, normal = FALSE,
                 main = "Nonparametric intervals")

## Restrict the display and customize the main axes and whiskers.
plot.subsample(smp, alpha = .05,
                 xvar.names = c("Solar.R", "Wind", "Temp"),
                 xlim = c(-.1, .8), las = 2, cex.axis = .75,
                 whisklty = 1, whisklwd = 1.5)
plot.data <- plot.subsample(smp, alpha = .05, show.plots = FALSE)
print(plot.data)

## Multivariate regression: one subsample bank, two response displays.
mv <- rfsrc(cbind(Ozone, Temp) ~ ., data = dta, ntree = 100,
             importance = "permute", block.size = 1)
mv.smp <- subsample(mv, B = 25, verbose = FALSE)
plot.subsample(mv.smp, m.target = "Ozone", alpha = .05,
                 main = "Ozone")
plot.subsample(mv.smp, m.target = "Temp", alpha = .05,
                 main = "Temperature")
print(extract.subsample(mv.smp, m.target = "Temp", alpha = .05)$var.sel.Z)

Survival Prediction Performance and Diagnostic Plots

Description

Plot survival estimates and calculate inverse-probability-of-censoring weighted Brier scores and cumulative/dynamic time-dependent AUC. plotBrierAUC also compares the forest Brier curve with a predictor-free Kaplan–Meier reference estimated from the training data. get.cindex calculates concordance-based prediction error from observed outcomes and risk predictions for right-censored survival or competing risks.

Usage

## S3 method for class 'rfsrc'
plot.survival(x, show.plots = TRUE, subset,
  collapse = FALSE, cens.model = c("km", "rfsrc"), ...)

get.brier.survival(o, subset = NULL,
  cens.model = c("km", "rfsrc"), papply = lapply,
  times = NULL, conf.int = FALSE, keep.matrix = TRUE)

get.auct.survival(o, subset = NULL,
  cens.model = c("km", "rfsrc"), papply = lapply,
  times = NULL, conf.int = FALSE)

plotBrierAUC(x, subset = NULL,
  cens.model = c("km", "rfsrc"), papply = lapply,
  times = NULL, conf.int = TRUE,
  plots = c("brier", "auct"), show.plots = TRUE,
  brier.null = TRUE, ...)

get.cindex(time, censoring, predicted, weight, fast, do.trace = FALSE)

Arguments

Details

plot.survival produces the following plots, going from top to bottom and left to right:

1. Forest estimated survival function for each individual. The thick red line is the ensemble mean survival and the thick green line is the marginal Nelson–Aalen survival estimate.

2. Brier score stratified by ensemble mortality. Stratification is into four groups corresponding to the 0–25, 25–50, 50–75 and 75–100 percentile ranges of mortality. The red line is the overall Brier score.

3. Cumulative integrated Brier score divided by elapsed time, labeled CRPS in the plot.

4. Mortality (Ishwaran et al., 2008) versus observed time. Blue points are events and black points are censored observations. Mortality is estimated risk calibrated to the scale of the number of events. For example, a mortality value of 100 means that if all individuals had the same covariate values, an average of 100 events would be expected.

Whenever possible, out-of-bag predictions are used. For a prediction object, the predicted test survival curves and test outcomes are used for scoring, while the censoring distribution is still estimated from the grow data. If a prediction object contains no outcomes, plot.survival displays only its predicted survival curves. Brier score, AUC, CRPS, and mortality-versus-time cannot be calculated without evaluation outcomes. Direct calls to get.brier.survival, get.auct.survival, or plotBrierAUC therefore require outcomes for the cases being evaluated.

For subject i at horizon t, the IPCW Brier contribution is

L_i(t) = \frac{I(T_i \le t, \Delta_i > 0)}{\widehat G(T_i-)} \widehat S_i(t)^2 + \frac{I(T_i > t)}{\widehat G(t)} \{1-\widehat S_i(t)\}^2,

with subject-specific values of \widehat G when cens.model="rfsrc". The Brier score is the sample mean of these contributions.

The time-dependent AUC uses a cumulative/dynamic definition. Cases at t are observed events satisfying T_i \le t; controls satisfy T_i > t; observations censored at or before t are not included in the case-control comparison. The time-specific risk score is 1-\widehat S_i(t), and comparisons are weighted by the same grow-data censoring distribution.

When confidence intervals are requested, the Brier standard error is obtained by deleting one subject-level IPCW loss at a time while holding the survival predictions and censoring weights fixed.

The AUC standard error uses a stratified delete-one jackknife. One observed case or control is deleted at a time, the remaining IPCW weights within that stratum are renormalized, and the resulting case-delete and control-delete variance components are added. With equal weights this calculation reduces exactly to the ordinary DeLong variance. The effective case and control sample sizes and the largest normalized IPCW weights are returned to diagnose horizons at which a few observations dominate the weighted comparison.

Both are pointwise conditional standard errors, and the reported intervals are pointwise normal approximations rather than simultaneous bands. The survival predictions and the estimated censoring distribution are treated as fixed. They do not account for uncertainty from training the survival forest or estimating the censoring model. For an independent prediction sample, the intervals have a direct conditional test-performance interpretation. For a grow object, shared out-of-bag fits and the shared censoring estimate induce dependence among subject-level quantities, so the intervals should be viewed as fixed-fit working intervals rather than repeated-training confidence intervals. A score and its interval are returned as NA at a horizon where a required censoring survival probability is zero. The AUC standard error is also NA if fewer than two cases or two controls are available, or if deleting one subject leaves no positive IPCW weight in its stratum.

plotBrierAUC is a base-R plotting helper. It calls a shared numerical engine so that Brier score, its null reference, and AUC use the same censoring model, then adds the requested pointwise confidence bands for the forest curves. The AUC panel retains its horizontal reference at 0.5.

Null reference for the Brier score

With brier.null=TRUE, the Brier panel includes a dashed predictor-free reference and a Forest/Null legend. Its survival prediction is the same for every evaluated subject:

\widehat S_0(t) = \prod_{u \le t}\left\{1-\frac{d_{\mathrm{grow}}(u)} {Y_{\mathrm{grow}}(u)}\right\},

where d_{\mathrm{grow}}(u) and Y_{\mathrm{grow}}(u) are the event count and risk-set size at grow-data event time u. All finite grow outcomes are used, including their frequencies at tied times. Events are processed before censorings at a tied time. This is the event-survival estimate; cens.model separately controls the censoring distribution used for scoring.

The null Brier curve is obtained by replacing \widehat S_i(t) with \widehat S_0(t) in the IPCW loss above. It uses the same horizons, evaluation subset, censoring weights, and finite forest-loss contributions as the forest curve. Censored observations at or before a horizon retain their zero contribution in both averages. A reference horizon is NA if the forest horizon is unsupported or a required null loss is nonfinite. Confidence bands are drawn only for the forest.

For test evaluation, the reference is estimated from the saved grow outcomes, not from the test outcomes. Changing subset changes the cases being scored, not the grow-data reference. For training evaluation, the reference uses the full training sample and is an in-sample benchmark, even when the forest predictions are OOB. The last Kaplan–Meier estimate is carried forward after its last observed time, subject to the existing censoring-support checks.

A forest curve below the null curve has lower estimated prediction error than this predictor-free reference at that horizon. The two curves are compared directly on the Brier scale; the score is not rescaled and no comparison test is added.

Only right-censored survival families are supported. Competing-risk analyses should use plot.competing.risk.

Concordance error

get.cindex accepts aligned outcome and risk-prediction values directly. Supply predicted.oob for OOB evaluation or predicted with the corresponding outcomes for new-data evaluation. It uses the supplied values without fitting a forest or selecting an OOB component automatically.

A status code greater than one selects the competing-risk branch. For that branch, J is the largest usable status code and prediction columns are indexed by event code, not matched by their names. Without weights, the calculation for event j uses observations censored or experiencing event j; observations with other events are excluded. With weights, it uses the native event-specific weighted concordance calculation instead.

Missing times, statuses, predictions, and supplied weights are excluded from the applicable calculation. The helper returns the native concordance error; it is distinct from the time-specific cumulative/dynamic AUC returned by get.auct.survival.

Value

get.brier.survival returns a list containing brier.score, the optional brier.matx, integrated scores crps and crps.std, the estimated censoring distribution, aligned grow and evaluation event information, survival predictions, mortality, and the selected subset. When confidence intervals are requested, brier.score also contains std.err, lower, upper, and n.eval.

get.auct.survival returns a corresponding list whose auct.score contains time, auct, n.case, n.control, effective sample sizes n.case.eff and n.control.eff, and largest normalized weights max.case.weight and max.control.weight. The columns std.err, lower, and upper are included when confidence intervals are requested.

plotBrierAUC invisibly returns a list containing the requested Brier and/or AUC result objects. When Brier is requested with brier.null=TRUE, the Brier result also contains a null list with brier.score (columns time and brier.score), integrated scores crps and crps.std, and the number of evaluated contributions n.eval. It also contains the common null survival vector on the scoring grid, n.train (finite grow outcomes), method="kaplan-meier", and source="grow". No subject-by-time null loss matrix or null confidence band is returned. The integrated null scores use the same trapezoidal rule and standardization as the forest scores. Direct numerical-helper calls retain their existing return structure.

With evaluation outcomes, plot.survival invisibly returns the mortality-stratified and overall Brier and cumulative integrated Brier curves. Without evaluation outcomes, it invisibly returns the predicted survival curves and their ensemble mean.

get.cindex returns a numeric concordance-error value for right-censored survival, or a numeric vector in event-code order 1, ..., J for competing risks. In the unweighted competing-risk calculation an event with fewer than two usable observations returns NA.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

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

DeLong E.R., DeLong D.M. and Clarke-Pearson D.L. (1988). Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach, Biometrics, 44:837–845.

Efron B. and Tibshirani R.J. (1993). An Introduction to the Bootstrap. Chapman and Hall, New York.

Heagerty P.J. and Zheng Y. (2005). Survival model predictive accuracy and ROC curves, Biometrics, 61:92–105.

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

Ishwaran H., Kogalur U.B., Blackstone E.H. and Lauer M.S. (2008). Random survival forests, Annals of Applied Statistics, 2:841–860.

See Also

plot.competing.risk.rfsrc, predict.rfsrc, rfsrc

Examples

## veteran data
data(veteran, package = "randomForestSRC")
plot.survival(rfsrc(Surv(time, status) ~ ., veteran),
              cens.model = "rfsrc")

## pbc data
data(pbc, package = "randomForestSRC")
pbc.obj <- rfsrc(Surv(days, status) ~ ., pbc)

## ------------------------------------------------------------
## Concordance error from the stored OOB risk predictions
## ------------------------------------------------------------
print(get.cindex(pbc.obj$yvar[, 1], pbc.obj$yvar[, 2],
                 pbc.obj$predicted.oob))

## standard survival diagnostics
plot.survival(pbc.obj)
plot.survival(pbc.obj, subset = c(3, 10), collapse = TRUE)

## Brier and AUCT helpers with pointwise intervals
brier.obj <- get.brier.survival(pbc.obj, conf.int = TRUE)
print(head(brier.obj$brier.score))
auct.obj <- get.auct.survival(pbc.obj, conf.int = TRUE)
print(head(auct.obj$auct.score))

## compare two grow-data censoring models
brier.km <- get.brier.survival(pbc.obj, cens.model = "km")
brier.rf <- get.brier.survival(pbc.obj, cens.model = "rfsrc")
plot(brier.km$brier.score$time,
     brier.km$brier.score$brier.score, type = "s", col = 2,
     xlab = "Time", ylab = "Brier Score")
lines(brier.rf$brier.score$time,
      brier.rf$brier.score$brier.score, type = "s", col = 4)
legend("bottomright",
       legend = c("cens.model = km", "cens.model = rfsrc"),
       col = c(2, 4), lty = 1)

## Brier and AUCT curves; Brier includes the grow-data KM null reference
perf <- plotBrierAUC(pbc.obj)
print(head(data.frame(
  time = perf$brier$time,
  forest = perf$brier$brier.score$brier.score,
  null = perf$brier$null$brier.score$brier.score
)))
plotBrierAUC(pbc.obj, plots = "auct", conf.int = 0.90)
plotBrierAUC(pbc.obj, plots = "brier", conf.int = 0.90,
             ylim = c(0, .4), null.lty = 3, legend.pos = "topleft")

## Omit the Brier reference, retaining the forest curve and its band
plotBrierAUC(pbc.obj, plots = "brier", brier.null = FALSE)

## Independent test evaluation: the null still uses training outcomes
set.seed(19)
pbc.complete <- na.omit(pbc)
trn <- sample(seq_len(nrow(pbc.complete)),
              size = floor(0.7 * nrow(pbc.complete)))
grow <- rfsrc(Surv(days, status) ~ ., pbc.complete[trn, ], ntree = 100)
test <- predict(grow, newdata = pbc.complete[-trn, ])
test.perf <- plotBrierAUC(test, plots = "brier")
print(c(forest = test.perf$brier$crps.std,
        null = test.perf$brier$null$crps.std))

## Obtain both curves without opening a graphics device
perf <- plotBrierAUC(test, plots = "brier", show.plots = FALSE)
print(head(perf$brier$null$brier.score))

Plot Marginal and Partial Dependence of Predictors

Description

Display ensemble predictions against predictors, either at their observed values (marginal plots) or after fixing one predictor at a sequence of values and averaging over the other predictors (partial plots). Plot data can be returned without drawing and reused for custom graphics.

Usage

## S3 method for class 'rfsrc'
plot.variable(x, xvar.names, target,
  m.target = NULL, time, surv.type = c("mort", "rel.freq",
  "surv", "years.lost", "cif", "chf"), class.type =
  c("prob", "bayes"), partial = FALSE, oob = TRUE,
  show.plots = TRUE, plots.per.page = 4, granule = 5, sorted = TRUE,
  nvar, npts = 25, smooth.lines = FALSE, subset, ...)

Arguments

Details

Marginal and partial plots

Marginal plots relate the ensemble's predictions to the observed value of a predictor. They preserve the observed combinations of predictors. Continuous predictors receive a smoothed curve; factors and numeric predictors with few distinct values receive grouped boxplots.

Partial plots fix the selected predictor at each grid value for every observation, retaining the other predictor coordinates, and summarize those predictions. For a grid value z, write the observation-level partial prediction as P_i(z)=\widehat f(z,X_{i,-j}). With complete predictions the curve is n^{-1}\sum_i P_i(z); missing predictions are omitted from the mean. When OOB output is requested, P_i(z) uses the OOB trees for row i.

Partial mode reuses the fitted trees but performs additional prediction work at each grid value. Restrict xvar.names or nvar and use a moderate npts to control the computation. For explicit grids or several fixed predictors, use partial.rfsrc.

Response selection

Regression displays the predicted response. Classification displays the selected class probability by default. The m.target argument chooses the response for a multivariate or mixed forest; target then chooses a class when that response is a factor.

Right-censored survival displays mortality, relative-frequency mortality, survival probability, or cumulative hazard. Competing risks displays expected years lost, cumulative incidence, or cause-specific cumulative hazard for the selected event. Survival probability, cumulative incidence, and cumulative hazard depend on time; mortality and years lost are scalar summaries and do not use the requested time.

The existing marginal display multiplies right-censored survival and cumulative hazard values by 100; partial displays retain their raw scales. Relative mortality also retains its existing mode-specific scaling: marginal mortality is divided by the maximum of the selected row count and the available full-vector mortality values; partial mortality is divided by the available row count at each grid value.

Subsets and missing values

The plotting calculation starts from the stored training predictors and overlays saved imputed predictor values when available. It does not grow a replacement forest. Define subsets relative to the stored row order. In partial mode, grid construction uses nonmissing values of the selected predictor within that subset; the same subset is passed to the partial-prediction extraction step.

A subset affects the population being summarized, not which observations trained the forest. Likewise, OOB selection changes the contributing trees for an observation, not the definition of the plotting subset.

Partial means and standard errors

A single calculation is used for regression predictions, class probabilities, and scalar survival or competing-risk predictions. At grid value z, let I_z contain the selected rows with finite partial predictions, and let n_z be their number. The mean and empirical row-averaging standard error are

\bar P(z)=\frac{1}{n_z}\sum_{i\in I_z}P_i(z),\qquad \mathrm{se}(z)=\left\{ \frac{\sum_{i\in I_z}[P_i(z)-\bar P(z)]^2}{n_z(n_z-1)} \right\}^{1/2}.

Thus se is the sample standard deviation divided by \sqrt{n_z}. Classification uses the observed variation among predicted probabilities, rather than a Bernoulli variance calculated from their mean. Both numerator and denominator use the same available rows. With no finite predictions the mean is NA; with fewer than two, the standard deviation and standard error are NA. Identical finite predictions have zero standard error when at least two are present.

This standard error describes reference-row averaging with a fixed prediction function under independent reference sampling. For the training/OOB plots here, it is a descriptive row-averaging scale: it does not account for fitting the forest or for dependence between overlapping OOB fits. The continuous display uses pointwise bands \bar P(z)\pm2\mathrm{se}(z), not refit-based confidence intervals or simultaneous confidence bands. The numeric class-code display uses this calculation on its code scale, not a probability scale.

Discrete partial displays

Discrete displays retain boxplots of the transformed values \bar P(z)+[P_i(z)-\bar P(z)]/\sqrt{n_z}. Their sample standard deviation is exactly se when n_z>1, so continuous and discrete displays use the same row-averaging scale. The boxes and whiskers remain ordinary boxplot summaries of these transformed values; their endpoints are not the mean plus or minus two standard errors. Notches likewise concern these display values, not uncertainty from refitting the forest. A single available value remains at its mean, with an unavailable standard error.

The returned plotthis tables retain these display values. partial.summary gives the means, standard deviations, standard errors, and available counts directly for every predictor, including factors. For unscaled observation-level predictions, use partial() and get.partial.plot.data(..., granule = TRUE).

Returning and reusing plot data

Assign the result of plot.variable(..., show.plots = FALSE) to retain the calculations. For partial mode, plotthis contains one data frame per predictor with columns x and yhat. Continuous curves normally have one row per grid value; discrete displays retain the shrunken values for each selected observation.

Call plot.variable(saved) to redraw without recalculating predictions. The saved object determines the marginal or partial mode, outcome, subset, grid, and display settings. Its plots.per.page and smooth.lines components can be changed before replotting. To change the prediction target, subset, or grid, create a new plotting object from the grow object.

Value

Invisibly returns an object of class (rfsrc, plot.variable, family). Common components include family, partial, event.info, target, ylabel, n, xvar.names, nvar, plots.per.page, granule, and smooth.lines.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Friedman J.H. (2001). Greedy function approximation: a gradient boosting machine, Ann. of Statist., 5:1189-1232.

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

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

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.

See Also

rfsrc, partial.rfsrc, predict.rfsrc

Examples

## ------------------------------------------------------------
##
## Regression: marginal versus partial display.
##
## ------------------------------------------------------------
set.seed(19)
dta <- na.omit(airquality)
fit <- rfsrc(Ozone ~ ., data = dta, ntree = 100, importance = TRUE)
plot.variable(fit, xvar.names = c("Wind", "Temp"), sorted = FALSE,
              plots.per.page = 2)
pv <- plot.variable(fit, xvar.names = c("Wind", "Temp"), sorted = FALSE,
                    partial = TRUE, npts = 15, plots.per.page = 2,
                    show.plots = FALSE)
print(head(pv$plotthis$Wind))
print(pv$partial.summary$Wind)
plot.variable(pv, main = "Partial dependence", lwd = 2)

## ------------------------------------------------------------
##
## Inspect the mean and its row-averaging standard error.
##
## ------------------------------------------------------------
ss <- pv$partial.summary$Wind
plot(ss$x, ss$mean, type = "b", xlab = "Wind", ylab = "Partial prediction")
lines(ss$x, ss$mean - 2 * ss$se, lty = 3)
lines(ss$x, ss$mean + 2 * ss$se, lty = 3)

## ------------------------------------------------------------
##
## Customize a stored continuous curve without another forest calculation.
##
## ------------------------------------------------------------
pv$plots.per.page <- 1
pv$smooth.lines <- TRUE
plot.variable(pv, xvar.names = "Temp", main = "Temperature")

## ------------------------------------------------------------
##
## Subset indices refer to the fitted object's stored rows.
##
## ------------------------------------------------------------
use <- fit$xvar$Solar.R < 200
ps <- plot.variable(fit, xvar.names = "Wind", partial = TRUE,
                    subset = use, npts = 15, show.plots = FALSE)
print(head(ps$plotthis$Wind))

## ------------------------------------------------------------
##
## Multivariate regression: choose a response for the plot.
##
## ------------------------------------------------------------
mv <- rfsrc(cbind(Ozone, Temp) ~ ., data = dta, ntree = 100)
mv.temp <- plot.variable(mv, m.target = "Temp", xvar.names = "Wind",
                         partial = TRUE, npts = 15, show.plots = FALSE)
print(mv.temp$partial.summary$Wind)
plot.variable(mv.temp)

## ------------------------------------------------------------
##
## Classification: choose a class label.
##
## ------------------------------------------------------------
iris.fit <- rfsrc(Species ~ ., data = iris, ntree = 100)
plot.variable(iris.fit, target = "versicolor", xvar.names = "Petal.Length",
              partial = TRUE, npts = 15)

## ------------------------------------------------------------
##
## A factor predictor receives a boxplot.
##
## ------------------------------------------------------------
car.data <- mtcars
car.data$cyl <- factor(car.data$cyl)
car.fit <- rfsrc(mpg ~ ., data = car.data, ntree = 100)
pc <- plot.variable(car.fit, xvar.names = "cyl", partial = TRUE,
                     show.plots = FALSE)
boxplot(yhat ~ x, data = pc$plotthis$cyl,
        xlab = "Cylinders", ylab = "Partial display values")

## ------------------------------------------------------------
##
## Mixed outcomes: m.target selects the response, target selects its class.
##
## ------------------------------------------------------------
car.data$am <- factor(car.data$am)
mix <- rfsrc(Multivar(mpg, am) ~ ., data = car.data, ntree = 100)
plot.variable(mix, m.target = "am", target = "1", xvar.names = "wt",
              partial = TRUE, npts = 15)

## ------------------------------------------------------------
##
## Right-censored survival at a stored evaluation time.
##
## ------------------------------------------------------------
data(veteran, package = "randomForestSRC")
vfit <- rfsrc(Surv(time, status) ~ ., data = veteran, ntree = 100)
when <- vfit$time.interest[ceiling(length(vfit$time.interest) / 2)]
plot.variable(vfit, xvar.names = "age", surv.type = "surv", time = when,
              partial = TRUE, npts = 15)

## ------------------------------------------------------------
##
## Competing risks: cumulative incidence for the second event.
##
## ------------------------------------------------------------
data(follic, package = "randomForestSRC")
crfit <- rfsrc(Surv(time, status) ~ ., data = follic, ntree = 100, nsplit = 3)
when <- crfit$time.interest[ceiling(length(crfit$time.interest) / 2)]
plot.variable(crfit, xvar.names = "age", target = 2, surv.type = "cif",
              time = when, partial = TRUE, npts = 15)

Prediction for Random Forests for Survival, Regression, and Classification

Description

Obtain predicted values using a forest. Also returns performance values if the test data contains y-outcomes.

Usage

## S3 method for class 'rfsrc'
predict(object,
  newdata,
  importance = c(FALSE, TRUE, "none", "anti", "permute", "random"),
  get.tree = NULL,
  block.size = if (any(is.element(as.character(importance),
                     c("none", "FALSE")))) NULL else 10,
  na.action = c("na.omit", "na.impute", "na.random"),
  outcome = c("train", "test"),
  perf.type = NULL,
  proximity = FALSE,
  forest.wt = FALSE,
  ptn.count = 0,
  distance = FALSE,
  var.used = c(FALSE, "all.trees", "by.tree"),
  split.depth = c(FALSE, "all.trees", "by.tree"),
  case.depth = FALSE,
  seed = NULL,
  do.trace = FALSE, membership = FALSE,
  marginal.xvar = NULL,  ...)

Arguments

Details

Predicted values are obtained by "dropping" the test data down the trained forest-i.e., the forest grown using the training data. If the test data includes y-outcome values, performance metrics are also returned. Variable importance (VIMP), including joint VIMP, is returned if requested.

If no test data is supplied, the function uses the original training data and enters "restore" mode. This allows users to extract outputs from the trained forest that were not requested during the original grow call.

When newdata is supplied, non-hybrid proximity and distance matrices compare test observations with one another. Setting either option to "hybrid" instead compares each retained test observation with each retained original training observation, producing an ntest by ntrain matrix without constructing a square matrix over the combined data. Hybrid mode should not be combined with marginal.xvar.

If outcome = "test", predictions are computed using y-outcomes from the test data (which must include outcome values). Terminal node statistics are recalculated using these outcomes, while the tree topology remains fixed from training. Error rates and VIMP are then computed by bootstrapping the test set and applying out-of-bagging to maintain unbiased estimates.

Set csv = TRUE to return case-specific VIMP, and cse = TRUE to return case-specific error rates. These apply to all families except survival. Both options can also be used at training time. Use get.mv.csvimp and get.mv.cserror to extract these stored values; see Extracting Multivariate Values (multivariate.values).

For prediction from a quantile forest, use quantreg.rfsrc with object and newdata. Its helpers include get.quantile for quantile extraction and get.pinball.error for mean pinball losses.

Value

An object of class (rfsrc, predict), which is a list with the following components:

call

The prediction call.

family

The family used in the analysis.

n

Sample size of the test data (after handling missing values).

ntree

Number of trees in the trained forest.

yvar

Y-outcome values from the test data or original grow data (if newdata is missing).

yvar.names

Character vector of response variable names.

xvar

Data frame of test set predictor variables.

xvar.names

Character vector of predictor variable names.

leaf.count

Vector of length ntree giving the number of terminal nodes per tree.

proximity

Proximity matrix computed on the test data. In hybrid mode, rows correspond to retained test observations and columns correspond to retained original training observations.

distance

Distance matrix. Hybrid mode uses the same test-by-training layout as hybrid proximity.

forest

The trained forest object.

forest.wt

Forest weight matrix for test cases.

ptn.membership

Matrix of pruned terminal node membership. Only returned if ptn.count > 0.

membership

Matrix of terminal node membership for test cases. Each column corresponds to one tree.

inbag

Matrix indicating how many times each case appears in the bootstrap sample for each tree.

var.used

Number of times each variable was used in splitting.

imputed.indv

Indices of test observations with missing values.

imputed.data

Imputed version of the test data. Columns are ordered with responses first, followed by predictors.

split.depth

NULL; split-depth calculations are currently disabled in this prediction interface.

err.rate

Prediction error, if y-outcomes are present. For ordinary prediction with newdata, this is test error; restoring the training data uses OOB error.

importance

Variable importance (VIMP) for the test data. May be NULL.

predicted

Predicted values for the test data.

predicted.oob

OOB predicted values. May be NULL depending on context.

quantile

Estimated quantile values at the requested probabilities (quantile regression only).

quantile.oob

OOB quantile values. May be NULL.

class

(Classification only) Predicted class labels.

class.oob

(Classification only) OOB predicted class labels.

regrOutput

(Multivariate only) Response-named list of predictions and performance measures for continuous outcomes.

classOutput

(Multivariate only) Response-named list of predictions and performance measures for categorical outcomes.

chf

(Survival or competing risks) Cumulative hazard function (CHF); cause-specific CHF for competing risks.

chf.oob

(Survival or competing risks) OOB CHF. May be NULL.

survival

(Survival only) Survival function estimates.

survival.oob

(Survival only) OOB survival function. May be NULL.

time.interest

(Survival or competing risks) Sorted unique event times.

ndead

(Survival or competing risks) Number of deaths observed.

cif

(Competing risks only) Cumulative incidence function (CIF) for each event type.

cif.oob

(Competing risks only) OOB CIF. May be NULL.

Note

The dimensions and contents of returned objects depend on the forest family and whether y-outcomes are available in the test data. In particular, performance-related components (e.g., error rate, VIMP) will be NULL if y-outcomes are missing.

For multivariate families, predicted values, VIMP, error rates, and performance metrics are stored in the lists regrOutput and classOutput. These can be accessed using the helper functions get.mv.predicted, get.mv.vimp, and get.mv.error. See Extracting Multivariate Values (multivariate.values) for OOB selection, standardization, block errors, and return formats.

Additional classification scores are documented under Classification Performance Metrics (classification.performance); get.imbalanced.performance provides binary-class summaries with selectable thresholds. For right-censored survival diagnostics, see plotBrierAUC.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

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

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

See Also

Classification Performance Metrics (classification.performance), Extracting Multivariate Values (multivariate.values), Fast Saving and Loading Objects (fast.saveload), holdout.vimp, imbalanced.rfsrc, impute.learn.rfsrc, plot.competing.risk.rfsrc, plot.rfsrc, plot.quantreg.rfsrc, plot.survival.rfsrc, plotBrierAUC, plot.variable.rfsrc, print.rfsrc, quantreg.rfsrc, rfsrc, rfsrc.fast, vimp.rfsrc

Examples

## ------------------------------------------------------------
## typical train/testing scenario
## ------------------------------------------------------------

data(veteran, package = "randomForestSRC")
train <- sample(1:nrow(veteran), round(nrow(veteran) * 0.80))
veteran.grow <- rfsrc(Surv(time, status) ~ ., veteran[train, ]) 
veteran.pred <- predict(veteran.grow, veteran[-train, ])
print(veteran.grow)
print(veteran.pred)

## ------------------------------------------------------------
## hybrid proximity and distance
## - rows are test cases; columns are training cases
## ------------------------------------------------------------

data(mtcars)
hybrid.train <- mtcars[1:24, ]
hybrid.test <- mtcars[25:32, ]
hybrid.obj <- rfsrc(mpg ~ ., data = hybrid.train, ntree = 50)

hybrid.proximity <- predict(hybrid.obj, hybrid.test,
                            proximity = "hybrid")$proximity
hybrid.distance <- predict(hybrid.obj, hybrid.test,
                           distance = "hybrid")$distance

## Both matrices have 8 test rows and 24 training columns.
print(dim(hybrid.proximity))
print(dim(hybrid.distance))


## ------------------------------------------------------------
## restore mode
## - if predict is called without specifying the test data
##   the original training data is used and the forest is restored
## ------------------------------------------------------------

## first train the forest
airq.obj <- rfsrc(Ozone ~ ., data = airquality)

## now we restore it and compare it to the original call
## they are identical
predict(airq.obj)
print(airq.obj)

## we can retrieve various outputs that were not asked for in
## in the original call

## here we extract the proximity matrix
prox <- predict(airq.obj, proximity = TRUE)$proximity
print(prox[1:10,1:10])

## here we extract the number of times a variable was used to grow
## the grow forest
var.used <- predict(airq.obj, var.used = "by.tree")$var.used
print(head(var.used))

## ------------------------------------------------------------
## prediction when test data has missing values
## ------------------------------------------------------------

data(pbc, package = "randomForestSRC")
trn <- pbc[1:312,]
tst <- pbc[-(1:312),]
o <- rfsrc(Surv(days, status) ~ ., trn)

## default imputation method used by rfsrc
print(predict(o, tst, na.action = "na.impute"))

## random imputation
print(predict(o, tst, na.action = "na.random"))

## ------------------------------------------------------------
## requesting different performance for classification
## ------------------------------------------------------------

## default performance is misclassification
o <- rfsrc(Species~., iris)
print(o)

## get (normalized) brier performance
print(predict(o, perf.type = "brier"))

## ------------------------------------------------------------
## vimp for each tree: illustrates get.tree 
## ------------------------------------------------------------

## regression analysis but no VIMP
o <- rfsrc(mpg~., mtcars)

## now extract VIMP for each tree using get.tree
vimp.tree <- do.call(rbind, lapply(1:o$ntree, function(b) {
     predict(o, get.tree = b, importance = TRUE)$importance
}))

## boxplot of tree VIMP
boxplot(vimp.tree, outline = FALSE, col = "cyan")
abline(h = 0, lty = 2, col = "red")

## summary information of tree VIMP
print(summary(vimp.tree))

## extract tree-averaged VIMP using importance=TRUE
## remember to set block.size to 1
print(predict(o, importance = TRUE, block.size = 1)$importance)

## use direct call to vimp() for tree-averaged VIMP
print(vimp(o, block.size = 1)$importance)

## ------------------------------------------------------------
## vimp for just a few trees
## illustrates how to get vimp if you have a large data set
## ------------------------------------------------------------

## survival analysis but no VIMP
data(pbc, package = "randomForestSRC")
o <- rfsrc(Surv(days, status) ~ ., pbc, ntree = 2000)

## get vimp for a small number of trees
print(predict(o, get.tree=1:250, importance = TRUE)$importance)


## ------------------------------------------------------------
## case-specific vimp
## returns VIMP for each case
## ------------------------------------------------------------

o <- rfsrc(mpg~., mtcars)
op <- predict(o, importance = TRUE, csv = TRUE)
csvimp <- get.mv.csvimp(op, standardize=TRUE)
print(csvimp)

## ------------------------------------------------------------
## case-specific error rate
## returns tree-averaged error rate for each case
## ------------------------------------------------------------

o <- rfsrc(mpg~., mtcars)
op <- predict(o, importance = TRUE, cse = TRUE)
cserror <- get.mv.cserror(op, standardize=TRUE)
print(cserror)


## ------------------------------------------------------------
## predicted probability and predicted class labels are returned
## in the predict object for classification analyses
## ------------------------------------------------------------

data(breast, package = "randomForestSRC")
breast.obj <- rfsrc(status ~ ., data = breast[(1:100), ])
breast.pred <- predict(breast.obj, breast[-(1:100), ])
print(head(breast.pred$predicted))
print(breast.pred$class)


## ------------------------------------------------------------
## unique feature of randomForestSRC
## cross-validation can be used when factor labels differ over
## training and test data
## ------------------------------------------------------------

## first we convert all x-variables to factors
data(veteran, package = "randomForestSRC")
veteran2 <- data.frame(lapply(veteran, factor))
veteran2$time <- veteran$time
veteran2$status <- veteran$status

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

## train the forest and use this to predict on test data
o.grow <- rfsrc(Surv(time, status) ~ ., veteran2[train, ]) 
o.pred <- predict(o.grow, veteran2[-train , ])
print(o.grow)
print(o.pred)

## even harder ... factor level not previously encountered in training
veteran3 <- veteran2[1:3, ]
veteran3$celltype <- factor(c("newlevel", "1", "3"))
o2.pred <- predict(o.grow, veteran3)
print(o2.pred)
## the unusual level is treated like a missing value but is not removed
print(o2.pred$xvar)

## ------------------------------------------------------------
## example illustrating the flexibility of outcome = "test"
## illustrates restoration of forest via outcome = "test"
## ------------------------------------------------------------

## first we train the forest
data(pbc, package = "randomForestSRC")
pbc.grow <- rfsrc(Surv(days, status) ~ ., pbc)

## use predict with outcome = TEST
pbc.pred <- predict(pbc.grow, pbc, outcome = "test")

## notice that error rates are the same!!
print(pbc.grow)
print(pbc.pred)

## note this is equivalent to restoring the forest
pbc.pred2 <- predict(pbc.grow)
print(pbc.grow)
print(pbc.pred)
print(pbc.pred2)

## similar example, but with na.action = "na.impute"
airq.obj <- rfsrc(Ozone ~ ., data = airquality, na.action = "na.impute")
print(airq.obj)
print(predict(airq.obj))
## ... also equivalent to outcome="test" but na.action = "na.impute" required
print(predict(airq.obj, airquality, outcome = "test", na.action = "na.impute"))

## classification example
iris.obj <- rfsrc(Species ~., data = iris)
print(iris.obj)
print(predict.rfsrc(iris.obj, iris, outcome = "test"))

## ------------------------------------------------------------
## another example illustrating outcome = "test"
## unique way to check reproducibility of the forest
## ------------------------------------------------------------

## training step
set.seed(542899)
data(pbc, package = "randomForestSRC")
train <- sample(1:nrow(pbc), round(nrow(pbc) * 0.50))
pbc.out <- rfsrc(Surv(days, status) ~ .,  data=pbc[train, ])

## standard prediction call
pbc.train <- predict(pbc.out, pbc[-train, ], outcome = "train")
##non-standard predict call: overlays the test data on the grow forest
pbc.test <- predict(pbc.out, pbc[-train, ], outcome = "test")

## check forest reproducibilility by comparing "test" predicted survival
## curves to "train" predicted survival curves for the first 3 individuals
Time <- pbc.out$time.interest
matplot(Time, t(pbc.train$survival[1:3,]), ylab = "Survival", col = 1, type = "l")
matlines(Time, t(pbc.test$survival[1:3,]), col = 2)

## ------------------------------------------------------------
## multivariate forest example
## ------------------------------------------------------------

## train the forest
trn <- 1:20
o <- rfsrc(cbind(mpg, disp)~.,mtcars[trn,])

## print training results for each outcome
print(o, outcome.target="mpg")
print(o, outcome.target="disp")

## print test results for each outcome
p <- predict(o, mtcars[-trn,])
print(p, outcome.target="mpg")
print(p, outcome.target="disp")

Print Summary Output of a RF-SRC Analysis

Description

Print summary output from a RF-SRC analysis. This is the default print method for the package.

Usage

## S3 method for class 'rfsrc'
print(x, outcome.target = NULL, ...,
                          quantreg.tau = NULL)

Arguments

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

Examples

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

Quantile Regression Forests

Description

Estimates conditional quantiles for continuous responses using a forest prediction and the training-residual distribution, or the native Greenwald–Khanna quantile algorithm. Returns quantiles, CDF values, and distribution summaries for training or new observations. Univariate, multivariate, and mixed-outcome forests are supported.

Usage

## S3 method for class 'rfsrc'
quantreg(formula, data, object, newdata,
  method = "exchangeable", splitrule = NULL, prob = NULL, prob.epsilon = NULL,
  oob = TRUE, fast = FALSE, maxn = 1e3, ...)

extract.quantile(o)
get.quantile(o, target.prob = NULL, pretty = TRUE)
get.quantile.stat(o, pretty = TRUE)
get.quantile.crps(o, pretty = TRUE, subset = NULL, standardize = TRUE)
get.pinball.error(o, tau = NULL, subset = NULL, m.target = NULL,
                 pretty = TRUE)

Arguments

Details

Training and prediction

Supply formula and data to grow a forest, or object to reuse its saved forest. newdata can be supplied with either form of call. Training responses and residuals are saved once and reused by subsequent calls. Response values in newdata are used for evaluation; they do not replace this training information.

OOB residuals and their exchangeable distribution

With oob = TRUE during growth and OOB predictions available, the saved residual for training observation j is

R_j^{\mathrm{OOB}} = Y_j - \widehat m_j^{\mathrm{OOB}},

where \widehat m_j^{\mathrm{OOB}} averages predictions from trees whose training samples exclude observation j. Thus, each residual compares the observed response with a prediction from trees that did not use that observation to grow.

The exchangeable method pools these residuals into the empirical CDF

\widehat F_R^{\mathrm{OOB}}(u) = \frac{1}{N_R}\sum_{j\in\mathcal I_R} I\{R_j^{\mathrm{OOB}}\leq u\},

where \mathcal I_R indexes the finite residuals and N_R=|\mathcal I_R|. This is the OOB exchangeable residual distribution: each finite residual receives equal weight, and the same residual distribution is used at every predictor row. Exchangeability is the working assumption motivating this common distribution; OOB specifies how its residuals are constructed. Its lower-inverse quantile is denoted by \widehat Q_R^{\mathrm{OOB}}(\tau).

For training-row output using OOB mean predictions, the quantile is

\widehat Q_j^{\mathrm{OOB}}(\tau) = \widehat m_j^{\mathrm{OOB}} + \widehat Q_R^{\mathrm{OOB}}(\tau).

The common residual distribution includes each observation's finite OOB residual. For a new predictor row x, the quantile is

\widehat Q(\tau\mid x) = \widehat m(x) + \widehat Q_R^{\mathrm{OOB}}(\tau),

where \widehat m(x) is the full-ensemble prediction from the saved forest. The OOB residual distribution is retained even though the new-row mean prediction uses the full ensemble.

With oob = FALSE during growth, the selected residuals instead use full-ensemble training predictions. Both residual quantities are retained when available, but the original selection is reused in subsequent calls. Changing oob in an object call changes the current training-row output, not that residual selection. When growth and newdata are requested in the same call, oob first selects the training residual bank and then ordinary predictions are used for the new rows.

If an OOB output component is absent, the existing full-ensemble fallback is used. Missing values within an available OOB component remain missing; nonfinite residuals are excluded from the empirical distribution. Each residual-method quantile component records residual.source = "oob" or "ensemble" separately from its oob setting for the current predictions. A new-data result can therefore have residual.source = "oob" and oob = FALSE. Its numerical training OOB residuals remain in forest$residual.oob; the suffix .oob identifies the OOB values, not a logical flag.

The forest method uses the same saved residual bank but assigns row-specific forest weights. For training rows it requests OOB weights when oob = TRUE and full-ensemble weights otherwise. For GK, oob selects the native OOB or full-ensemble quantile output; GK quantiles are not calculated from residuals.

Response data and residuals in the returned object

The usual xvar and yvar components describe the observations for which predictions were returned. In a grow object they contain training data. In a new-data prediction object they contain the retained prediction data, with yvar available when responses were retained for evaluation. The saved forest's xvar and yvar continue to describe the training data. Quantile reporting grids and the null CRPS reference use these training responses, whereas performance scores use the current object's responses.

For a univariate grow object o, o$residual is o$yvar - o$predicted, and o$residual.oob is o$yvar - o$predicted.oob. For multivariate and mixed outcomes, these vectors are stored beside the corresponding predictions in o$regrOutput[[response.name]]. No residual is substituted for an unavailable OOB prediction.

The saved forest retains residual and residual.oob for use with future prediction rows. These are vectors for univariate regression and matrices for multivariate or mixed outcomes, with training observations in rows and named continuous responses in columns. An unavailable component is NULL in the grow output; its column in a partially available saved bank is NA. A wholly unavailable bank is NULL.

The named setting forest$quantreg$residual.source records the original residual selection for each continuous response. Residual-method quantile components also report their residual.source. Neither setting replaces a numerical residual.oob component. A prediction object retains the training residuals under forest; it does not place these training-length vectors beside its current predictions. Retain the original grow object for further quantreg calls.

Residual quantiles

For one continuous response, write R_j for its saved residuals, constructed as described above, and \widehat m(x) for the current forest mean prediction.

The default "exchangeable" method adds quantiles of the common empirical residual distribution to the current prediction:

\widehat Q(\tau\mid x) = \widehat m(x) + \widehat Q_R(\tau).

Its CDF is

\widehat F(s\mid x) = \frac{1}{N_R}\sum_{j\in\mathcal I_R} I\{R_j \leq s-\widehat m(x)\},

where \mathcal I_R indexes its N_R finite residuals. Each row uses the same residual-distribution shape, translated by its own forest prediction. This is the locally adjusted residual-CDF construction associated with Zhang et al. (2019).

The "forest" method weights the residuals using forest weights for the evaluated row:

\widehat F(s\mid x) = \sum_j w_j(x) I\{R_j \leq s-\widehat m(x)\}.

These weights allow the residual distribution to vary with the predictor row. They are applied to residuals, rather than directly to training responses as in Meinshausen (2006).

Both methods use the lower generalized inverse. After ordering residuals and their weights, the residual quantile is the first residual whose cumulative weight reaches \tau. The current forest prediction is then added. Equal weights give the inverse empirical-CDF convention, corresponding to a type-1 sample quantile. Quantiles can extend beyond the observed training-response range.

The mean and standard deviation are computed from the same residual distribution. With \bar R(x)=\sum_j w_j(x)R_j, they are

\widehat\mu(x)=\widehat m(x)+\bar R(x),\qquad \widehat\sigma(x)=\left\{\sum_j w_j(x) [R_j-\bar R(x)]^2\right\}^{1/2}.

Equal weights are used for "exchangeable". The residuals are not recentered to force their weighted mean to zero, so the distribution mean need not equal the forest prediction.

Nonfinite residuals and their associated weight columns are removed together. Remaining forest weights are normalized separately for each row. A row with a nonfinite current prediction, nonfinite retained weights, or no positive usable weight has unavailable quantiles and summaries. No usable residuals also yields unavailable output. Negative weights and incompatible weight dimensions are rejected.

Greenwald–Khanna

The "gk" method requests quantiles from the native Greenwald–Khanna algorithm. It avoids requesting the explicit forest-weight matrix. The returned quantiles are retained directly; a reporting CDF is then reconstructed from their values and probability levels. At each grid value this CDF uses the largest requested probability whose quantile is at or below that value, and reaches one at the upper training-response endpoint. Consequently, prob affects both native quantiles and the resolution of this reconstructed CDF. maxn only controls its evaluation grid.

CDF reporting grid

The grid, yunq, consists of sorted unique finite training responses, separately for each continuous response. When there are more than maxn values, an approximately equally spaced selection of their order positions is retained. maxn = Inf retains all values.

For the residual methods, cdf evaluates the residual-based CDF at these grid values. These evaluations do not determine the primary quantiles or moments. In particular, the last CDF value can be less than one when a shifted residual lies above the reporting grid; the primary quantiles and moments still use that residual.

density contains CDF increments. Its first column is the CDF at the first grid value, and later columns are successive differences. Each increment is probability in a reporting interval, rather than a density divided by grid spacing. Its row sum equals the last reported CDF value. These increments are not renormalized to represent omitted upper-tail mass.

Quantile storage is proportional to the number of evaluated rows times length(prob). CDF and increment storage are proportional to that row count times the grid size. Forest weights additionally require one column per training observation. Reducing maxn does not reduce the forest-weight matrix or the number of native GK quantiles requested.

Extracting results

extract.quantile(o) returns a named list with one component per returned continuous response, including univariate objects.

get.quantile(o) extracts stored quantiles. With target.prob = NULL, all stored levels are returned. Requested probabilities are sorted and deduplicated, then matched to the nearest stored levels. Column labels use the requested probabilities. Quantiles are not recomputed by this helper. To calculate additional levels, call quantreg on the original grow object with those levels in prob.

get.quantile.stat(o) returns columns mean, median, and std. For the residual methods, the mean and standard deviation use the direct residual summaries described above. The median retains the nearest-stored lookup at probability 0.5; include 0.5 in prob to obtain that exact level. For GK, the mean and standard deviation retain the reporting-grid moment calculation from the CDF increments, so these summaries depend on the grid and its represented probability mass.

Matrices retain their row and column dimensions when only one row, probability, or grid value is present.

CRPS

get.quantile.crps(o) returns a score curve with columns y and crps. Observed responses must be present.

At each grid value s_k, the helper averages [I\{Y_i\leq s_k\}-\widehat F(s_k\mid X_i)]^2 over the selected rows with finite response and CDF values. It accumulates the trapezoidal integral from the first grid value to s_k. With standardize = TRUE, each integral is divided by its integration width. Zero-width standardized entries are NA; the corresponding raw integrals are zero. The last entry summarizes the integral over the retained grid. Integration outside that grid is not included.

Pinball loss

get.pinball.error(o) returns the mean pinball loss at each requested level; smaller values are better. For observed response Y_i and predicted quantile \widehat Q_i(\tau), the loss is

\rho_\tau\{Y_i-\widehat Q_i(\tau)\},\qquad \rho_\tau(u)=u\{\tau-I(u<0)\}.

Averaging uses the selected rows with finite responses and quantiles, separately at each level; no usable rows gives NA.

For "exchangeable" and "forest", scoring uses exact stored levels or quantiles calculated from the saved residual distribution, independently of the reporting grid. Thus tau need not have been included in prob. GK retains interpolation of its reporting CDF. Scoring preserves the object's prediction mode and residual selection without refitting or calling prediction.

Printed performance

For univariate quantile output with observed responses and performance output, print.rfsrc reports raw and standardized grid-integrated CRPS and pinball losses. A selected continuous response from a multivariate object can be printed using print(object, outcome.target = "response.name"). To extract that response, use its name in the list returned by get.quantile(object, pretty = FALSE). These selections act on the returned results; restoration and prediction retain all responses. Without outcome.target, the multivariate print displays the existing performance summary, with the mean error followed by the response-specific errors, rather than response-specific quantile scores. Quantile-score OOB labels follow the quantile calculation for that response. The requested regression error retains its own provenance.

Pinball levels can be supplied for one print call using print(object, quantreg.tau = c(.2, .5, .8)), or set for the session using options(quantreg.tau = c(.2, .5, .8)). Explicit levels take precedence over the session option. When neither is supplied, the former rfsrc.pinball.taus option is accepted for compatibility, followed by the default c(0.1, 0.5, 0.9). These are reporting settings, not arguments to quantreg, and do not modify the fitted object. Levels must be finite and lie strictly between zero and one; they need not have been included in prob. Invalid levels give an explanatory error rather than silently suppressing the scores.

plot.quantreg(object, quantreg.tau = c(.2, .5, .8)) adds the same losses to the main plot as an annotation. This optional annotation leaves the plotted interval and the CRPS inset unchanged.

Value

An rfsrc grow or prediction object with the additional class quantreg. Existing classes are retained. For a univariate regression response, the component quantreg contains:

For multivariate and mixed-outcome forests, quantreg is a named list of these components for the returned continuous responses. Classification responses have no quantile component.

Grow objects also provide numerical residual and residual.oob components, at the top level for univariate regression and within each continuous response's regrOutput component otherwise. The saved forest retains the training residual vectors or named matrices and the residual-source settings described in Details. The usual xvar, yvar, and other forest output components retain their meanings.

Helper functions

extract.quantile always returns a response-named list of the quantile components described above. The other helpers return the following for one continuous response when pretty = TRUE, and a response-named list otherwise:

get.quantile

An observation-by-probability matrix, with column names such as q.50.

get.quantile.stat

A data frame with columns mean, median, and std, one row per observation.

get.quantile.crps

A data frame with columns y and crps, one row per reporting-grid value.

get.pinball.error

A numeric vector of mean losses, with names such as "tau=0.2", one entry per requested level.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

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

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

See Also

rfsrc, predict.rfsrc

Examples

## ------------------------------------------------------------
## A basic analysis using default settings
## ------------------------------------------------------------
o <- quantreg(Temp ~ ., data = na.omit(airquality))
plot.quantreg(o)


## ------------------------------------------------------------
## Exchangeable OOB residual quantiles for wine alcohol
## ------------------------------------------------------------
data(wine, package = "randomForestSRC")
set.seed(17)
prob <- c(.05, .25, .50, .75, .95)
o <- quantreg(alcohol ~ ., data = wine, method = "exchangeable",
              oob = TRUE, prob = prob, ntree = 100)
print(head(get.quantile(o)))
print(head(get.quantile(o, c(.25, .50, .75))))
print(head(get.quantile.stat(o)))

## The common residual quantiles are added to each OOB mean prediction.
r.oob <- as.numeric(o$residual.oob)
r.oob <- r.oob[is.finite(r.oob)]
q.res <- as.numeric(quantile(r.oob, probs = prob, type = 1))
q.oob <- outer(as.numeric(o$predicted.oob), q.res, "+")
print(all.equal(unname(get.quantile(o)), unname(q.oob)))
print(o$quantreg$residual.source)

## Inspect the stored levels and reporting-grid probability mass.
print(o$quantreg$prob)
print(summary(rowSums(o$quantreg$density)))
crps <- get.quantile.crps(o)
print(crps)
plot(crps, type = "l")

## Pinball losses at nondefault levels: no refitting is needed.
print(o, quantreg.tau = c(.2, .5, .8))
print(get.pinball.error(o, tau = c(.2, .5, .8)))
plot.quantreg(o, quantreg.tau = c(.2, .5, .8))

## Optionally use the same reporting levels throughout the session.
op <- options(quantreg.tau = c(.2, .5, .8))
print(getOption("quantreg.tau"))
print(o)
options(op)

## ------------------------------------------------------------
## Forest-weighted residuals and predictor-only deployment
## ------------------------------------------------------------
set.seed(23)
train <- sample.int(nrow(wine), floor(.7 * nrow(wine)))
o <- quantreg(alcohol ~ ., data = wine[train, ],
              method = "forest", prob = prob, ntree = 100)
o.test <- quantreg(object = o, newdata = wine[-train, ],
                   method = "forest")
print(head(get.quantile(o.test)))
print(tail(get.quantile.crps(o.test, standardize = FALSE), 1))

x.test <- wine[-train, setdiff(names(wine), "alcohol"), drop = FALSE]
o.predict <- quantreg(object = o, newdata = x.test, method = "forest")
print(head(get.quantile(o.predict)))

## Exchangeable prediction uses the same saved OOB residual bank.
o.exchangeable <- quantreg(object = o, newdata = x.test,
                            method = "exchangeable")
print(head(get.quantile(o.exchangeable)))
print(identical(o.exchangeable$forest$residual.oob, o$residual.oob))
print(o.exchangeable$quantreg$residual.source)  # "oob": training residuals
print(o.exchangeable$quantreg$oob)              # FALSE: new-row predictions
print(c(evaluated.rows = o.exchangeable$n,
        training.rows = length(o.exchangeable$forest$residual.oob)))

## Calculate a single level using the original grow object.
o.median <- quantreg(object = o, newdata = x.test,
                     method = "forest", prob = .5)
print(head(get.quantile(o.median)))

## ------------------------------------------------------------
## Grid resolution does not determine residual quantiles or moments
## ------------------------------------------------------------
coarse <- quantreg(object = o, newdata = x.test, method = "exchangeable",
                  prob = prob, maxn = 2, seed = 37)
fine <- quantreg(object = o, newdata = x.test, method = "exchangeable",
                prob = prob, maxn = Inf, seed = 37)
print(all.equal(get.quantile(coarse), get.quantile(fine)))
print(all.equal(get.quantile.stat(coarse), get.quantile.stat(fine)))

## ------------------------------------------------------------
## Multivariate and mixed outcomes
## ------------------------------------------------------------
dta <- na.omit(airquality)
mv <- quantreg(cbind(Ozone, Temp) ~ ., data = dta,
               splitrule = "mahalanobis", prob = prob, ntree = 100)
q.mv <- get.quantile(mv, pretty = FALSE)
print(names(q.mv))
print(head(q.mv$Ozone))
print(head(q.mv$Temp))
print(head(mv$regrOutput$Temp$residual.oob))
print(identical(mv$regrOutput$Temp$residual.oob,
                as.numeric(mv$forest$residual.oob[, "Temp"])))

## Restore all responses; select Temp only for extraction and printing.
mv.restore <- quantreg(object = mv, prob = .5)
q.restore <- get.quantile(mv.restore, pretty = FALSE)
print(names(q.restore))
print(head(q.restore$Temp))
print(get.mv.error(mv.restore))
print(mv.restore, outcome.target = "Temp")

dta$Month <- factor(dta$Month)
mixed <- quantreg(cbind(Ozone, Temp, Month) ~ ., data = dta,
                  prob = prob, ntree = 100)
print(names(extract.quantile(mixed)))

## ------------------------------------------------------------
## Native GK quantiles and alternative regression splitting
## ------------------------------------------------------------
gk <- quantreg(alcohol ~ ., data = wine, method = "gk",
               prob = c(.05, .25, .50, .75, .95),
               prob.epsilon = .01, maxn = 100, ntree = 100)
print(head(get.quantile(gk)))

mse <- quantreg(alcohol ~ ., data = wine, method = "exchangeable",
                splitrule = "mse", prob = c(.05, .50, .95), ntree = 100)
print(head(get.quantile(mse)))

## ------------------------------------------------------------
## Larger data set; iowa housing
## ------------------------------------------------------------

data(housing,  package = "randomForestSRC")

## the original data contains lots of missing data; use fast imputation
iowa <- housing
iowa$PID <- NULL
iowa$SalePrice <- log(iowa$SalePrice)
iowa <- impute(SalePrice ~. , iowa, splitrule = "random", nimpute = 1)

## use fewer trees and shallow trees for speed
o <- quantreg(SalePrice ~., iowa, ntree = 50, nodesize = 20)
plot.quantreg(o, prbL=.05, prbU=.95)

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

Description

Fast OpenMP-parallel implementation of random forests (Breiman, 2001) for regression, classification, survival analysis (Ishwaran et al., 2008), competing risks (Ishwaran et al., 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).

The package supports both deterministic and randomized splitting rules (Geurts et al., 2006; Ishwaran, 2015) across all families. Multiple types of variable importance (VIMP) are available, including holdout VIMP and confidence regions (Ishwaran and Lu, 2019), for both individual and grouped variables. Variable selection can be performed using minimal depth (Ishwaran et al., 2010, 2011). Fast interfaces for missing data imputation are provided using several forest-based algorithms (Tang and Ishwaran, 2017).

Highlighted updates:

1. For survival and competing risk analysis, concordance-based performance is now computed using Uno inverse-probability-of-censoring weighting (Uno et al. 2011). This affects all survival performance values derived from the concordance index, including out-of-bag and test error rates and variable importance (VIMP). To revert to unweighted Harrell concordance, set use.uno = FALSE when fitting a survival forest. Concordance calculations now use an efficient O(n \log n) algorithm based on a binary indexed tree (Fenwick 1994, Therneau 2024), replacing the naive O(n^2) pairwise computation for large n.

2. For variable selection, we recommend using VarPro, an R package for model-independent variable selection using rule-based variable priority. It supports regression, classification, survival analysis, and includes a new mode for unsupervised learning. See https://www.varprotools.org for more information.

3. For computational speed, the default VIMP method has changed from "permute" (Breiman-Cutler permutation) to "anti" (importance = "anti" or importance = TRUE). While faster, this may be less accurate in settings such as highly imbalanced classification. To revert to permutation VIMP, use importance = "permute".

This is the main entry point to the randomForestSRC package. For more information on OpenMP support and the package as a whole, see package?randomForestSRC.

Usage

rfsrc(formula, data, ntree = 500,
  mtry = NULL, ytry = NULL,
  nodesize = NULL, nodedepth = NULL,
  splitrule = NULL, nsplit = NULL,
  importance = c(FALSE, TRUE, "none", "anti", "permute", "random"),
  block.size = if (any(is.element(as.character(importance),
                     c("none", "FALSE")))) NULL else 10,
  bootstrap = c("by.root", "none", "by.user"),
  samptype = c("swor", "swr"), samp = NULL, membership = FALSE,
  sampsize = if (samptype == "swor") function(x){x * .632} else function(x){x},
  na.action = c("na.omit", "na.impute"), nimpute = 1,
  ntime = 150, cause,
  perf.type = NULL,
  proximity = FALSE, distance = FALSE, forest.wt = FALSE,
  xvar.wt = NULL, yvar.wt = NULL, split.wt = NULL, case.wt  = NULL,
  case.depth = FALSE,
  forest = TRUE,
  use.uno = TRUE, save.memory = FALSE,
  var.used = c(FALSE, "all.trees", "by.tree"),
  split.depth = c(FALSE, "all.trees", "by.tree"),
  seed = NULL,
  do.trace = FALSE,
  ...)

## convenient interface for growing a CART tree
rfsrc.cart(formula, data, ntree = 1, mtry = ncol(data),
   bootstrap = "none", nsplit = 0, ...)

Arguments

Details

1. Types of forests

   The type of forest is automatically inferred from the outcome and formula. Supported forest types include:

   • Regression forests for continuous outcomes.

   • Classification forests for factor outcomes.

   • Multivariate forests for continuous, categorical, or mixed outcomes.

   • Unsupervised forests when no outcome is specified.

   • Survival forests for right-censored time-to-event data.

   • Competing risk forests for multi-event survival settings.

2. Splitting

   1. Splitting rules are set using the splitrule option.

   2. Random splitting is invoked via splitrule = "random".

   3. Use nsplit to enable randomized splitting and improve speed; see Improving computational speed.

3. Available splitting rules

   • Regression

     1. "mse" (default): weighted mean squared error (Breiman et al., 1984).

     2. "quantile.regr": quantile regression via check-loss; see quantreg.rfsrc for quantile estimation methods and extraction helpers.

     3. "la.quantile.regr": local adaptive quantile regression.

   • Classification

     1. "gini" (default): Gini index.

     2. "auc": AUC-based splitting; appropriate for imbalanced data.

     3. "entropy": entropy-based splitting.

   • Survival

     1. "logrank" (default): log-rank splitting.

     2. "bs.gradient": Brier score gradient splitting. Uses 90th percentile of observed times by default, or set prob.

     3. "logrankscore": log-rank score splitting.

   • Competing risks (see Ishwaran et al., 2014)

     1. "logrankCR" (default): Gray's test-based weighted log-rank splitting.

     2. "logrank": cause-specific weighted log-rank; use cause to target specific events.

   • Multivariate

     1. Default: normalized composite splitting (Tang and Ishwaran, 2017).

     2. "mahalanobis": Mahalanobis splitting with optional covariance matrix; for multivariate regression.

   • Unsupervised Splitting uses pseudo-outcomes and the composite rule. See sidClustering for advanced unsupervised analysis.

   • Custom splitting Custom rules can be defined using splitCustom.c. Up to 16 rules per family are allowed. Use "custom", "custom1", etc. Compilation required.

4. Improving computational speed

   See rfsrc.fast. Strategies include:

   • Increase nodesize.

   • Set save.memory = TRUE for large survival or competing risk models.

   • Set block.size = NULL to avoid repeated cumulative error computation.

   • Use perf.type = "none" to disable VIMP and C-index calculations.

   • Set nsplit to a small integer (e.g., 1-10).

   • Reduce bootstrap size with sampsize, samptype.

   • Set ntime to a coarse grid (e.g., 50) for survival models.

   • Pre-filter variables; use max.subtree for fast variable selection.

5. Prediction Error

   Error is computed using OOB data:

   • Regression: mean squared error.

   • Classification: misclassification rate, Brier score, or G-mean-based error, as selected by perf.type.

   • Survival: C-error = 1 - concordance, using Uno weighting when use.uno = TRUE and unweighted Harrell concordance otherwise.

   If bootstrap = "none", OOB error is unavailable. Use predict.rfsrc for cross-validation error instead.

   See Classification Performance Metrics (classification.performance) for additional classification scores and get.imbalanced.performance for binary-class performance summaries and threshold choices. For right-censored survival, plotBrierAUC displays Brier and AUC curves, including a Kaplan–Meier reference for Brier score; get.cindex calculates concordance-based error.

6. Variable Importance (VIMP)

   VIMP methods:

   • "permute": permutation VIMP (Breiman-Cutler).

   • "random": randomized left/right assignment.

   • "anti" (default): anti-split assignment.

   The block.size option controls granularity. For confidence intervals and the extract.subsample and extract.bootsample helpers, see subsample.rfsrc; use plot.subsample.rfsrc to display the results. Also see holdout.vimp for a more conservative variant.

7. Multivariate Forests

   Use:

   rfsrc(Multivar(y1, ..., yd) ~ ., data)

   or

   rfsrc(cbind(y1, ..., yd) ~ ., data)

   Use get.mv.formula to construct a multivariate formula. For get.mv.predicted, get.mv.error, get.mv.vimp, block errors, and case-specific extraction helpers, see Extracting Multivariate Values (multivariate.values).

8. Unsupervised Forests

   Use:

   rfsrc(data = X)

   or

   rfsrc(Unsupervised(ytry) ~ ., data = X)

   Random subsets of ytry pseudo-responses are used for each mtry variable. No performance metrics are computed.

9. Survival, Competing Risks

   • Survival: use Surv(time, status) ~ .. Status must be 0 (censored) or 1 (event).

   • Competing risks: status = 0 (censored), 1-J (event types). Use cause to target specific events.

   • Larger nodesize is typically needed for competing risks.

10. Missing data imputation

    Use na.action = "na.impute". Iteration with nimpute > 1 replaces missing values using OOB predictions. Observations or variables with all missing values are removed. See impute.rfsrc for standalone imputation and impute.learn.rfsrc for learning an imputer to apply to new data.

11. Allowable data types and factors

    Variables must be numeric, integer, factor, or logical. Non-factors are coerced to numeric. For unordered factors, all complementary subsets are considered for splits.

    Factor levels are mapped to ensure consistency across training/test data. Consider converting factors to numeric for high-dimensional settings.

Value

An object of class (rfsrc, grow) with the following components:

call

The original call to rfsrc.

family

The family used in the analysis.

n

Sample size after applying na.action.

ntree

Number of trees grown.

mtry

Number of variables randomly selected at each node.

nodesize

Minimum terminal node size.

nodedepth

Maximum depth allowed for each tree.

splitrule

Splitting rule used.

nsplit

Number of random split points.

yvar

Response values.

yvar.names

Character vector of response variable names.

xvar

Data frame of predictor variables.

xvar.names

Character vector of predictor variable names.

xvar.wt

Non-negative weights specifying the selection probability of each variable.

split.wt

Non-negative weights adjusting each variable's split statistic.

cause.wt

Weights for composite competing risk splitting.

leaf.count

Number of terminal nodes per tree. A value of 0 indicates a rejected tree (may occur with missing data); a value of 1 indicates a stump.

proximity

Proximity matrix indicating how often case pairs fall in the same terminal node.

forest

Forest object, returned if forest=TRUE. Required for prediction and most wrappers.

forest.wt

Forest weight matrix.

membership

Terminal node membership matrix (rows: cases; columns: trees).

inbag

Inbag count matrix (rows: cases; columns: trees).

var.used

Number of times each variable is used to split a node.

imputed.indv

Indices of individuals with missing values.

imputed.data

Imputed dataset with responses followed by predictors.

split.depth

Matrix or array recording minimal split depth of variables by case and tree.

err.rate

Cumulative OOB error rate.

err.block.rate

Cumulative error per ensemble block (size defined by block.size). If block.size = 1, error per tree.

importance

Variable importance (VIMP) for each predictor.

predicted

In-bag predicted values.

predicted.oob

Out-of-bag (OOB) predicted values.

class

(Classification) In-bag predicted class labels.

class.oob

(Classification) OOB predicted class labels.

regrOutput

(Multivariate) Response-named list of predictions and performance results for continuous outcomes.

classOutput

(Multivariate) Response-named list of predictions and performance results for categorical outcomes.

survival

(Survival) In-bag survival functions.

survival.oob

(Survival) OOB survival functions.

chf

(Survival or competing risks) In-bag cumulative hazard function.

chf.oob

(Survival or competing risks) OOB cumulative hazard function.

time.interest

(Survival or competing risks) Unique sorted event times.

ndead

(Survival or competing risks) Total number of observed events.

cif

(Competing risks) In-bag cumulative incidence function by cause.

cif.oob

(Competing risks) OOB cumulative incidence function by cause.

Note

Values returned by the forest depend on the family:

• Regression: predicted and predicted.oob are vectors of predicted values.

• Classification: predicted and predicted.oob are matrices of class probabilities. VIMP and performance metrics are returned as a matrix with J+1 columns (J = number of classes). The first column ("all") gives unconditional results; remaining columns give class-conditional results.

• Survival: predicted contains mortality estimates (Ishwaran et al., 2008). These are calibrated to the number of expected events under identical covariate profiles. Also returned are matrices of the survival function and CHF for each individual over time.interest.

• Competing risks: predicted contains expected life years lost by cause (Ishwaran et al., 2013). Also returned are three-dimensional arrays for CIF and CSCHF indexed by case, time, and event type.

• Multivariate: Predictions, VIMP, and error rates are returned in regrOutput and classOutput. Use get.mv.predicted, get.mv.vimp, and get.mv.error to extract results. Their OOB selection, standardization, and return formats are described in Extracting Multivariate Values (multivariate.values).

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Breiman L., Friedman J.H., Olshen R.A. and Stone C.J. (1984). Classification and Regression Trees, Belmont, California.

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

Cutler A. and Zhao G. (2001). PERT-Perfect random tree ensembles. Comp. Sci. Statist., 33: 490-497.

Dietterich, T. G. (2000). An experimental comparison of three methods for constructing ensembles of decision trees: bagging, boosting, and randomization. Machine Learning, 40, 139-157.

Fenwick, P.M. (1994). A new data structure for cumulative frequency tables. Software: Practice and Experience, 24(3):327–336.

Gray R.J. (1988). A class of k-sample tests for comparing the cumulative incidence of a competing risk, Ann. Statist., 16: 1141-1154.

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.

Harrell et al. F.E. (1982). Evaluating the yield of medical tests, J. Amer. Med. Assoc., 247:2543-2546.

Hothorn T. and Lausen B. (2003). On the exact distribution of maximally selected rank statistics, Comp. Statist. Data Anal., 43:121-137.

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

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

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.

Lin, Y. and Jeon, Y. (2006). Random forests and adaptive nearest neighbors. J. Amer. Statist. Assoc., 101(474), 578-590.

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

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.

LeBlanc M. and Crowley J. (1993). Survival trees by goodness of split, J. Amer. Statist. Assoc., 88:457-467.

Loh W.-Y and Shih Y.-S (1997). Split selection methods for classification trees, Statist. Sinica, 7:815-840.

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.

Mogensen, U.B, Ishwaran H. and Gerds T.A. (2012). Evaluating random forests for survival analysis using prediction error curves, J. Statist. Software, 50(11): 1-23.

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

Segal M.R. (1988). Regression trees for censored data, Biometrics, 44:35-47.

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.

Therneau, T.M. (2024). Concordance. Vignette in the survival R package.

Uno, H., Cai, T., Pencina, M.J., D'Agostino, Ralph B. and Wei, L-J. (2011). On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data. Statistics in Medicine, 30:1105-1117.

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

See Also

Classification Performance Metrics (classification.performance), Extracting Multivariate Values (multivariate.values), Fast Saving and Loading Objects (fast.saveload),

get.tree.rfsrc,

holdout.vimp,

imbalanced.rfsrc, impute.rfsrc, impute.learn.rfsrc,

max.subtree.rfsrc,

partial.rfsrc, plot.competing.risk.rfsrc, plot.rfsrc, plot.quantreg.rfsrc, plot.subsample.rfsrc, plot.survival.rfsrc, plotBrierAUC, plot.variable.rfsrc, predict.rfsrc, print.rfsrc,

quantreg.rfsrc,

rfsrc, rfsrc.anonymous, rfsrc.cart, rfsrc.fast,

sidClustering.rfsrc,

subsample.rfsrc,

tune.rfsrc,

vimp.rfsrc

Examples

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

## veteran data
## randomized trial of two treatment regimens for lung cancer
data(veteran, package = "randomForestSRC")
v.obj <- rfsrc(Surv(time, status) ~ ., data = veteran, block.size = 1)

## plot tree number 3
plot(get.tree(v.obj, 3))

## print results of trained forest
print(v.obj)

## plot results of trained forest
plot(v.obj)

## plot survival curves for first 10 individuals -- direct way
matplot(v.obj$time.interest, 100 * t(v.obj$survival.oob[1:10, ]),
    xlab = "Time", ylab = "Survival", type = "l", lty = 1)

## plot survival curves for first 10 individuals
## using function "plot.survival" 
plot.survival(v.obj, subset = 1:10)

## Brier and AUCT helpers with pointwise intervals
## see plot.survival for more details
brier.obj <- get.brier.survival(v.obj, conf.int = TRUE)
head(brier.obj$brier.score)
auct.obj <- get.auct.survival(v.obj, conf.int = TRUE)
head(auct.obj$auct.score)

## Brier and AUCT curves with confidence bands
## see plot.survival for more details
plotBrierAUC(v.obj, xlim = c(0, 600))
plotBrierAUC(v.obj, plots = "brier", xlim = c(0, 600))
plotBrierAUC(v.obj, plots = "auct", xlim = c(0, 600))

## plot CRPS (continuous rank probability score) as function of time
## here's how to calculate the CRPS for every time point
trapz <- randomForestSRC:::trapz
time <- v.obj$time.interest
bscore <- brier.obj$brier.score
crps <- sapply(1:length(time), function(j) {
  trapz(time[1:j], bscore[1:j, 2] / diff(range(time[1:j])))
})
plot(time, crps, ylab = "CRPS", type = "s", lwd = 2)


## fast nodesize optimization for veteran data
## optimal nodesize in survival is larger than other families
## see the function "tune" for more examples
tune.nodesize(Surv(time,status) ~ ., veteran)


## Primary biliary cirrhosis (PBC) of the liver
data(pbc, package = "randomForestSRC")
pbc.obj <- rfsrc(Surv(days, status) ~ ., pbc)
print(pbc.obj)


## save.memory example for survival
## growing many deep trees creates memory issue without this option!
data(pbc, package = "randomForestSRC")
print(rfsrc(Surv(days, status) ~ ., pbc, splitrule = "random",
            ntree = 25000, nodesize = 1, save.memory = TRUE))



##------------------------------------------------------------
## trees can be plotted for any family 
## see get.tree for details and more examples
##------------------------------------------------------------

## survival where factors have 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))

## classification
iris.obj <- rfsrc(Species ~., data = iris)
plot(get.tree(iris.obj, 25, class.type = "bayes"))
plot(get.tree(iris.obj, 25, target = "setosa"))
plot(get.tree(iris.obj, 25, target = "versicolor"))
plot(get.tree(iris.obj, 25, target = "virginica"))

## ------------------------------------------------------------
## simple example of VIMP using iris classification
## ------------------------------------------------------------

## directly from trained forest
print(rfsrc(Species~.,iris,importance=TRUE)$importance)

## VIMP (and performance) use misclassification error by default
## but brier prediction error can be requested
print(rfsrc(Species~.,iris,importance=TRUE,perf.type="brier")$importance)

## example using vimp function (see vimp help file for details)
iris.obj <- rfsrc(Species ~., data = iris)
print(vimp(iris.obj)$importance)
print(vimp(iris.obj,perf.type="brier")$importance)

## example using hold out vimp (see holdout.vimp help file for details)
print(holdout.vimp(Species~.,iris)$importance)
print(holdout.vimp(Species~.,iris,perf.type="brier")$importance)

## ------------------------------------------------------------
## confidence interval for vimp using subsampling
## compare with holdout vimp
## ------------------------------------------------------------

## new York air quality measurements
o <- rfsrc(Ozone ~ ., data = airquality)
so <- subsample(o)
plot(so)

## compare with holdout vimp
print(holdout.vimp(Ozone ~ ., data = airquality)$importance)


##------------------------------------------------------------
## example of imputation in survival analysis
##------------------------------------------------------------

data(pbc, package = "randomForestSRC")
pbc.obj2 <- rfsrc(Surv(days, status) ~ ., pbc, na.action = "na.impute")

## same as above but iterate the missing data algorithm
pbc.obj3 <- rfsrc(Surv(days, status) ~ ., pbc,
         na.action = "na.impute", nimpute = 3)

## fast way to impute data (no inference is done)
## see impute for more details
pbc.imp <- impute(Surv(days, status) ~ ., pbc, splitrule = "random")

##------------------------------------------------------------
## compare RF-SRC to Cox regression
## Illustrates C-error and Brier score measures of performance
## assumes "pec" and "survival" libraries are loaded
##------------------------------------------------------------

if (library("survival", logical.return = TRUE)
    & library("pec", logical.return = TRUE)
    & library("prodlim", logical.return = TRUE))

{
  ##prediction function required for pec
  predictSurvProb.rfsrc <- function(object, newdata, times, ...){
    ptemp <- predict(object,newdata=newdata,...)$survival
    pos <- sindex(jump.times = object$time.interest, eval.times = times)
    p <- cbind(1,ptemp)[, pos + 1]
    if (NROW(p) != NROW(newdata) || NCOL(p) != length(times))
      stop("Prediction failed")
    p
  }

  ## data, formula specifications
  data(pbc, package = "randomForestSRC")
  pbc.na <- na.omit(pbc)  ##remove NA's
  surv.f <- as.formula(Surv(days, status) ~ .)
  pec.f <- as.formula(Hist(days,status) ~ 1)

  ## run cox/rfsrc models
  ## for illustration we use a small number of trees
  cox.obj <- coxph(surv.f, data = pbc.na, x = TRUE)
  rfsrc.obj <- rfsrc(surv.f, pbc.na, ntree = 150)

  ## compute bootstrap cross-validation estimate of expected Brier score
  ## see Mogensen, Ishwaran and Gerds (2012) Journal of Statistical Software
  set.seed(17743)
  prederror.pbc <- pec(list(cox.obj,rfsrc.obj), data = pbc.na, formula = pec.f,
                        splitMethod = "bootcv", B = 50)
  print(prederror.pbc)
  plot(prederror.pbc)

  ## compute out-of-bag C-error for cox regression and compare to rfsrc
  ## use uno weights throughout
  rfsrc.obj <- rfsrc(surv.f, pbc.na)
  uno.wt <- rfsrc.obj$forest$uno.weights$weight
  cat("out-of-bag Cox Analysis ...", "\n")
  cox.err <- sapply(1:100, function(b) {
    if (b%%10 == 0) cat("cox bootstrap:", b, "\n")
    train <- sample(1:nrow(pbc.na), nrow(pbc.na), replace = TRUE)
    cox.obj <- tryCatch({coxph(surv.f, pbc.na[train, ])}, error=function(ex){NULL})
    if (!is.null(cox.obj)) {
      get.cindex(pbc.na$days[-train],
                 pbc.na$status[-train],
                 predict(cox.obj, pbc.na[-train, ]),
                 weight=uno.wt[-train])
    } else NA
  })
  cat("\n\tOOB error rates (with Uno-IPCW):\n")
  cat("\tRSF            : ", get.mv.error(rfsrc.obj), "\n")
  cat("\tCox regression : ", mean(cox.err, na.rm = TRUE), "\n")
}

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

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

data(wihs, package = "randomForestSRC")
wihs.obj <- rfsrc(Surv(time, status) ~ ., wihs, nsplit = 3, ntree = 100)
plot.competing.risk(wihs.obj)
cif <- wihs.obj$cif.oob
Time <- wihs.obj$time.interest
idu <- wihs$idu
cif.haart <- cbind(apply(cif[,,1][idu == 0,], 2, mean),
                   apply(cif[,,1][idu == 1,], 2, mean))
cif.aids  <- cbind(apply(cif[,,2][idu == 0,], 2, mean),
                   apply(cif[,,2][idu == 1,], 2, mean))
matplot(Time, cbind(cif.haart, cif.aids), type = "l",
        lty = c(1,2,1,2), col = c(4, 4, 2, 2), lwd = 3,
        ylab = "Cumulative Incidence")
legend("topleft",
       legend = c("HAART (Non-IDU)", "HAART (IDU)", "AIDS (Non-IDU)", "AIDS (IDU)"),
       lty = c(1,2,1,2), col = c(4, 4, 2, 2), lwd = 3, cex = 1.5)


## illustrates the various splitting rules
## illustrates event specific and non-event specific variable selection
if (library("survival", logical.return = TRUE)) {

  ## use the pbc data from the survival package
  ## events are transplant (1) and death (2)
  data(pbc, package = "survival")
  pbc$id <- NULL

  ## modified Gray's weighted log-rank splitting
  ## (equivalent to cause=c(1,1) and splitrule="logrankCR")
  pbc.cr <- rfsrc(Surv(time, status) ~ ., pbc)
 
  ## log-rank cause-1 specific splitting and targeted VIMP for cause 1
  pbc.log1 <- rfsrc(Surv(time, status) ~ ., pbc, 
              splitrule = "logrankCR", cause = c(1,0), importance = TRUE)

  ## log-rank cause-2 specific splitting and targeted VIMP for cause 2
  pbc.log2 <- rfsrc(Surv(time, status) ~ ., pbc, 
              splitrule = "logrankCR", cause = c(0,1), importance = TRUE)

  ## extract VIMP from the log-rank forests: event-specific
  ## extract minimal depth from the Gray log-rank forest: non-event specific
  var.perf <- data.frame(md = max.subtree(pbc.cr)$order[, 1],
                         vimp1 = 100 * pbc.log1$importance[ ,1],
                         vimp2 = 100 * pbc.log2$importance[ ,2])
  print(var.perf[order(var.perf$md), ], digits = 2)

}

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

## new York air quality measurements
airq.obj <- rfsrc(Ozone ~ ., data = airquality, na.action = "na.impute")

# partial plot of variables (see plot.variable for more details)
plot.variable(airq.obj, partial = TRUE, smooth.lines = TRUE)

## motor trend cars
mtcars.obj <- rfsrc(mpg ~ ., data = mtcars)

## ------------------------------------------------------------
## regression with custom bootstrap
## ------------------------------------------------------------

ntree <- 25
n <- nrow(mtcars)
s.size <- n / 2
swr <- TRUE
samp <- randomForestSRC:::make.sample(ntree, n, s.size, swr)
o <- rfsrc(mpg ~ ., mtcars, bootstrap = "by.user", samp = samp)

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

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

## wisconsin prognostic breast cancer data
data(breast, package = "randomForestSRC")
breast.obj <- rfsrc(status ~ ., data = breast, block.size=1)
plot(breast.obj)

## ------------------------------------------------------------
## big data set, reduce number of variables using simple method
## ------------------------------------------------------------

## use Iowa housing data set
data(housing, package = "randomForestSRC")

## original data contains lots of missing data, use fast imputation
## however see impute for other methods
housing2 <- impute(data = housing, fast = TRUE)

## run shallow trees to find variables that split any tree
xvar.used <- rfsrc(SalePrice ~., housing2, ntree = 250, nodedepth = 4,
                   var.used="all.trees", mtry = Inf, nsplit = 100)$var.used

## now fit forest using filtered variables
xvar.keep  <- names(xvar.used)[xvar.used >= 1]
o <- rfsrc(SalePrice~., housing2[, c("SalePrice", xvar.keep)])
print(o)

## ------------------------------------------------------------
## imbalanced classification data
## see the "imbalanced" function for further details
##
## (a) use balanced random forests with undersampling of the majority class
## Specifically let n0, n1 be sample sizes for majority, minority
## cases.  We sample 2 x n1 cases with majority, minority cases chosen
## with probabilities n1/n, n0/n where n=n0+n1
##
## (b) balanced random forests using "imbalanced"
##
## (c) q-classifier (RFQ) using "imbalanced" 
##
## ------------------------------------------------------------

## Wisconsin breast cancer example
data(breast, package = "randomForestSRC")
breast <- na.omit(breast)

## balanced random forests - brute force
y <- breast$status
obdirect <- rfsrc(status ~ ., data = breast, nsplit = 10,
            case.wt = randomForestSRC:::make.wt(y),
            sampsize = randomForestSRC:::make.size(y))
print(obdirect)
print(get.imbalanced.performance(obdirect))

## balanced random forests - using "imbalanced" 
ob <- imbalanced(status ~ ., data = breast, method = "brf")
print(ob)
print(get.imbalanced.performance(ob))

## q-classifier (RFQ) - using "imbalanced" 
oq <- imbalanced(status ~ ., data = breast)
print(oq)
print(get.imbalanced.performance(oq))

## q-classifier (RFQ) - with auc splitting
oqauc <- imbalanced(status ~ ., data = breast, splitrule = "auc")
print(oqauc)
print(get.imbalanced.performance(oqauc))

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

## two equivalent ways to implement unsupervised forests
mtcars.unspv <- rfsrc(Unsupervised() ~., data = mtcars)
mtcars2.unspv <- rfsrc(data = mtcars)


## illustration of sidClustering for the mtcars data
## see sidClustering for more details
mtcars.sid <- sidClustering(mtcars, k = 1:10)
print(split(mtcars, mtcars.sid$cl[, 3]))
print(split(mtcars, mtcars.sid$cl[, 10]))


## ------------------------------------------------------------
## bivariate regression using Mahalanobis splitting
## also illustrates user specified covariance matrix
## ------------------------------------------------------------

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

  ## load boston housing data, specify the bivariate regression
  data(BostonHousing)
  f <- formula("Multivar(lstat, nox) ~.")
  
  ## Mahalanobis splitting  
  bh.mreg <- rfsrc(f, BostonHousing, importance = TRUE, splitrule = "mahal")
  
  ## performance error and vimp
  vmp <- get.mv.vimp(bh.mreg)
  pred <- get.mv.predicted(bh.mreg)

  ## standardized error and vimp
  err.std <- get.mv.error(bh.mreg, standardize = TRUE)
  vmp.std <- get.mv.vimp(bh.mreg, standardize = TRUE)

  ## same analysis, but with user specified covariance matrix
  sigma <- cov(BostonHousing[, c("lstat","nox")])
  bh.mreg2 <- rfsrc(f, BostonHousing, splitrule = "mahal", sigma = sigma)
  
}

## ------------------------------------------------------------
## multivariate mixed forests (nutrigenomic study)
## study effects of diet, lipids and gene expression for mice
## diet, genotype and lipids used as the multivariate y
## genes used for the x features
## ------------------------------------------------------------

## load the data (data is a list)
data(nutrigenomic, package = "randomForestSRC")

## assemble the multivariate y data
ydta <- data.frame(diet = nutrigenomic$diet,
                   genotype = nutrigenomic$genotype,
                   nutrigenomic$lipids)

## multivariate mixed forest call
## uses "get.mv.formula" for conveniently setting formula
mv.obj <- rfsrc(get.mv.formula(colnames(ydta)),
                data.frame(ydta, nutrigenomic$genes),
		importance=TRUE, nsplit = 10)

## print results for diet and genotype y values	    
print(mv.obj, outcome.target = "diet")
print(mv.obj, outcome.target = "genotype")

## extract standardized VIMP
svimp <- get.mv.vimp(mv.obj, standardize = TRUE)

## plot standardized VIMP for diet, genotype and lipid for each gene
boxplot(t(svimp), col = "bisque", cex.axis = .7, las = 2,
        outline = FALSE,
        ylab = "standardized VIMP",
        main = "diet/genotype/lipid VIMP for each gene")
	

## ------------------------------------------------------------
## illustrates yvar.wt which sets the probability of selecting 
## the response variables in multivariate regression
## ------------------------------------------------------------

## use mtcars: add fake responses
mult.mtcars  <- cbind(mtcars, mtcars$mpg, mtcars$mpg)
names(mult.mtcars) = c(names(mtcars), "mpg2", "mpg3")

## noise up the fake responses
mult.mtcars$mpg2 <- sample(mtcars$mpg)
mult.mtcars$mpg3 <- sample(mtcars$mpg)

formula = as.formula(Multivar(mpg, mpg2, mpg3)  ~ .)

## select 2 of the 3 responses randomly at each split with an associated weight vector.
## choose the noisy y responses which should degrade performance
yvar.wt = c(0.000001, 0.5, 0.5)
ytry = 2

mult.grow <- rfsrc(formula = formula, data = mult.mtcars, ytry = ytry, yvar.wt = yvar.wt)

print(mult.grow)
print(get.mv.error(mult.grow))

## Also, compare the following two results, as they should be similar:
yvar.wt = c(1.0, 00000.1, 00000.1)
ytry = 1

result1 = rfsrc(formula = formula, data = mult.mtcars, ytry = ytry, yvar.wt = yvar.wt)
result2 = rfsrc(mpg ~ ., mtcars)

print(get.mv.error(result1))
print(get.mv.error(result2))

## ------------------------------------------------------------
## custom splitting using the pre-coded examples
## ------------------------------------------------------------

## motor trend cars
mtcars.obj <- rfsrc(mpg ~ ., data = mtcars, splitrule = "custom")

## iris analysis
iris.obj <- rfsrc(Species ~., data = iris, splitrule = "custom1")

## WIHS analysis
wihs.obj <- rfsrc(Surv(time, status) ~ ., wihs, nsplit = 3,
                  ntree = 100, splitrule = "custom1")

Anonymous Random Forests

Description

Anonymous random forests is carefully modified to ensure that the original training data is not retained. This enables users to share the trained forest with others without disclosing the underlying data.

Usage

rfsrc.anonymous(formula, data, forest = TRUE, ...)

Arguments

Details

This function calls rfsrc and returns a forest object with the original training data removed. This enables users to share their forest while preserving the privacy of their data.

To enable prediction on new (test) data, certain minimal information from the training data must still be retained. This includes:

• Names of the original variables.

• For factor variables, the levels of each factor.

• Summary statistics used for imputation: the mean for continuous variables and the most frequent class for factors.

• Tree topology, including split points used to grow the trees.

For maximal privacy, users are strongly encouraged to replace variable names with non-identifiable labels and convert all variables to continuous format when possible. If factor variables are used, their levels should also be anonymized. However, the user is solely responsible for de-identifying the data and verifying that privacy is maintained. We provide NO GUARANTEES regarding data confidentiality.

Missing data handling: Anonymous forests do not support imputation of training data. The option na.action = "na.impute" is automatically downgraded to "na.omit". If training data contain missing values, we recommend pre-imputing them using impute.

Test data, however, can be imputed at prediction time:

• na.action = "na.impute" performs a fast imputation by replacing missing values with the training mean (for numeric variables) or most frequent class (for factors).

• na.action = "na.random" uses a fast random draw from training distributions for imputation.

Although anonymous forests are compatible with many package functions, they are only guaranteed to work with functions that do not explicitly require access to the original training data.

Value

An object of class (rfsrc, grow, anonymous).

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

See Also

rfsrc

Examples

## ------------------------------------------------------------
## regression
## ------------------------------------------------------------
print(rfsrc.anonymous(mpg ~ ., mtcars))

## ------------------------------------------------------------
## plot anonymous regression tree (using get.tree)
## TBD CURRENTLY NOT IMPLEMENTED 
## ------------------------------------------------------------
## plot(get.tree(rfsrc.anonymous(mpg ~ ., mtcars), 10))

## ------------------------------------------------------------
## classification
## ------------------------------------------------------------
print(rfsrc.anonymous(Species ~ ., iris))

## ------------------------------------------------------------
## survival
## ------------------------------------------------------------
data(veteran, package = "randomForestSRC")
print(rfsrc.anonymous(Surv(time, status) ~ ., data = veteran))

## ------------------------------------------------------------
## competing risks
## ------------------------------------------------------------
data(wihs, package = "randomForestSRC")
print(rfsrc.anonymous(Surv(time, status) ~ ., wihs, ntree = 100))

## ------------------------------------------------------------
## unsupervised forests
## ------------------------------------------------------------
print(rfsrc.anonymous(data = iris))

## ------------------------------------------------------------
## multivariate regression
## ------------------------------------------------------------
print(rfsrc.anonymous(Multivar(mpg, cyl) ~., data = mtcars))

## ------------------------------------------------------------
## prediction on test data with missing values using pbc data
## cases 1 to 312 have no missing values
## cases 313 to 418 having missing values
## ------------------------------------------------------------
data(pbc, package = "randomForestSRC")
pbc.obj <- rfsrc.anonymous(Surv(days, status) ~ ., pbc)
print(pbc.obj)

## mean value imputation
print(predict(pbc.obj, pbc[-(1:312),], na.action = "na.impute"))

## random imputation
print(predict(pbc.obj, pbc[-(1:312),], na.action = "na.random"))

## ------------------------------------------------------------
## train/test setting but tricky because factor labels differ over
## training and test data
## ------------------------------------------------------------

# first we convert all x-variables to factors
data(veteran, package = "randomForestSRC")
veteran.factor <- data.frame(lapply(veteran, factor))
veteran.factor$time <- veteran$time
veteran.factor$status <- veteran$status

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

# grow the forest on the training data and predict on the test data
v.grow <- rfsrc.anonymous(Surv(time, status) ~ ., veteran.factor[train, ]) 
v.pred <- predict(v.grow, veteran.factor[-train, ])
print(v.grow)
print(v.pred)

Fast Random Forests

Description

Fast approximate random forests using subsampling with forest options set to encourage computational speed. Applies to all families.

Usage

rfsrc.fast(formula, data,
  ntree = 500,
  nsplit = 10,
  bootstrap = "by.root",
  sampsize = function(x){min(x * .632, max(150, x ^ (3/4)))},
  samptype = "swor",
  samp = NULL,
  ntime = 50,
  forest = FALSE,
  save.memory = TRUE,
  ...)

Arguments

Details

Calls rfsrc by choosing options (like subsampling) to encourage computational speeds. This will provide a good approximation but will not be as good as default settings of rfsrc.

Value

An object of class (rfsrc, grow).

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

See Also

rfsrc

Examples

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

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

## do quick and *dirty* imputation
housing <- impute(SalePrice ~ ., housing,
         ntree = 50, nimpute = 1, splitrule = "random")

## grow a fast forest
o1 <- rfsrc.fast(SalePrice ~ ., housing)
o2 <- rfsrc.fast(SalePrice ~ ., housing, nodesize = 1)
print(o1)
print(o2)

## grow a fast bivariate forest
o3 <- rfsrc.fast(cbind(SalePrice,Overall.Qual) ~ ., housing)
print(o3)

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

data(wine, package = "randomForestSRC")
wine$quality <- factor(wine$quality)
o <- rfsrc.fast(quality ~ ., wine)
print(o)

## ------------------------------------------------------------
## grow fast random survival forests without C-calculation
## use brier score to assess model performance
## compare pure random splitting to logrank splitting
## ------------------------------------------------------------

data(peakVO2, package = "randomForestSRC")
f <- as.formula(Surv(ttodead, died)~.)
o1 <- rfsrc.fast(f, peakVO2, perf.type = "none")
o2 <- rfsrc.fast(f, peakVO2, perf.type = "none", splitrule = "random")
bs1 <- get.brier.survival(o1, cens.model = "km")
bs2 <- get.brier.survival(o2, cens.model = "km")
plot(bs2$brier.score, type = "s", col = 2)
lines(bs1$brier.score, type = "s", col = 4)
legend("bottomright", legend = c("random", "logrank"), fill = c(2,4))

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

data(wihs, package = "randomForestSRC")
o <- rfsrc.fast(Surv(time, status) ~ ., wihs)
print(o)

## ------------------------------------------------------------
## class imbalanced data using gmean performance
## ------------------------------------------------------------

data(breast, package = "randomForestSRC")
breast <- na.omit(breast)
f <- as.formula(status ~ .)
o <- rfsrc.fast(f, breast, perf.type = "gmean")
print(o)

## ------------------------------------------------------------
## class imbalanced data using random forests quantile-classifer (RFQ)
## fast=TRUE => rfsrc.fast
## see imbalanced function for further details
## ------------------------------------------------------------

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

Show the NEWS file

Description

Show the NEWS file of the randomForestSRC package.

Usage

rfsrc.news(...)

Arguments

Value

None.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

sidClustering using SID (Staggered Interaction Data) for Unsupervised Clustering

Description

Clustering of unsupervised data using SID (Mantero and Ishwaran, 2021). Also implements the artificial two-class approach of Breiman (2003).

Usage

## S3 method for class 'rfsrc'
sidClustering(data,
  method = "sid",
  k = NULL,
  reduce = TRUE,
  ntree = 500,
  ntree.reduce = function(p, vtry){100 * p / vtry},
  fast = FALSE,
  x.no.sid = NULL,
  use.sid.for.x = TRUE,
  x.only = NULL, y.only = NULL,
  dist.sharpen = TRUE, ...)

sid.perf.metric(truth, cluster, mode = c("entropy", "gini"))

Arguments

Details

Given an unsupervised dataset, random forests is used to compute a distance matrix measuring dissimilarity between all pairs of observations. By default, hierarchical clustering is applied to this distance matrix, although users may apply any other clustering algorithm. See the examples below for alternative workflows.

The default method, method = "sid", implements SID clustering (sidClustering). The algorithm begins by enhancing the original feature space using Staggered Interaction Data (SID). This transformation creates:

• SID main features: shifted and staggered versions of the original features that are made strictly positive and mutually non-overlapping in range;

• SID interaction features: pairwise multiplicative interactions formed between all SID main features.

A multivariate random forest is trained to predict SID main features using the SID interaction features as predictors. The rationale is that if a feature is informative for distinguishing clusters, it will exhibit systematic variation across the data space. Because each interaction feature is uniquely defined by the features it is formed from, node splits on interaction terms are able to capture and separate such variation, thus effectively identifying the clusters. See Mantero and Ishwaran (2021) for further details.

Since SID includes all pairwise interactions, the dimensionality of the feature space grows quadratically with the number of original variables (or worse when factor variables are present). As such, the reduction step using holdout variable importance (VIMP) is strongly recommended (enabled by default). This step can be disabled using reduce = FALSE, but only when the original feature space is of manageable size.

A second approach, proposed by Breiman (2003) and refined by Shi and Horvath (2006), transforms the unsupervised task into a two-class supervised classification problem. The first class consists of the original data, while the second class is generated artificially. The goal is to separate real data from synthetic data. A proximity matrix is constructed from this supervised model, and the proximity values for the original class are extracted and converted into a distance matrix (distance = 1 - proximity) for clustering.

Artificial data can be generated using two modes:

• mode 1 (default): draws random values from the empirical distribution of each feature;

• mode 2: draws uniformly between the observed minimum and maximum of each feature.

This method is invoked by setting method = "sh", "sh1", or "sh2". Mantero and Ishwaran (2021) found that while this approach works in certain settings, it can fail when clusters exist in lower-dimensional subspaces (e.g., when defined by interactions or involving both factors and continuous variables). Among the two modes, mode 1 is generally more robust.

The third method, method = "unsupv", trains a multivariate forest using the data both as predictors and as responses. The multivariate splitting rule is applied at each node. This method is fast and simple but may be less accurate compared to SID clustering.

Evaluating a clustering

sid.perf.metric measures how well a clustering separates known classes. It evaluates the mixture of true classes within each cluster: a cluster containing only one true class has zero impurity, while a cluster containing several classes has positive impurity. The true labels are supplied for evaluation after clustering. The helper can evaluate SID or any other clustering method.

Supply truth and cluster for the same observations in the same order. Cluster labels are arbitrary identifiers; they need not match the names or numerical values in truth. The number of clusters may differ from the number of true classes. For a matrix of clustering solutions, evaluate one column at a time. Use complete label vectors with unused factor levels removed.

Let K be the number of observed true classes and J the number of nonempty clusters. In the cross-tabulation table(truth, cluster), let n_{kj} count observations from true class k assigned to cluster j. Write n_j=\sum_{k=1}^{K}n_{kj} for the size of cluster j, n=\sum_{j=1}^{J}n_j for the total sample size, and p_{kj}=n_{kj}/n_j for the within-cluster class proportions.

For mode = "entropy" (the default), the impurity of cluster j is

H_j=-\sum_{k=1}^{K}p_{kj}\log_2(p_{kj}),

with zero-probability terms contributing zero. Entropy measures the remaining uncertainty about the true class within that cluster, in bits. For mode = "gini", the impurity is

G_j=1-\sum_{k=1}^{K}p_{kj}^{2}.

This is the probability that two independent draws from the cluster's empirical class distribution have different labels.

The returned result is the cluster-size-weighted average,

H=\sum_{j=1}^{J}\frac{n_j}{n}H_j \quad\mbox{or}\quad G=\sum_{j=1}^{J}\frac{n_j}{n}G_j.

Weighting by n_j/n makes a cluster's contribution proportional to the number of observations it represents, rather than giving a small cluster the same weight as a large one. The weighted entropy is the empirical conditional entropy of the true class given cluster membership. Smaller values indicate purer clusters in either mode.

Normalizing and interpreting the clustering score

The maximum impurity of a cluster with K possible true classes occurs when those classes are equally represented. The maxima are \log_2(K) for entropy and 1-1/K for Gini. Accordingly, normalized_measure is

H_{\mathrm{norm}}=\frac{H}{\log_2(K)} \quad\mbox{or}\quad G_{\mathrm{norm}}=\frac{G}{1-1/K}.

For K>1, the normalized score lies between zero and one. Zero means every cluster contains only one true class; one means every cluster contains the same proportion 1/K of each true class. The score is a fraction of the uniform-class maximum impurity, rather than a misclassification rate. With only one true class, the raw score is zero and the normalized score is undefined.

For example, suppose there are ten observations from two true classes, A and B. One cluster contains four A observations and is pure; a second contains three A and three B observations and is evenly mixed. Their weights are 0.4 and 0.6. The entropy impurities are zero and one, giving result = 0.6 and normalized_measure = 0.6. The Gini impurities are zero and 0.5, giving result = 0.3 and normalized_measure = 0.6 after division by 0.5.

The normalization uses the uniform-class maximum, not the impurity of the observed marginal class proportions or an expected score under random clustering. Consequently, when the true classes are imbalanced, a clustering that reproduces the overall class proportions in every cluster can provide no class separation and still have a normalized score below one.

These measures assess within-cluster purity. Dividing one true class into several pure clusters still gives zero impurity, and assigning each observation to its own cluster also gives zero. Thus a low score does not by itself select an appropriate number of clusters or require a one-to-one correspondence with the true classes. For method comparisons, use the same observations and true labels, and compare solutions with the same number of clusters when isolating differences in class separation. See Mantero and Ishwaran (2021) for theoretical background and empirical benchmarking.

Value

sidClustering returns a list with the following components:

Clustering performance

sid.perf.metric returns a list with three components:

result

The cluster-size-weighted entropy or Gini impurity.

measure

The selected mode, "entropy" or "gini".

normalized_measure

The weighted impurity divided by its maximum for equally represented true classes, as described in Details.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

Breiman, L. (2003). Manual on setting up, using and understanding random forest, V4.0. University of California Berkeley, Statistics Department, Berkeley.

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

Shi, T. and Horvath, S. (2006). Unsupervised learning with random forest predictors. Journal of Computational and Graphical Statistics, 15(1):118-138.

See Also

rfsrc, rfsrc.fast

Examples

## ------------------------------------------------------------
## mtcars example
## ------------------------------------------------------------

## default SID method 
o1 <- sidClustering(mtcars)
print(split(mtcars, o1$cl[, 10]))

## using artifical class approach
o1.sh <- sidClustering(mtcars, method = "sh")
print(split(mtcars, o1.sh$cl[, 10]))


## ------------------------------------------------------------
## glass data set
## ------------------------------------------------------------

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

  ## this is a supervised problem, so we first strip the class label
  data(Glass)
  glass <- Glass
  y <- Glass$Type
  glass$Type <- NULL

  ## default SID call 
  o2 <- sidClustering(glass, k = 6)
  print(table(y, o2$cl))
  print(sid.perf.metric(y, o2$cl))
  print(sid.perf.metric(y, o2$clustering, mode = "gini"))

  ## compare with Shi-Horvath mode 1 
  o2.sh <- sidClustering(glass, method = "sh1", k = 6)
  print(table(y, o2.sh$cl))
  print(sid.perf.metric(y, o2.sh$cl))

  ## plain-vanilla unsupervised analysis
  o2.un <- sidClustering(glass, method = "unsupv", k = 6)
  print(table(y, o2.un$cl))
  print(sid.perf.metric(y, o2.un$cl))

}

## ------------------------------------------------------------
## vowel data set
## ------------------------------------------------------------

if (library("mlbench", logical.return = TRUE) &&
    library("cluster", logical.return = TRUE)) {

  ## strip the class label
  data(Vowel)
  vowel <- Vowel
  y <- Vowel$Class
  vowel$Class <- NULL

  ## SID 
  o3 <- sidClustering(vowel, k = 11)
  print(table(y, o3$cl))
  print(sid.perf.metric(y, o3$cl))

  ## compare to Shi-Horvath which performs poorly in
  ## mixed variable settings
  o3.sh <- sidClustering(vowel, method = "sh1", k = 11)
  print(table(y, o3.sh$cl))
  print(sid.perf.metric(y, o3.sh$cl))

  ## Shi-Horvath improves with PAM clustering
  ## but still not as good as SID
  o3.sh.pam <- pam(o3.sh$dist, k = 11)$clustering
  print(table(y, o3.sh.pam))
  print(sid.perf.metric(y, o3.sh.pam))

  ## plain-vanilla unsupervised analysis
  o3.un <- sidClustering(vowel, method = "unsupv", k = 11)
  print(table(y, o3.un$cl))
  print(sid.perf.metric(y, o3.un$cl))

}

## ------------------------------------------------------------
##  two-d V-shaped cluster (y=x, y=-x) sitting in 12-dimensions 
##  illustrates superiority of SID to Breiman/Shi-Horvath
## ------------------------------------------------------------

p <- 10
m <- 250
n <- 2 * m
std <- .2

x <- runif(n, 0, 1)
noise <- matrix(runif(n * p, 0, 1), n)
y <- rep(NA, n)
y[1:m] <- x[1:m] + rnorm(m, sd = std)
y[(m+1):n] <- -x[(m+1):n] + rnorm(m, sd = std)
vclus <- data.frame(clus = c(rep(1, m), rep(2,m)), x = x, y = y, noise)

## SID
o4 <- sidClustering(vclus[, -1], k = 2)
print(table(vclus[, 1], o4$cl))
print(sid.perf.metric(vclus[, 1], o4$cl))

## Shi-Horvath
o4.sh <- sidClustering(vclus[, -1], method = "sh1", k = 2)
print(table(vclus[, 1], o4.sh$cl))
print(sid.perf.metric(vclus[, 1], o4.sh$cl))

## plain-vanilla unsupervised analysis
o4.un <- sidClustering(vclus[, -1], method = "unsupv", k = 2)
print(table(vclus[, 1], o4.un$cl))
print(sid.perf.metric(vclus[, 1], o4.un$cl))


## ------------------------------------------------------------
##  two-d V-shaped cluster using fast random forests
## ------------------------------------------------------------

o5 <- sidClustering(vclus[, -1], k = 2, fast = TRUE)
print(table(vclus[, 1], o5$cl))
print(sid.perf.metric(vclus[, 1], o5$cl))

Subsampling Inference for Variable Importance and Prediction Error

Description

Grow forests on repeated subsamples of the training data to estimate uncertainty in variable importance (VIMP) and out-of-bag prediction error. The saved subsample estimates can be extracted, printed, and plotted at different confidence levels without repeating the forest fits. An optional double-bootstrap procedure is also available.

Usage

## S3 method for class 'rfsrc'
subsample(obj,
  B = 100,
  block.size = 1,
  importance,
  subratio = NULL,
  stratify = TRUE,
  performance = FALSE,
  performance.only = FALSE,
  joint = FALSE,
  xvar.names = NULL,
  bootstrap = FALSE,
  verbose = TRUE)

extract.subsample(obj, alpha = .05, target = 0, m.target = NULL,
                  standardize = TRUE, raw = FALSE)

extract.bootsample(obj, alpha = .05, target = 0, m.target = NULL,
                   standardize = TRUE, raw = FALSE)

## S3 method for class 'subsample'
print(x, alpha = .05, standardize = TRUE, ...)

## S3 method for class 'bootsample'
print(x, alpha = .05, standardize = TRUE, ...)

Arguments

Details

What is resampled

For ordinary subsampling, subsample draws a smaller set of observations without replacement and grows a new forest on that data set. It repeats this operation B times. The original forest supplies the full-data reference VIMP and error statistics; the new forests supply their subsample counterparts.

The new grow calls reuse the stored formula, number of trees, mtry, node-size and depth settings, splitting rule, nsplit, predictor and split weights, sampling-size function, sampling type, and performance settings collected from obj. For survival families, the stored time grid and the choice of Uno weighting are also used. The procedure requires original root-level sampling, bootstrap = "by.root".

Outer subsampling and within-forest tree sampling have different roles. subratio determines the outer data set. The stored sampsize function is then applied to the size of that data set when growing its trees.

Subsample refits use the grow function's default missing-value omission. They do not repeat the original forest's imputation, avoiding the cost of additional imputation within each replicate. When stored training tables contain missing values, the effective refit sample can therefore be smaller than the drawn subsample. The settings listed above define the refit; specialized grow options outside that list are not recovered.

Choosing and reusing VIMP

Requesting VIMP in the original grow call avoids a separate calculation on the saved trees. For example, specify importance = "permute" there to obtain uncertainty for permutation VIMP. Existing VIMP takes precedence over a different importance or block.size supplied to ordinary subsample.

joint = TRUE adds a row named joint; it measures the effect of perturbing the selected predictors together. performance = TRUE adds a row named err; performance.only = TRUE retains only that error row. These rows use the same interval machinery as the individual VIMP rows. A positive lower bound for an error rate is not a variable-selection result.

Subsample size and stratification

The default fraction is constant through n=1000; it decreases in proportion to n^{-1/2} only above that threshold. With stratification, integer sample sizes are allocated separately within each stratum. Small strata can receive no observations; stratification does not guarantee that every class or event type appears in each new forest. Missing stratum labels are excluded from the allocation.

For survival and competing risks, each drawn subsample must retain at least one observed event after the usual missing-value omission. If it does not, the procedure stops with a message recommending a larger subratio. No class-presence check is added for classification. The allocation rule is unchanged; the procedure does not force one observation into every stratum or redraw a rejected subsample.

The extractor uses the nominal size m=nr, where r is subratio, for scaling. The actual integer sample count, which can be smaller after stratification or missing-data processing, is not recorded in the returned object.

Extraction, confidence levels, and standardization

For an ordinary subsampling object, extract.subsample selects one outcome and returns all three interval summaries: nonparametric subsampling, normal intervals using the subsampling standard error, and normal intervals using the jackknife standard error. For multivariate or mixed outcomes, m.target selects an outcome in the stored results. It is not a new forest-prediction request.

The extractor and print method default to alpha = .05; plot.subsample defaults to alpha = .01. Set the same alpha explicitly to compare printed and plotted intervals. Each interval has nominal confidence level 1-\alpha; the calculations are made separately for each statistic, without a multiple-testing adjustment.

Extraction, printing, and plotting use standardize = TRUE by default. For a regression outcome, both full-data and subsample statistics are divided by the variance of that outcome in the full training data. The same variance is used for every replicate, and error rows are scaled in the same way as VIMP rows. Classification and survival statistics retain their original scale. Use standardize = FALSE for unscaled values. Regression standardization is defined for a finite, positive response variance.

raw = TRUE in extract.subsample requests additional output, including the replicate estimates and interval matrices. It does not turn off standardization.

The three subsampling intervals

The same subsample fits provide two variance estimators: the subsampling estimator, which centers the subsample estimates at their own mean, and the delete-d jackknife estimator, which centers them at the full-data estimate. A subsample retaining m observations from n can also be viewed as deleting d=n-m observations. The three interval summaries are normal intervals using either variance estimator and nonparametric intervals using the centered subsampling distribution.

For one statistic, let \widehat\theta be the full-data estimate and \theta_b^* its estimate in subsample b. Write r=m/n, let \bar\theta_* be the mean of the available subsample estimates, and let s_* be their sample standard deviation.

The subsampling standard error, returned as se.Z, is

\mathrm{SE}_{\mathrm{sub}}=\sqrt{r}\,s_*.

It measures dispersion about \bar\theta_*. The jackknife standard error, returned as se.jk.Z, uses squared deviations from \widehat\theta:

\mathrm{SE}_{\mathrm{jk}}= \left\{r\,\mathrm{mean}_b (\theta_b^*-\widehat\theta)^2\right\}^{1/2}.

This full-data centering retains both the subsample dispersion and the squared difference between the subsample mean and the full-data estimate:

\mathrm{mean}_b(\theta_b^*-\widehat\theta)^2 =\mathrm{mean}_b(\theta_b^*-\bar\theta_*)^2 +(\bar\theta_*-\widehat\theta)^2.

The squared mean displacement supplies the nonnegative bias adjustment discussed in Ishwaran and Lu (2019), equations (10) and (11).

Both normal intervals use \widehat\theta\pm z_{1-\alpha/2}\mathrm{SE}. The subsampling version is returned in var.sel.Z; the jackknife version is returned in var.jk.sel.Z.

The formulas above describe the implementation's scaling: se.jk.Z uses m/n, whereas the delete-d jackknife variance in equation (11) uses m/(n-m). These scaling factors are asymptotically equivalent when m/n\to0. Also, for K available estimates, se.Z uses the sample-variance denominator K-1, while se.jk.Z averages squared deviations with denominator K. With these finite-replicate conventions, the returned jackknife standard error need not exceed the returned subsampling standard error.

For the nonparametric interval, define the roots Z_b^*=\sqrt{m}(\theta_b^*-\widehat\theta) and let q_p be their empirical quantile, calculated using the default quantile interpolation. Its endpoints are

\left[\widehat\theta-q_{1-\alpha/2}/\sqrt{n}, \widehat\theta-q_{\alpha/2}/\sqrt{n}\right].

Thus the nonparametric interval uses reversed, centered-root quantiles rather than the unscaled percentiles of the replicate VIMP values. These formulas describe the subsampling path; double-bootstrap objects use their own extractor.

Double-bootstrap output

With bootstrap = TRUE, each outer sample draws n observations with replacement. The per-tree samples are then generated from that bootstrap multiset. Under samptype = "swor", sampling is without replacement from multiset positions, so repeated original observation identities remain possible.

extract.subsample dispatches double-bootstrap objects to extract.bootsample. This returns percentile intervals and normal intervals centered at the bootstrap mean, with the ordinary bootstrap standard deviation. The nonparametric bootstrap table has no p-value column. The centered-root p-value above is specific to subsampling without replacement.

Reading extracted summaries

The default extraction returns se.Z, se.jk.Z, var.sel, var.sel.Z, and var.jk.sel.Z. The three tables contain interval endpoints and a signif column indicating that the lower endpoint exceeds zero. var.sel reports an interval median; the normal tables report the full-data estimate in their mean column.

The normal tables report one-sided upper-tail p-values for zero VIMP. Their positive-lower-bound rule therefore corresponds to a one-sided level of \alpha/2, not \alpha. The nonparametric pvalue uses the same centered roots as its interval. It is the proportion of available roots satisfying Z_b^*>\sqrt{n}\,\widehat\theta, a one-sided upper-tail comparison for zero importance. Ties are excluded by the strict inequality. No Monte Carlo adjustment is added. The interval still uses the default quantile interpolation, so empirical p-values and interpolated interval endpoints need not yield identical decisions at every finite-sample boundary.

With raw = TRUE, the extractor also returns vmp, vmpS, ci, ci.Z, ci.jk.Z, and boxplot.dta. The replicate matrix vmpS has one row per statistic and one column per replicate. The five rows of each interval matrix contain its lower endpoint, lower quartile, center, upper quartile, and upper endpoint. boxplot.dta contains separately centered and rescaled replicate values.

Printing interval summaries

print(x, alpha = .05, standardize = TRUE) displays interval matrices for every stored outcome and statistic column, without fitting more forests. Subsampling output includes the nonparametric, normal subsampling, and normal jackknife intervals. Double-bootstrap output includes its nonparametric and normal intervals. The displayed values are rounded to three decimal places.

print.subsample.rfsrc and print.bootsample.rfsrc are equivalent to print.subsample and print.bootsample, respectively, with the same arguments. Use the extraction helpers to select one outcome or to obtain numerical results at full precision.

Value

subsample returns a list containing the full-data reference and the replicate estimates. Interval summaries are obtained by extraction, printing, or plotting.

The returned class includes subsample for subsampling or bootsample for double bootstrapping, together with the retained forest classes.

Extraction helpers

Both helpers return a list for one outcome and statistic column. Summary tables have one row per predictor or requested joint/error statistic. They reuse the saved replicates without fitting forests.

extract.subsample

For ordinary subsampling, returns se.Z, se.jk.Z, var.sel, var.sel.Z, and var.jk.sel.Z. With raw = TRUE, also returns vmp, vmpS, ci, ci.Z, ci.jk.Z, and boxplot.dta, as described in Details. For a double-bootstrap object, returns the extract.bootsample result below.

extract.bootsample

Returns se, the bootstrap standard errors; var.sel, the percentile-interval table; and var.sel.Z, the normal-interval table. Both tables contain lower, mean, upper, and signif; only var.sel.Z includes a pvalue column. With raw = TRUE, also returns the bootstrap means vmp, replicate matrix vmpS, and five-row interval matrices ci and ci.Z. Both interval matrices use the bootstrap mean in their middle row.

The vmpS matrix has statistics in rows and replicates in columns for either procedure. All extracted estimates use the requested standardize setting.

The print methods are used for their display side effect. Use extract.subsample or extract.bootsample for the numerical interval summaries.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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.

Geyer, C.J. (2013). 5601 Notes: The Subsampling Bootstrap. July 11, 2013. Section 5, equation (12).

Politis, D.N. and Romano, J.P. (1994). Large sample confidence regions based on subsamples under minimal assumptions. The Annals of Statistics, 22(4):2031-2050.

Shao, J. and Wu, C.J. (1989). A general theory for jackknife variance estimation. The Annals of Statistics, 17(3):1176-1197.

See Also

plot.subsample.rfsrc, rfsrc, vimp.rfsrc, holdout.vimp.rfsrc

Examples

## Small settings keep these examples manageable. Increase B for more
## stable interval endpoints, especially at high confidence levels.
set.seed(19)
dta <- na.omit(airquality)

## ------------------------------------------------------------
## Regression: grow once, then reuse the saved replicate estimates
## ------------------------------------------------------------
reg.o <- rfsrc(Ozone ~ ., data = dta, ntree = 100,
               importance = "permute", block.size = 1)
reg.smp <- subsample(reg.o, B = 25, verbose = FALSE)

ans <- extract.subsample(reg.smp, alpha = .05)
print(ans$var.sel.Z)
print(ans$var.jk.sel.Z)
print(reg.smp, alpha = .05)
plot.subsample(reg.smp, alpha = .05)
plot.subsample(reg.smp, alpha = .05, normal = FALSE)

## Detailed output and unstandardized VIMP are separate requests.
raw <- extract.subsample(reg.smp, alpha = .05,
                         standardize = FALSE, raw = TRUE)
print(raw$ci.Z)
print(dim(raw$vmpS))

## ------------------------------------------------------------
## Joint VIMP and prediction-error uncertainty
## ------------------------------------------------------------
with.error <- subsample(reg.o, B = 25, performance = TRUE,
                         joint = TRUE, xvar.names = c("Day", "Month"),
                         verbose = FALSE)
ans <- extract.subsample(with.error, alpha = .05)
print(ans$var.sel.Z[c("joint", "err"), , drop = FALSE])

error.only <- subsample(reg.o, B = 25, performance.only = TRUE,
                         verbose = FALSE)
print(extract.subsample(error.only, alpha = .05)$var.sel.Z)

## ------------------------------------------------------------
## Multivariate regression: select an outcome after subsampling
## ------------------------------------------------------------
mv <- rfsrc(cbind(Ozone, Temp) ~ ., data = dta, ntree = 100,
             importance = "permute", block.size = 1)
mv.smp <- subsample(mv, B = 25, verbose = FALSE)
print(extract.subsample(mv.smp, m.target = "Temp", alpha = .05)$var.sel.Z)
plot.subsample(mv.smp, m.target = "Ozone", alpha = .05)
plot.subsample(mv.smp, m.target = "Temp", alpha = .05,
                jknife = TRUE)

## ------------------------------------------------------------
## Classification with a different performance measure
## ------------------------------------------------------------
if (requireNamespace("caret", quietly = TRUE)) {
  cls.data <- caret::twoClassSim(500, linearVars = 5, noiseVars = 5)
  cls <- rfsrc(Class ~ ., data = cls.data, ntree = 100,
                importance = "permute", block.size = 1,
                perf.type = "brier")
  cls.smp <- subsample(cls, B = 25, verbose = FALSE)
  print(extract.subsample(cls.smp, alpha = .05)$var.sel.Z)
  plot.subsample(cls.smp, alpha = .05, cex.axis = .7)
}

## ------------------------------------------------------------
## Survival
## ------------------------------------------------------------
data(pbc, package = "randomForestSRC")
pbc.complete <- na.omit(pbc)
srv <- rfsrc(Surv(days, status) ~ ., data = pbc.complete,
               ntree = 100, importance = "permute", block.size = 1)
srv.smp <- subsample(srv, B = 25, verbose = FALSE)
plot.subsample(srv.smp, alpha = .05)

## ------------------------------------------------------------
## Optional double bootstrap (separate resampling procedure)
## ------------------------------------------------------------
reg.boot <- subsample(reg.o, B = 25, bootstrap = TRUE, verbose = FALSE)
print(reg.boot, alpha = .05)
print(extract.bootsample(reg.boot, alpha = .05)$var.sel.Z)

Tune Random Forest for optimal mtry and nodesize

Description

Finds the optimal mtry and nodesize for a random forest using out-of-bag (OOB) error. Two search strategies are supported: a grid-based search and a golden-section search with noise control. Works for all response families supported by rfsrc.fast.

Usage

## S3 method for class 'rfsrc'
tune(formula, data,
  mtry.start = ncol(data) / 2,
  nodesize.try = c(1:9, seq(10, 100, by = 5)), ntree.try = 100,
  sampsize = function(x) { min(x * .632, max(150, x^(3/4))) },
  nsplit = 1, step.factor = 1.25, improve = 1e-3, strikeout = 3, max.iter = 25,
  method = c("grid", "golden"),
  final.window = 5, reps.initial = 2, reps.final = 3,
  trace = FALSE, do.best = TRUE, seed = NULL, ...)

## S3 method for class 'rfsrc'
tune.nodesize(formula, data,
  nodesize.try = c(1:9, seq(10, 150, by = 5)), ntree.try = 100,
  sampsize = function(x) { min(x * .632, max(150, x^(4/5))) },
  nsplit = 1, method = c("grid", "golden"),
  final.window = 5, reps.initial = 2, reps.final = 3, max.iter = 50,
  trace = TRUE, seed = NULL, ...)

Arguments

Details

Error estimate. If 2 * ssize < n, a disjoint holdout of size ssize is used for evaluation; otherwise OOB error is used.

Subsample used during tuning. Both functions derive a single integer ssize from sampsize and pass it to rfsrc.fast for all tuning fits. This improves stability and comparability across candidates. When do.best = TRUE in tune, the final forest is fit with the user-supplied sampsize exactly as provided.

Grid search. tune performs a step-out search over mtry for each nodesize in nodesize.try, using step.factor, improve, strikeout, and max.iter. tune.nodesize evaluates the supplied nodesize.try grid directly.

Golden search. Uses a guarded golden-section line search with noise control. For each one-dimensional search (over nodesize or mtry), the routine probes a small left-anchor grid 1:9, iterates golden shrinkage until the bracket width is at most final.window, then runs a short local sweep with reps.final replicates. In tune the searches over nodesize and mtry alternate in a simple coordinate loop, with improve and strikeout as stopping controls.

Value

For tune:

• results: matrix with columns nodesize, mtry, err.

• optimal: named numeric vector c(nodesize = ..., mtry = ...).

• rf: fitted forest at the optimum if do.best = TRUE.

For tune.nodesize:

• nsize.opt: optimal nodesize.

• err: data frame with columns nodesize and err.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

See Also

rfsrc.fast

Examples

## ------------------------------------------------------------
## White wine classification example
## ------------------------------------------------------------
data(wine, package = "randomForestSRC")
wine$quality <- factor(wine$quality)

## Fixed seed makes tuning reproducible
set.seed(1)

## Full tuner over nodesize and mtry (grid)
o1 <- tune(quality ~ ., wine, sampsize = 100, method = "grid")
print(o1$optimal)

## Golden search alternative
o2 <- tune(quality ~ ., wine, sampsize = 100, method = "golden",
           reps.initial = 2, reps.final = 3, seed = 1)
print(o2$optimal)

## visualize the nodesize/mtry surface
if (library("interp", logical.return = TRUE)) {

  plot.tune <- function(o, linear = TRUE) {
    x <- o$results[, 1]
    y <- o$results[, 2]
    z <- o$results[, 3]
    so <- interp(x = x, y = y, z = z, linear = linear)
    idx <- which.min(z)
    x0 <- x[idx]; y0 <- y[idx]
    filled.contour(x = so$x, y = so$y, z = so$z,
                   xlim = range(so$x, finite = TRUE) + c(-2, 2),
                   ylim = range(so$y, finite = TRUE) + c(-2, 2),
                   color.palette = colorRampPalette(c("yellow", "red")),
                   xlab = "nodesize", ylab = "mtry",
                   main = "error rate for nodesize and mtry",
                   key.title = title(main = "OOB error", cex.main = 1),
                   plot.axes = {
                     axis(1); axis(2)
                     points(x0, y0, pch = "x", cex = 1, font = 2)
                     points(x, y, pch = 16, cex = .25)
                   })
  }

  plot.tune(o1)
  plot.tune(o2)
}

## ------------------------------------------------------------
## nodesize only: grid vs golden
## ------------------------------------------------------------
o3 <- tune.nodesize(quality ~ ., wine, sampsize = 100, method = "grid",
                    trace = TRUE, seed = 1)
o4 <- tune.nodesize(quality ~ ., wine, sampsize = 100, method = "golden",
                    reps.initial = 2, reps.final = 3, trace = TRUE, seed = 1)
plot(o3$err, type = "s", xlab = "nodesize", ylab = "error")

## ------------------------------------------------------------
## Tuning for class imbalance (rfq with geometric mean performance)
## ------------------------------------------------------------
data(breast, package = "randomForestSRC")
breast <- na.omit(breast)
o5 <- tune(status ~ ., data = breast, rfq = TRUE, perf.type = "gmean",
           method = "golden", seed = 1)
print(o5$optimal)

## ------------------------------------------------------------
## Competing risks example (nodesize only)
## ------------------------------------------------------------
data(wihs, package = "randomForestSRC")
plot(tune.nodesize(Surv(time, status) ~ ., wihs, trace = TRUE)$err, type = "s")

van de Vijver Microarray Breast Cancer

Description

Gene expression profiling for predicting clinical outcome of breast cancer (van't Veer et al., 2002). Microarray breast cancer data set of 4707 expression values on 78 patients with survival information.

References

van't Veer L.J. et al. (2002). Gene expression profiling predicts clinical outcome of breast cancer. Nature, 12, 530–536.

Examples

data(vdv, package = "randomForestSRC")

Veteran's Administration Lung Cancer Trial

Description

Randomized trial of two treatment regimens for lung cancer. This is a standard survival analysis data set.

Source

Kalbfleisch and Prentice, The Statistical Analysis of Failure Time Data.

References

Kalbfleisch J. and Prentice R, (1980) The Statistical Analysis of Failure Time Data. New York: Wiley.

Examples

data(veteran, package = "randomForestSRC")

VIMP for Single or Grouped Variables

Description

Calculate variable importance (VIMP) for a single variable or group of variables for training or test data.

Usage

## S3 method for class 'rfsrc'
vimp(object, xvar.names, 
  importance = c("anti", "permute", "random"), block.size = 10,
  joint = FALSE, seed = NULL, do.trace = FALSE, ...)

Arguments

Details

Using a previously trained forest, this function calculates variable importance (VIMP) for the specified variables in xvar.names. By default, VIMP is computed using the original training data, but the user may supply a new test set via the newdata argument. See rfsrc for further details on how VIMP is computed.

If joint = TRUE, joint VIMP is returned. This is defined as the importance of a group of variables when the entire group is perturbed simultaneously.

Setting csv = TRUE returns case-specific VIMP, which provides VIMP estimates at the individual observation level. This applies to all families except survival. See examples below.

Value

An object of class (rfsrc, predict) containing importance values.

Author(s)

Hemant Ishwaran and Udaya B. Kogalur

References

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

See Also

holdout.vimp.rfsrc, rfsrc

Examples

## ------------------------------------------------------------
## classification example
## showcase different vimp
## ------------------------------------------------------------

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

## anti vimp (default)
print(vimp(iris.obj)$importance)

## anti vimp using brier prediction error
print(vimp(iris.obj, perf.type = "brier")$importance)

## permutation vimp
print(vimp(iris.obj, importance = "permute")$importance)

## random daughter vimp
print(vimp(iris.obj, importance = "random")$importance)

## joint anti vimp 
print(vimp(iris.obj, joint = TRUE)$importance)

## paired anti vimp
print(vimp(iris.obj, c("Petal.Length", "Petal.Width"), joint = TRUE)$importance)
print(vimp(iris.obj, c("Sepal.Length", "Petal.Width"), joint = TRUE)$importance)

## ------------------------------------------------------------
## survival example
## anti versus permute VIMP with different block sizes
## ------------------------------------------------------------

data(pbc, package = "randomForestSRC")
pbc.obj <- rfsrc(Surv(days, status) ~ ., pbc)

print(vimp(pbc.obj)$importance)
print(vimp(pbc.obj, block.size=1)$importance)
print(vimp(pbc.obj, importance="permute")$importance)
print(vimp(pbc.obj, importance="permute", block.size=1)$importance)

## ------------------------------------------------------------
## imbalanced classification example
## see the imbalanced function for more details
## ------------------------------------------------------------

data(breast, package = "randomForestSRC")
breast <- na.omit(breast)
f <- as.formula(status ~ .)
o <- rfsrc(f, breast, ntree = 2000)

## permutation vimp
print(100 * vimp(o, importance = "permute")$importance)

## anti vimp using gmean performance
print(100 * vimp(o, perf.type = "gmean")$importance[, 1])

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

airq.obj <- rfsrc(Ozone ~ ., airquality)
print(vimp(airq.obj))

## ------------------------------------------------------------
## regression example where vimp is calculated on test data
## ------------------------------------------------------------

set.seed(100080)
train <- sample(1:nrow(airquality), size = 80)
airq.obj <- rfsrc(Ozone~., airquality[train, ])

## training data vimp
print(airq.obj$importance)
print(vimp(airq.obj)$importance)

## test data vimp
print(vimp(airq.obj, newdata = airquality[-train, ])$importance)

## ------------------------------------------------------------
## case-specific vimp
## returns VIMP for each case
## ------------------------------------------------------------

o <- rfsrc(mpg~., mtcars)
v <- vimp(o, csv = TRUE)
csvimp <- get.mv.csvimp(v, standardize=TRUE)
print(csvimp)

## ------------------------------------------------------------
## case-specific joint vimp
## returns joint VIMP for each case
## ------------------------------------------------------------

o <- rfsrc(mpg~., mtcars)
v <- vimp(o, joint = TRUE, csv = TRUE)
csvimp <- get.mv.csvimp(v, standardize=TRUE)
print(csvimp)

## ------------------------------------------------------------
## case-specific joint vimp for multivariate regression
## returns joint VIMP for each case, for each outcome
## ------------------------------------------------------------

o <- rfsrc(Multivar(mpg, cyl) ~., data = mtcars)
v <- vimp(o, joint = TRUE, csv = TRUE)
csvimp <- get.mv.csvimp(v, standardize=TRUE)
print(csvimp)

Women's Interagency HIV Study (WIHS)

Description

Competing risk data set involving AIDS in women.

Format

A data frame containing:

Source

Study included 1164 women enrolled in WIHS, who were alive, infected with HIV, and free of clinical AIDS on December, 1995, when the first protease inhibitor (saquinavir mesylate) was approved by the Federal Drug Administration. Women were followed until the first of the following occurred: treatment initiation, AIDS diagnosis, death, or administrative censoring (September, 2006). Variables included history of injection drug use at WIHS enrollment, whether an individual was African American, age, and CD4 nadir prior to baseline.

References

Bacon M.C, von Wyl V., Alden C., et al. (2005). The Women's Interagency HIV Study: an observational cohort brings clinical sciences to the bench, Clin Diagn Lab Immunol, 12(9):1013-1019.

Examples

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

White Wine Quality Data

Description

The inputs include objective tests (e.g. PH values) and the output is based on sensory data (median of at least 3 evaluations made by wine experts) of white wine. Each expert graded the wine quality between 0 (very bad) and 10 (very excellent).

References

Cortez, P., Cerdeira, A., Almeida, F., Matos T. and Reis, J. (2009). Modeling wine preferences by data mining from physicochemical properties. In Decision Support Systems, Elsevier, 47(4):547-553.

Examples

## load wine and convert to a multiclass problem
data(wine, package = "randomForestSRC")
wine$quality <- factor(wine$quality)