varpro: Model-Independent Variable Selection via the Rule-Based...

In varPro: Model-Independent Variable Selection via the Rule-Based Variable Priority

Model-Independent Variable Selection via the Rule-Based Variable Priority (VarPro)

Description

Identifies predictors associated with a response by comparing forest rule regions with their near-miss sets. Supports regression, multivariate regression, classification, and right-censored survival. The returned object contains the rule-generating forest and rule-level scores; importance summarizes these scores and cv.varpro selects an importance cutoff using predictive performance.

Usage

varpro(formula, data, nvar = 30, ntree = 500,
       split.weight = TRUE, split.weight.method = NULL, sparse = TRUE,
       nodesize = NULL, max.rules.tree = 150, max.tree = min(150, ntree),
       parallel = TRUE, cores = get.mc.cores(),
       verbose = FALSE, seed = NULL, ...)

Arguments

formula	Formula specifying the response and predictors, such as y ~ ., cbind(y1, y2) ~ ., or Surv(time, status) ~ .. Use a factor response for classification and numeric responses for regression. Survival status is coded 0 for censoring and 1 for an event.
data	Data frame containing the response and predictors. Incomplete observations are removed, unused factor levels are dropped, and categorical predictors are hot-encoded. If observations are omitted, a warning reports the input, omitted, and retained counts. See Details.
nvar	Maximum number of processed predictor columns retained by split-weight screening before rule generation. A categorical predictor can contribute several columns. Applies when computed or custom split-weights are used; split.weight = FALSE without custom weights uses all processed predictors.
ntree	Number of trees in the rule-generating forest.
split.weight	Should preliminary split-weights guide rule generation? Positive-weight predictors are screened to at most nvar columns, and candidate split variables are sampled according to their retained weights.
split.weight.method	Character string or vector selecting preliminary weighting methods: "lasso", "tree", and "vimp". Methods can be combined, for example c("lasso", "tree"). The default NULL selects a combination automatically. See Details.
sparse	Should preliminary weighting concentrate more strongly on promising predictors? Useful when relatively few predictors are expected to carry signal. See Details.
nodesize	Minimum terminal node size for the rule-generating forest. Larger values produce broader regions with more observations for the local comparisons. The default is chosen automatically from the data dimensions.
max.rules.tree	Maximum number of tree branches sampled per selected tree. A branch can yield several comparisons, one for each variable released from its rule.
max.tree	Maximum number of trees sampled for rule extraction. Increasing this value provides more trees for importance aggregation, up to the number grown.
parallel	Should the preliminary lasso cross-validation folds run in parallel?
cores	Number of cores for lasso-fold parallelism. The default respects options(mc.cores = ...) or MC_CORES. Forest parallelism is controlled separately by randomForestSRC.
verbose	Should progress messages and preliminary split-weights be printed?
seed	Seed supplied to the main forest calculations. Use set.seed before the call to also control R-level sampling, lasso folds, and auxiliary calculations. See Note.
...	Additional named controls for split-weight construction, rule generation, and external survival estimation. The supported controls used most often are described in Details. Unnamed controls, unrecognized names, and duplicate names cause an error.

Details

What VarPro measures

VarPro studies a predictor's contribution by comparing responses in a rule region and its near-miss set. In regression, the target is the conditional mean of the response; in classification, it is the vector of conditional class probabilities. Survival analysis uses a forest estimate of mortality or restricted mean survival time as its response summary.

The workflow has three steps. varpro prepares the data, generates rules, and stores their release comparisons. importance combines the comparisons into variable-level scores. cv.varpro uses the scores to form candidate variable sets and chooses among them using forest prediction error.

Rules and near-miss comparisons

A root-to-terminal-node path defines a region R by combining its split conditions. Releasing predictor s removes every condition involving that predictor while retaining the conditions on the other predictors. This produces the enlarged release region R^{(s)}. The code compares responses in R with those in the near-miss set

C_s = R^{(s)} \setminus R = R^c \cap R^{(s)}.

