subsample.rfsrc.RdEstimate standard errors and confidence intervals for variable importance (VIMP) and out-of-bag prediction error by fitting forests to repeated subsamples of the training data. Extract, print, or plot the results at different confidence levels from the same subsample estimates. An optional double-bootstrap procedure is also available.
# 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, ...)For subsample, a forest grow object of class
(rfsrc, grow) with saved forest information and training data.
Grow the forest with forest = TRUE and
bootstrap = "by.root". Anonymous forests and nonuniform
case weights are not supported. For extract.subsample, an
object returned by subsample, including double-bootstrap
output. For extract.bootsample, an object returned by
subsample(..., bootstrap = TRUE).
Number of subsamples, or double-bootstrap replicates when
bootstrap = TRUE. Each replicate requires a new forest fit.
Supply a positive integer; at least two replicates are needed to
calculate a sample standard deviation.
Number of trees per block for calculating VIMP
when it must first be computed for obj. If obj already
contains VIMP, its block size is used instead.
VIMP perturbation method, such as "anti",
"permute", or "random". For subsampling, existing VIMP
in obj determines the method. If VIMP is absent, it is
computed using this argument, or importance = TRUE when
omitted. Ignored when performance.only = TRUE.
Fraction of training observations to draw for each
subsample. Supply a finite number greater than zero and at most one;
for inference, choose a fraction below one that gives enough
observations to fit each forest. The default is
\(\exp(-1)\) for \(n\leq1000\) and
\(\exp(-1)\sqrt{1000/n}\) for larger \(n\).
This controls the data set for each forest, not its per-tree
sampsize. Ignored when bootstrap = TRUE.
Sample separately within response classes for
univariate classification, or within event-status groups, including
censoring, for survival and competing risks. Regression,
multivariate, and mixed-outcome forests use unstratified subsampling.
Ignored when bootstrap = TRUE.
Include uncertainty estimates for OOB prediction
error as well as VIMP. The error measure is determined by the
original forest's perf.type.
Compute uncertainty for OOB prediction error
only. Sets performance = TRUE and turns off individual and
joint VIMP.
Include joint VIMP for the predictor group specified by
xvar.names, in addition to individual VIMP. Joint VIMP uses
the default anti perturbation of vimp, independently of the
method used for individual VIMP.
Predictors to perturb together for joint VIMP.
The default NULL uses all predictors. Individual VIMP is
still computed for all predictors. To select results for display,
use xvar.names in plot.subsample.
Use the double-bootstrap procedure instead of subsampling without replacement. See Details.
Print progress information during resampling.
Significance level for extraction or printing, strictly
between zero and one. Intervals have nominal confidence level
\(1-\alpha\); the default .05 gives 95 percent intervals.
Statistic to extract within the selected response.
Supply a nonnegative integer: 0 for regression or
right-censored survival; 0 for overall classification or
1, 2, and so on for its classes; and 0,
1, and so on for competing-risk events. Classes and events
follow the order in the forest. For competing risks, the first event
uses target = 0 here but target = 1 in
plot.subsample.
Name of one response to extract from a multivariate or mixed-outcome forest. The default selects the first regression response with results, or the first classification response with results if none is available for regression. Not needed for univariate forests.
For regression, divide estimates by the response
variance in the full training data. This applies to VIMP, prediction
error, and their uncertainty summaries. Other families are unchanged.
Set FALSE for the original scale.
Include replicate estimates and interval matrices in the
extraction result. This controls the amount of output, not its
scale; use standardize = FALSE for unscaled values.
An object returned by subsample. Use
print.subsample for subsampling and print.bootsample
for double-bootstrap output.
Additional arguments to the print methods; currently unused.
Use subsample() to add standard errors and confidence intervals
to a forest's VIMP estimates. The original forest provides the
full-data estimate for each predictor. Forests grown on repeated
subsamples show how these estimates vary when the training
observations change. The same approach can estimate uncertainty in
OOB prediction error.
The analysis has two stages: subsample() fits the replicate
forests, then extract.subsample(), print(), or
plot.subsample() uses their saved estimates to construct
intervals. Once resampling is complete, you can change the confidence
level, interval method, or response without fitting more forests.
By default, the analysis computes uncertainty for individual VIMP.
Specify the desired VIMP method in the original grow call, for
example importance = "permute". When obj already
contains VIMP, subsampling uses that method and block size, even
if different values are supplied to subsample(). Otherwise,
VIMP is first calculated on the saved trees using the requested
importance and block.size.
Two additional statistics can be included:
joint = TRUE adds a row named
joint, measuring the effect of perturbing the predictors
in xvar.names together. Individual VIMP is also included.
performance = TRUE adds a row
named err for OOB prediction error. Use
performance.only = TRUE to compute only this row.
Each replicate draws training observations without replacement and
grows a new forest on that subset. B controls the number of
replicates; subratio controls their nominal sample size.
Increase B for more stable interval endpoints, especially at
high confidence levels.
Subsampling the data and sampling observations for individual trees
are separate operations. subratio determines the data set
available to each new forest. The original forest's sampsize
function is then applied to that data set's size to determine its
per-tree sample size. The forest settings used for these fits are
listed in the Note.
With stratify = TRUE, univariate classification samples
within classes, and survival and competing risks sample within
event-status groups. Integer allocation can leave small groups
unrepresented; missing group labels are excluded. Other families
sample without stratification.
For survival and competing risks, a subsample must contain at least
one observed event after missing-value omission. Otherwise the
procedure stops and asks for a larger subratio; it does not
draw a replacement subsample.
After smp <- subsample(obj),
extract.subsample(smp, alpha = .05) returns all three
subsampling interval summaries for one response: nonparametric,
normal subsampling, and normal jackknife. For multivariate or mixed
outcomes, use m.target to choose the response. For
classification or competing risks, target selects its
overall, class-specific, or event-specific statistic as described
under Arguments.
print(smp, alpha = .05) displays intervals for every response
and statistic column, rounded to three decimal places.
plot.subsample(smp, alpha = .05) draws the normal subsampling
intervals; its normal and jknife arguments select the
other methods. Use extraction for numerical results at full
precision, and raw = TRUE for the replicate estimates and
interval matrices listed under Value.
Extraction and printing default to 95 percent intervals
(alpha = .05); plotting defaults to 99 percent
(alpha = .01). Supply the same alpha when comparing
them. Intervals are calculated separately for each statistic,
without a multiple-testing adjustment.
Extraction, printing, and plotting use standardize = TRUE
by default. For regression, the full-data and replicate estimates
are divided by the response variance in the full training data.
This same divisor is used for every replicate and for both VIMP
and prediction error. Confidence limits and standard errors are
then on the corresponding standardized scale.
Use standardize = FALSE for the original scale. Classification
and survival statistics are unchanged by this option. Regression
standardization requires a finite, positive response variance.
For one statistic, let \(\widehat\theta\) be its full-data estimate
and \(\theta_b^*\) its estimate in subsample \(b\). Let \(n\)
be the training sample size, \(r\) the value of subratio,
and \(m=nr\) the nominal subsample size. Write \(\bar\theta_*\)
and \(s_*\) for the mean and sample standard deviation of the
available subsample estimates.
The subsampling standard
error measures dispersion about the subsample mean:
$$\mathrm{SE}_{\mathrm{sub}}=\sqrt{r}\,s_*.$$
It is returned as se.Z. The interval is
\(\widehat\theta\pm z_{1-\alpha/2}\mathrm{SE}_{\mathrm{sub}}\)
and is reported in var.sel.Z.
A subsample of size \(m\)
also deletes \(d=n-m\) observations. The delete-\(d\)
jackknife calculation measures dispersion about the full-data
estimate:
$$\mathrm{SE}_{\mathrm{jk}}=
\left\{r\,\mathrm{mean}_b
(\theta_b^*-\widehat\theta)^2\right\}^{1/2}.$$
It is returned as se.jk.Z. The interval is
\(\widehat\theta\pm z_{1-\alpha/2}\mathrm{SE}_{\mathrm{jk}}\)
and is reported in var.jk.sel.Z.
These use the
empirical distribution of the centered roots
\(Z_b^*=\sqrt{m}(\theta_b^*-\widehat\theta)\). If \(q_p\)
is their empirical quantile, computed with the default
quantile interpolation, the endpoints are
$$\left[\widehat\theta-q_{1-\alpha/2}/\sqrt{n},
\widehat\theta-q_{\alpha/2}/\sqrt{n}\right].$$
The interval is reported in var.sel. Centering, scaling,
and reversing the root quantiles converts subsample variation
to an interval for the full-data statistic; this is not a
percentile interval of the unscaled VIMP replicates.
The jackknife calculation includes both subsample variation and the 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 final term is the nonnegative bias adjustment discussed by Ishwaran and Lu (2019), equations (10) and (11). The Note explains the scaling and finite-replicate denominators used here.
Each table gives interval endpoints and a signif column
indicating whether the lower endpoint exceeds zero. For VIMP,
this identifies positive importance under the selected interval
method. An err row instead describes uncertainty in prediction
error; its positive lower endpoint is not a variable-selection
result. The normal tables report the full-data estimate in
mean; the nonparametric table reports an interval median.
Normal-table p-values are one-sided upper-tail tests of zero VIMP. A positive lower endpoint of a two-sided \(1-\alpha\) normal interval corresponds to a one-sided level of \(\alpha/2\). The nonparametric p-value is the proportion of available roots satisfying \(Z_b^*>\sqrt{n}\,\widehat\theta\). Ties are excluded, and no Monte Carlo adjustment is added. Because interval endpoints use interpolated quantiles, the empirical p-value and interval can give different decisions at finite-sample boundaries.
Set bootstrap = TRUE to draw \(n\) observations with
replacement for each outer replicate, then generate the per-tree
samples from that bootstrap data set. With samptype = "swor",
trees sample without replacement from its rows; an original
observation can still appear more than once because the bootstrap
data set can contain duplicates.
For these objects, extract.subsample() calls
extract.bootsample(). It returns percentile intervals and
normal intervals centered at the bootstrap mean, using the
bootstrap standard deviation. Printing displays both methods.
The percentile table has no p-value column; the centered-root
p-value above applies only to subsampling without replacement.
Subsample forests use the original 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. Survival forests also use
the saved time grid and choice of Uno weighting. Specialized grow
options outside these settings are not reproduced.
Each new fit omits rows with missing values rather than repeating
the original imputation. Missing-value omission and integer
allocation within strata can reduce the actual subsample size.
Interval scaling uses the nominal size \(m=nr\), where \(r\)
is subratio; actual sample counts are not saved in the result.
se.jk.Z uses the factor \(m/n\); equation (11) of Ishwaran
and Lu (2019) uses \(m/(n-m)\). These factors are asymptotically
equivalent when \(m/n\to0\). With \(K\) available estimates,
se.Z uses the sample-variance denominator \(K-1\), whereas
se.jk.Z averages squared deviations with denominator \(K\).
Consequently, the returned jackknife standard error can be smaller
than the subsampling standard error for a finite set of replicates.
subsample() returns the reference forest and replicate
estimates. Use extraction, printing, or plotting to obtain intervals.
Its class includes subsample for subsampling or
bootsample for double bootstrapping, together with the
original forest classes.
Reference forest with full-data statistics. If VIMP was initially absent, it has been calculated on the saved trees.
List of full-data statistic matrices, one per outcome.
Rows correspond to predictors and any requested joint or
err statistic. Classification and competing risks can have
several columns. Unavailable appended statistics are NA
where their outcome and column can be matched. These values are
not standardized.
For subsampling, a length-B list of replicate
results, each organized by outcome. The tables use the outcomes,
row names, and column names in vmp; unavailable entries are
NA.
Double-bootstrap replicate estimates, in place of
vmpS.
Nominal subsample fraction, or NULL for
double-bootstrap objects.
Whether only prediction error was requested.
Each helper returns a list for one response and statistic column.
Summary tables have one row per predictor or requested joint/error
statistic. All extracted values use the requested
standardize setting.
extract.subsample()For subsampling, returns
se.Z and se.jk.Z, the subsampling and jackknife
standard errors; var.sel, the nonparametric table;
var.sel.Z, the normal subsampling table; and
var.jk.sel.Z, the normal jackknife table. For a
double-bootstrap object, it returns the
extract.bootsample() result.
extract.bootsample()Returns se, the
bootstrap standard errors; var.sel, the percentile
table; and var.sel.Z, the normal table. Both tables
contain lower, mean, upper, and
signif; only var.sel.Z includes pvalue.
With raw = TRUE, the helpers also return:
vmpFull-data estimates for subsampling, or bootstrap means for double bootstrapping.
vmpSMatrix with statistics in rows and replicates in columns, for either procedure.
ci, ci.Z, ci.jk.ZFive-row interval
matrices for the nonparametric, normal subsampling, and normal
jackknife methods. Rows give the lower endpoint, lower
quartile, center, upper quartile, and upper endpoint.
Double-bootstrap output includes only ci and
ci.Z, both with the bootstrap mean in the middle row.
boxplot.dtaCentered and rescaled replicate values used for subsampling plots. Not returned for double bootstrapping.
The print methods display interval summaries. Use the extraction
helpers for numerical results. print.subsample.rfsrc and
print.bootsample.rfsrc have the same arguments as
print.subsample and print.bootsample, respectively.
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.
# \donttest{
## 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)
# }