quantreg.rfsrc: Quantile Regression Forests

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

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

formula	Model formula. Required when growing a new forest and ignored when object is supplied. At least one response must be continuous.
data	Training data containing the response and predictor variables. Required when growing a new forest. Input is converted to a plain data frame before formula processing.
object	The original quantreg grow object, including its saved forest, training responses, and residuals. Retain this object for subsequent prediction calls.
newdata	Optional data containing the training predictors. Input is converted to a plain data frame. Response columns can also be supplied for performance evaluation. Use quantreg, rather than a direct predict call, to obtain quantile output for new observations. Omitted or NULL newdata restores predictions for the training observations.
method	Quantile calculation method. The default "exchangeable" treats the training residuals as exchangeable, pools their finite values with equal weights, and adds quantiles of this common residual distribution to the forest prediction. With the default oob = TRUE during growth, this is an OOB exchangeable residual distribution. "forest" uses a forest-weighted residual distribution for each evaluated row. "gk" requests native Greenwald–Khanna quantiles; "GK", "G-K", and "g-k" are also accepted. The former name "local" is retained as an alias for "exchangeable". The method applies to the current call. Specify it again in later calls to retain a nondefault method.
splitrule	Splitting rule used to grow the forest. The default for univariate regression is "la.quantile.regr", local adaptive quantile regression splitting. Multivariate forests use their default splitting rule when either "la.quantile.regr" or "quantile.regr" is specified. Other applicable rules can be supplied, including "mse" for regression and "mahalanobis" for multivariate regression.
prob	Nonempty numeric vector of finite probabilities strictly between zero and one. Levels are sorted; duplicates are retained. Can be supplied at training or prediction time. With NULL, the residual methods use percentiles 1 through 99 during training, and an object call uses the saved training probabilities. For GK, the training defaults also depend on prob.epsilon; specify prob explicitly to control the number of native quantiles.
prob.epsilon	Finite probability approximation tolerance in (0,1) for Greenwald–Khanna. Larger values permit coarser summaries. A prediction-time NULL uses the saved training tolerance. When both probability arguments are NULL during GK growth, the default requests 2n quantiles, where n is the processed training sample size.
oob	Selects the residual distribution and the current training-row predictions. During growth, TRUE selects residuals computed as the observed response minus its OOB mean prediction, when that component is available. These residuals define the default exchangeable distribution and are also used by the forest method. FALSE during growth selects full-ensemble residuals instead. Both numerical quantities are retained as residual and residual.oob when their predictions are available. For training-row output, oob also selects the current mean predictions, forest weights, or native GK quantiles. Changing it in a later object call leaves the original residual selection unchanged. Supplied newdata uses ordinary full-ensemble predictions together with the selected training residuals. See the OOB residual discussion below.
fast	Use rfsrc.fast instead of rfsrc when growing a new forest? Does not regrow an existing object.
maxn	Positive integer, or Inf, limiting the number of ordered, unique training-response values in the CDF reporting grid. Controls the resolution of cdf, density, and the grid-integrated CRPS curve. For a fixed forest and prediction call, it does not change the returned quantiles or the residual-method means and standard deviations.
o	A grow or prediction object returned by quantreg. CRPS and pinball scoring also require observed responses in o$yvar.
target.prob	Probabilities to extract with get.quantile. NULL returns all stored levels. Otherwise, supply finite values in [0,1]; each is matched to the nearest stored level.
pretty	Simplify the helper result when there is one continuous response? TRUE returns a matrix, data frame, or numeric vector, as described in Value. FALSE retains a response-named list.
subset	Rows to evaluate for CRPS or pinball loss. Supply positive integer indices or a nonmissing logical vector with one entry per row of the returned object. NULL uses all rows. The selection must be nonempty; repeated indices include a row repeatedly.
standardize	Divide each cumulative CRPS integral by its integration width? The default is TRUE.
tau	Finite pinball-loss levels strictly between zero and one. Order and duplicates are retained. NULL uses the session option quantreg.tau, with the fallbacks described under Printed performance. Explicit levels apply to this call only.
m.target	For get.pinball.error, the name of one continuous response to evaluate. NULL evaluates all available continuous responses. Selects stored results without making a prediction call.
...	Additional forest options. Training options are selected from those recognized by rfsrc; object calls forward options to the prediction method. Options controlled by quantreg, including the quantile method and associated forest weights, take precedence.

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:

quantiles	Matrix with evaluated observations in rows and requested probability levels in columns.
prob	Sorted probability levels corresponding to quantiles.
cdf	CDF matrix evaluated at yunq.
density	Successive CDF increments, with the same dimensions as cdf.
yunq	Ordered training-response reporting grid.
mean, std	Direct distribution mean and standard deviation vectors for the residual methods. NULL for GK, whose moment helper uses the reporting grid.
method	Canonical method name for this call: "exchangeable", "forest", or "gk". Calls using the former "local" name return "exchangeable".
oob	Whether the current mean predictions or native GK quantiles came from the OOB component. This setting does not change the saved residual selection.
prediction	Current mean predictions used to shift the residuals. NULL for GK.
residual.source	The selected training residual distribution: "oob" or "ensemble" for the residual methods. NULL for GK, whose quantiles do not use residuals.

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)

randomForestSRC documentation built on Sept. 16, 2026, 5:06 p.m.