For example, a rule x_1 \leq 0.5,\ x_2 > 0.3 becomes x_2 > 0.3 when x_1 is released. Its near-miss set contains observations with x_1 > 0.5,\ x_2 > 0.3. Comparing responses in the original region and the near-miss set measures the local contribution of x_1 while preserving the restriction on x_2.

Rules are extracted from the forest. At most max.tree trees and max.rules.tree branches per selected tree contribute to the analysis. Original-region observations are out-of-bag for the corresponding tree. Repeated comparisons allow a predictor to contribute differently in different regions of the feature space.

Guided rule generation and the lasso

Preliminary split-weights help the forest form informative rules. Three sources of information are available:

"lasso"

Absolute coefficients from a cross-validated penalized regression on standardized predictors. This favors predictors with strong linear effects or categorical contrasts.

"tree"

Split frequencies from a shallow forest, which can identify nonlinear effects and interactions.

"vimp"

Positive permutation importance from a preliminary forest, measuring a predictor's contribution to prediction accuracy.

Combining methods uses complementary evidence when forming the weights. With split.weight.method = NULL, the combination is chosen automatically. Set an explicit method to control which sources are used.

The lasso penalty is selected by cross-validation. The weights use the more regularized solution within one standard error of the minimum error. Coefficient magnitudes are combined across responses or classes when needed. For survival, the lasso uses the response summary described below.

With sparse = TRUE, stronger preliminary weights receive greater emphasis. Use FALSE to spread the emphasis more broadly. At most nvar positive-weight processed columns are retained for rule generation. Increase nvar to explore more predictors, or use split.weight = FALSE to generate rules using all processed predictors without preliminary weighting.

The retained weights guide candidate-variable sampling in the forest. Final VarPro importance is calculated from the resulting release comparisons.

Importance scores and selection

Call importance(o) to summarize a varpro object. Its default local.std = TRUE recomputes locally standardized comparisons from the region memberships. For regression, the score comparing R and its near-miss set C_s is the absolute Welch statistic

T = \frac{|\bar y_R - \bar y_{C_s}|} {\sqrt{s_R^2/n_R + s_{C_s}^2/n_{C_s}}}.

The same calculation is applied separately to each numeric response, including external survival summaries. For classification, the overall comparison is the square root of a Pearson chi-squared statistic divided by its degrees of freedom, with continuity correction for two-class tables. Class-specific comparisons use absolute standardized differences in proportions.

Within each tree, comparisons for a predictor are averaged with weights equal to their original-region sample sizes. Class-specific summaries use the corresponding class counts. The tree-level scores are then averaged after winsorization at the lower and upper 10 percent quantiles by default. With local standardization, this mean is the reported z score. With local.std = FALSE, importance summarizes the stored rule scores and divides their winsorized mean by their winsorized standard deviation across trees.

Larger scores indicate stronger response changes under release. The default importance cutoff of 0.79 is used to highlight variables in its plots. To choose a cutoff from predictive performance, use cv.varpro. Its optional cv.folds argument assesses the complete selection procedure on held-out folds and appends stability summaries to the returned importance tables. See importance.varpro for the summary and plotting controls.

Categorical predictors and outcome families

Numeric predictors retain their values. Two-level categorical predictors become binary columns; predictors with more levels become indicator columns. Screening, forest rules, and importance(o) operate on these processed columns. o$x stores the complete processed predictor data, while o$xvar.names identifies the columns retained for rule generation.

get.orgvimp(o) summarizes scores at the original-variable level using the largest score among a predictor's encoded columns. For multivariate regression it also takes the maximum across outcomes. cv.varpro uses this original-variable summary for its candidate sets. get.vimp(o) returns named scores for the represented processed columns. With pretty = FALSE, get.vimp and get.orgvimp return all processed or original predictors, respectively, filling absent scores with zero.

In binary classification, the working response labels are 0 for the majority class and 1 for the minority class. For imbalanced classes, the default handling uses AUC splitting and RFQ with geometric-mean performance for the preliminary permutation-importance calculation. The original response labels remain in y.org and the rule-generating forest. The controls use.rfq and iratio.threshold modify this behavior.

Survival targets

A preliminary survival forest estimates a numeric response summary for each training observation. By default, this is forest mortality. Supplying rmst = tau requests restricted mean survival time, obtained by integrating the estimated survival curve up to tau. A vector of horizons supplies one numeric response per horizon.

These full-ensemble training summaries are used for preliminary weighting and the default locally standardized importance. With mortality or one RMST horizon, rules are generated using the original survival outcome. With multiple RMST horizons, rule generation uses multivariate regression on the RMST responses. Thus the horizons determine the survival features being prioritized. The controls ntree.external, nodesize.external, and ntime.external govern the preliminary survival forest.

Additional controls and computation

Useful arguments supplied through ... include nfolds for lasso cross-validation, rmst for the survival target, and ntree.external for the preliminary survival forest.

Custom split-weights can be supplied as a named numeric vector through split.weight.custom. Names must refer to processed predictor columns; omitted columns receive zero weight. Use get.splitweight.custom(formula, data) to obtain a correctly named starting vector, then assign positive weights to the predictors of interest.

Reducing ntree, max.tree, or max.rules.tree reduces computation. Increasing nodesize produces larger comparison groups with less local detail. The parallel and cores arguments control lasso-fold parallelism; forest threading is controlled by randomForestSRC.

Value

An object of class "varpro" with the following components:

rf	The rule-generating rfsrc forest.
split.weight	Named weights for the retained predictor columns, or NULL when computed and custom weighting are both absent.
split.weight.raw	List of untransformed preliminary components actually calculated, named lasso, tree, or vimp. Empty when preliminary estimation is skipped.
results	Rule-level data frame with tree and branch identifiers, released-variable index, original-region sample size n.oob, and importance columns. Indices refer to xvar.names. Classification includes class-specific counts and scores; multivariate regression includes one score per response. importance supplies the variable-level summary.
x	Complete-case processed predictor data, including columns removed by preliminary screening. Encoding information is stored as attributes.
xvar.names	Processed predictor names retained for rule generation.
xvar.org.names	Original predictor names before encoding.
y	Working response: numeric response(s), recoded class labels, or external survival summary values.
y.org	Original response values for the retained observations. For survival, these are the observed times and event indicators.
yvar.names	Names of the working response columns.
family	Family of the rule-generating forest: "regr", "regr+", "class", or "surv". Multiple RMST horizons use "regr+".
max.rules.tree, max.tree	Stored rule-extraction limits used by subsequent importance calculations.
model.info	Compact description of the analysis. Its observations component records input, retained, and omitted row counts for the data supplied to this varpro call.

Note

For reproducible serial analyses, call set.seed immediately before varpro and use parallel = FALSE. The seed argument controls only the forest calls to which it is passed. Rule resampling and lasso folds also use random numbers. Calling importance again can resample rules; retain a summary for repeated inspection of the same calculation.

The ... controls are consumed by specific stages. General rfsrc arguments such as na.action, splitrule, and perf.type are not forwarded wholesale. Perform any desired imputation before calling varpro.

Author(s)

Min Lu and Hemant Ishwaran

References

Lu, M. and Ishwaran, H. (2024). Model-independent variable selection via the rule-based variable priority. arXiv:2409.09003. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.48550/arXiv.2409.09003")}.

Friedman, J., Hastie, T. and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1), 1–22. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v033.i01")}.

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

Ishwaran, H. (2025). Multivariate Statistics: Classical Foundations and Modern Machine Learning. Chapman and Hall/CRC.

See Also

importance.varpro, cv.varpro, predict.varpro, ivarpro, uvarpro, partialpro, outpro, isopro, randomForestSRC::rfsrc, glmnet::cv.glmnet.

Examples

## Small regression example.
o <- varpro(mpg ~ ., mtcars, ntree = 1)


## Nonlinear regression with five signal and five noise predictors.
set.seed(137)
n <- 300
x <- matrix(runif(n * 10), nrow = n)
colnames(x) <- sprintf("x%02d", seq_len(ncol(x)))
mu <- 10 * sin(pi * x[, 1] * x[, 2]) +
      20 * (x[, 3] - 0.5)^2 + 10 * x[, 4] + 5 * x[, 5]
d <- data.frame(y = mu + rnorm(n), x)
signal <- colnames(x)[1:5]

## Shallow-tree-guided rule generation.
o <- varpro(y ~ ., d, ntree = 25, ntree.reduce = 25,
            split.weight.method = "tree", nodesize = 10,
            max.tree = 25, max.rules.tree = 20, parallel = FALSE)
imp <- importance(o)
print(data.frame(variable = rownames(imp), z = imp$z,
                 signal = rownames(imp) %in% signal))

## Automatic weighting with more trees and comparisons.
set.seed(137)
o <- varpro(y ~ ., d, ntree = 150, max.tree = 100,
            max.rules.tree = 50, nfolds = 5, parallel = FALSE)
imp <- importance(o, plot.it = TRUE)
print(imp)
print(get.vimp(o))
print(get.vimp(o, pretty = FALSE))
print(get.topvars(o))

## Choose a variable set using predictive performance.
set.seed(137)
## Add cv.folds = 5 to assess selection and obtain stability summaries.
cv <- cv.varpro(y ~ ., d, ntree = 100,
                zcut = seq(0.1, 3, length.out = 12), nblocks = 6,
                max.tree = 75, max.rules.tree = 40,
                nfolds = 5, parallel = FALSE)
print(cv$imp)
print(cv$imp.conserve)

## Hot-encoding: compare processed-column and original-variable scores.
## Several predictors are converted to factors before analysis.
data(BostonHousing, package = "mlbench")
Boston <- BostonHousing
Boston$zn <- factor(Boston$zn)
Boston$chas <- factor(Boston$chas)
Boston$lstat <- factor(round(0.2 * Boston$lstat))
Boston$nox <- factor(round(20 * Boston$nox))
Boston$rm <- factor(round(Boston$rm))
set.seed(137)
b <- varpro(medv ~ ., Boston, ntree = 100,
            max.tree = 75, max.rules.tree = 40,
            nfolds = 5, parallel = FALSE)

## Calculate the processed-column importance
set.seed(137)
b.imp <- importance(b)
print(b.imp)

## Map max scores back to original predictors.
set.seed(137)
b.org <- get.orgvimp(b)
print(b.org)
print(get.orgvimp(b, vmp = b.imp)) ## should be identical; no seed required

## Include every original predictor, filling absent scores with zero.
set.seed(137)
print(get.orgvimp(b, pretty = FALSE))


## Multiclass classification: overall and class-specific importance.
data(wine, package = "randomForestSRC")
wine$quality <- factor(wine$quality)
set.seed(137)
w <- varpro(quality ~ ., wine, ntree = 100,
            max.tree = 75, max.rules.tree = 40,
            nfolds = 5, parallel = FALSE)
w.imp <- importance(w)
print(w.imp$unconditional)
print(w.imp$conditional.z)

## Survival: canonical example
data(peakVO2, package = "randomForestSRC")
s <- varpro(Surv(ttodead, died)~., peakVO2)
print(importance(s))


## Survival: prioritize variables affecting RMST through 500 days.
data(pbc, package = "randomForestSRC")
pbc <- na.omit(pbc)
set.seed(137)
s <- varpro(Surv(days, status) ~ ., pbc, rmst = 500,
            ntree = 100, max.tree = 75, max.rules.tree = 40,
            nfolds = 5, parallel = FALSE)
print(importance(s))
print(get.orgvimp(s))
print(get.vimp(s, pretty = FALSE))

## Multiple horizons give one importance summary per horizon.
set.seed(137)
s.multi <- varpro(Surv(days, status) ~ ., pbc,
                  rmst = c(500, 1000), ntree = 100,
                  max.tree = 75, max.rules.tree = 40,
                  nfolds = 5, parallel = FALSE)
print(importance(s.multi))

varPro documentation built on Sept. 22, 2026, 5:09 p.m.