quantreg.rfsrc: Quantile Regression Forests

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

Quantile Regression Forests

Description

Estimates conditional response distributions and their quantiles using random forests. The main method estimates a distribution of training residuals and shifts it by the forest mean prediction. The default method = "exchangeable" gives residuals equal weight; method = "forest" uses forest weights. An optional Greenwald–Khanna approximation is available for large data sets. Supports continuous responses in univariate, multivariate, and mixed-outcome forests.

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 responses and predictors. Required when growing a new forest. Converted to a plain data frame.
object	The original training object returned by quantreg, including its saved forest, training responses, and residuals.
newdata	Data for new observations, containing the same predictor variables as the training data. Converted to a plain data frame. Response columns are optional and are used for performance evaluation. If omitted or NULL, calculate results for the training observations.
method	The residual method adds residual quantiles to the forest mean prediction. Choose how the residuals are weighted: "exchangeable" Equal weights (default). "forest" Forest weights that depend on the predictor values. For large data sets, "gk" uses the Greenwald–Khanna quantile approximation instead of the residual method. It also accepts "GK", "G-K", and "g-k". The former name "local" is an alias for "exchangeable". To reuse a nondefault method, specify it in each call to quantreg(object = fit, ...).
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 rule when either "la.quantile.regr" or "quantile.regr" is specified. Other applicable rules include "mse" for regression and "mahalanobis" for multivariate regression.
prob	Quantile levels to calculate, as a nonempty numeric vector of finite probabilities strictly between zero and one. Can be supplied during fitting or prediction; levels are sorted and duplicates retained. With NULL, the residual methods calculate percentiles 1 through 99 when fitting. With object supplied, prob = NULL reuses the saved levels. GK training defaults also depend on prob.epsilon; supply prob to control the number of quantiles.
prob.epsilon	Greenwald–Khanna approximation tolerance on the probability scale, a finite number in (0,1). Larger values allow coarser summaries. With object supplied, NULL reuses the saved tolerance. When both prob and prob.epsilon are NULL during GK fitting, the default requests 2n quantiles, where n is the processed training sample size.
oob	Use out-of-bag (OOB) predictions for training observations if TRUE (default), or predictions from all trees if FALSE. When fitting, this also selects the training residuals for the residual methods. Calls with object reuse the saved residual choice. Predictions for newdata always use all trees. See OOB predictions in Details.
fast	If TRUE, grow a new forest with rfsrc.fast instead of rfsrc. Ignored when object is supplied.
maxn	Maximum number of response values in the CDF reporting grid, as a positive integer or Inf. Controls the resolution of cdf, density, and grid-integrated CRPS. For a fixed forest and prediction call, it does not change the quantiles or the residual-method means and standard deviations.
o	A training or prediction object returned by quantreg. CRPS and pinball scoring require observed responses in o$yvar.
target.prob	Quantile levels to extract with get.quantile. NULL returns all stored levels. Otherwise, supply finite probabilities in [0,1]; each is matched to the nearest stored level.
pretty	If TRUE, simplify a helper result for one continuous response to a matrix, data frame, or numeric vector, as described in Value. If FALSE, return a response-named list.
subset	Rows to evaluate for CRPS or pinball loss. Supply positive integer indices or a logical vector with one nonmissing entry per row of the returned object. NULL uses all rows. The selection must be nonempty; repeated indices include a row more than once.
standardize	If TRUE, divide each cumulative CRPS integral by its integration width.
tau	Quantile levels at which to calculate pinball loss, as finite probabilities strictly between zero and one. Scores follow the supplied order, including repeated levels. NULL uses the session option quantreg.tau, with the defaults described under Evaluating predictions. Supplied 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. Extracts results from o without recalculating predictions.
...	Additional options for rfsrc when growing a forest, or for its prediction method when object is supplied. quantreg sets the quantile method and any required forest-weight options.

Details

quantreg() estimates the conditional cumulative distribution function (CDF) of a response Y given predictor values X=x,

F(s\mid x)=P(Y\leq s\mid X=x),

and obtains conditional quantiles from this distribution. This is different than nonparametric regression which targets the conditional mean E(Y\mid X=x). Intervals formed from the conditional quantiles describe individual response values; they are not confidence intervals for the conditional mean.

The residual method estimates this conditional distribution using forest mean predictions and the empirical distribution of training residuals. For training observation j, the residual is

R_j = Y_j - \widehat m_j,

where \widehat m_j is the training prediction selected when fitting the forest. For an observation with predictor values x, adding each residual to the forest mean prediction \widehat m(x) gives adjusted response values \widehat m(x)+R_j. Their weighted empirical distribution estimates F(s\mid x).

Let \mathcal I_R index the finite training residuals, and let w_j(x) be their nonnegative weights, summing to one. The residual CDF is

\widehat F_R(u\mid x) = \sum_{j\in\mathcal I_R} w_j(x) I\{R_j\leq u\}.

Its quantile \widehat Q_R(\tau\mid x) is the smallest residual whose cumulative weight is at least \tau. Adding the forest mean prediction gives the response quantile and shifts the CDF:

\widehat Q(\tau\mid x) = \widehat m(x) + \widehat Q_R(\tau\mid x),\qquad \widehat F(s\mid x) = \widehat F_R\{s-\widehat m(x)\mid x\}.

Predicted quantiles can fall outside the range of the training responses.

The default method = "exchangeable" uses w_j(x)=1/N_R, where N_R is the number of finite training residuals. This uses a common residual distribution for all predictor values. The estimated conditional response distributions therefore have the same shape and spread, with location determined by \widehat m(x). With the default OOB residuals, prediction intervals follow the signed OOB-error construction in Zhang et al. (2020, Section 2.3). Equal weights give the same residual quantiles as quantile(..., type = 1). When calculating quantiles for a training observation, its own finite residual remains in the pool.

With method = "forest", the same residuals receive forest weights that depend on x. This allows both the shape and spread of the predicted distribution to vary with the predictors.

Meinshausen (2006) estimates the conditional CDF by weighting the training responses directly:

\widehat F_{\mathrm{QRF}}(s\mid x) = \sum_{j=1}^n v_j(x) I\{Y_j\leq s\},

where v_j(x) are forest weights. The "forest" option here instead weights residuals and then shifts their distribution by \widehat m(x). It assigns probability to the adjusted responses Y_j-\widehat m_j+\widehat m(x). The two estimators generally differ even if the same weights are used.

OOB predictions

When fitting a forest, oob determines which trees supply the training predictions used to calculate residuals. With oob = TRUE (default), \widehat m_j averages predictions from trees whose training samples exclude observation j. With oob = FALSE, it uses all trees. The fitted forest retains this choice of training residuals.

When the original fitted object is supplied as object, oob selects OOB or full-ensemble mean predictions, forest weights, or GK quantiles for the training observations. It does not change the saved residuals. For example, quantreg(object = fit, oob = FALSE) uses all trees for the training predictions while keeping the residuals selected when fit was created.

For newdata, predictions and weights always use all trees. The residual methods reuse the saved training residuals, so a default fit uses OOB residuals for new observations as well. When fitting and prediction are combined in one call, oob selects the training residuals before predicting the new data.

If an entire OOB prediction vector, weight matrix, or GK quantile matrix is unavailable (NULL), the corresponding result from all trees is used. Individual missing entries are not replaced.

Greenwald–Khanna method for large data

For large data sets, method = "gk" uses the Greenwald–Khanna (2001) algorithm for compact approximate quantile summaries. This option obtains quantiles from the forest without using training residuals or requesting a forest-weight matrix. Specify prob to control the number of returned quantiles. The argument prob.epsilon controls the approximation tolerance on the probability scale; larger values allow coarser summaries.

The CDF is reconstructed from the returned quantiles. At each grid value, it takes the largest requested probability whose quantile is no greater than that value. The CDF is zero below all requested quantiles and reaches one at the largest finite training response. Thus prob controls the quantiles and the resolution of the reconstructed CDF, while maxn controls where that CDF is evaluated. GK means and standard deviations use the CDF increments and depend on the grid and the probability mass it represents.

Fitting and results

Supply formula and data to fit a forest, and save the result as fit. The argument prob specifies the quantile levels; for example, prob = c(.05, .50, .95) requests the 5th percentile, median, and 95th percentile. Use quantreg(object = fit, newdata = testdata) to predict new observations, or quantreg(object = fit, prob = c(.1, .5, .9)) to calculate different quantiles for the training observations. These calls reuse the fitted trees and require the original training object.

To fit and predict in one call, supply formula, data, and newdata together. New data can contain predictors alone. Include response columns to calculate performance scores or use plot.quantreg(). These responses are used only for evaluation; residual distributions and the CDF reporting grid use training data.

The residual methods calculate the mean and standard deviation of the predicted distribution from the same residuals and weights. With \bar R(x)=\sum_{j\in\mathcal I_R} w_j(x)R_j, they are

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

The exchangeable method uses equal weights. Residuals are used without recentering, so the distribution mean can differ from the forest mean prediction.

To report the CDF, the function constructs a grid yunq from the sorted unique finite training responses, separately for each continuous response. When there are more than maxn values, it selects approximately equally spaced positions in this ordered set. Use maxn = Inf to include all values.

For the residual methods, cdf evaluates the predicted CDF at these grid values. Its last value may be less than one when the shifted residual distribution extends above the grid. Quantiles, means, and standard deviations use the full residual distribution, including values outside the grid, and do not depend on maxn. The component density contains interval probabilities: its first column is the CDF at the first grid value, and subsequent columns are successive CDF differences. These probabilities are not divided by interval widths. Their row sum equals the last reported CDF value, without adjustment for probability above the grid.

Use get.quantile(o) to extract quantiles and get.quantile.stat(o) to obtain the mean, median, and standard deviation for each observation. In get.quantile(), each requested target.prob is matched to the nearest stored probability, while column labels show the requested levels. To calculate new levels, use quantreg(object = fit, prob = ...) with the original training object. The median returned by get.quantile.stat() uses the stored level nearest 0.5; include 0.5 in prob to obtain that exact level. Use extract.quantile(o) to obtain all quantile components in a named list, one entry per continuous response.

With method = "forest", the weight matrix has one row per evaluated observation and one column per training observation, which can require substantial memory. Reducing maxn reduces the cdf and density matrices; it does not reduce the forest-weight matrix or the number of GK quantiles.

Evaluating predictions

get.quantile.crps(o) compares the predicted CDF with observed responses using the continuous ranked probability score (CRPS). At each reporting-grid value s_k, it averages [I\{Y_i\leq s_k\}-\widehat F(s_k\mid X_i)]^2 over the selected observations with finite response and CDF values. It integrates these errors by the trapezoidal rule from the first grid value to s_k, returning a curve with columns y and crps. Smaller scores indicate less prediction error.

The last entry integrates over the full grid range, so the score depends on both the grid range and resolution. With standardize = TRUE, each integral is divided by its integration width. The first entry is then NA because its width is zero; with standardize = FALSE, the first entry is zero.

get.pinball.error(o) evaluates individual quantiles. For response Y_i and predicted quantile \widehat Q_i(\tau), the pinball loss is

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

The helper returns the mean loss at each requested level, using the selected observations with finite responses and quantiles. Smaller values indicate better quantile predictions. If no observations are available at a level, its mean loss is NA.

For the residual methods, pinball scoring uses stored quantiles at exactly matching levels, or calculates them from the saved residual distribution. Levels in tau need not appear in prob, and the score does not depend on the CDF reporting grid. GK scoring interpolates the reported CDF. Scoring uses the predictions and residual choice saved in o, without refitting the forest or recalculating its mean predictions.

For a univariate quantile object with observed responses and performance output, print.rfsrc reports raw and standardized grid-integrated CRPS and mean pinball losses. For a multivariate object, use print(object, outcome.target = "response.name") to display quantile scores for one continuous response. Without outcome.target, the display gives mean and response-specific regression or classification errors. Selecting a response for printing or extraction does not recalculate predictions. OOB labels refer to the predictions used in each calculation.

Set pinball levels for one print call with print(object, quantreg.tau = c(.2, .5, .8)), or for the session with options(quantreg.tau = c(.2, .5, .8)). Explicit levels override the session option. Otherwise the defaults are taken, in order, from quantreg.tau, the former option rfsrc.pinball.taus, or c(.1, .5, .9). These defaults also apply to get.pinball.error(). Invalid levels cause an error. These settings control reporting without modifying the fitted object.

To display pinball losses on a plot, use plot.quantreg(object, quantreg.tau = c(.2, .5, .8)). The losses appear in the lower-right corner of the main panel; the quantile interval and CRPS inset are unchanged. Pinball labels are omitted unless quantreg.tau is supplied in the plotting call, even when the session option is set.

Value

quantreg() returns an rfsrc training (grow) or prediction object with the additional class quantreg. For a univariate regression response, o$quantreg contains:

quantiles	Matrix with evaluated observations in rows and requested probability levels in columns, including when there is only one observation or level.
prob	Sorted probability levels corresponding to quantiles.
cdf	CDF matrix with evaluated observations in rows and yunq grid values in columns.
density	CDF increments, with the same dimensions as cdf.
yunq	Ordered training-response grid used to report the CDF.
mean, std	Mean and standard deviation of the estimated conditional response distribution, as vectors for the residual methods. The standard deviation describes response variability, not the standard error of a mean estimate. NULL for GK; use get.quantile.stat() for its grid-based summaries.
method	Method used in this call: "exchangeable", "forest", or "gk". The former name "local" is reported as "exchangeable".
oob	Whether the reported mean predictions or GK quantiles use OOB output.
prediction	Forest mean predictions used to shift the residual distribution. NULL for GK.
residual.source	Source of the saved training residuals: "oob" or "ensemble" for the residual methods. NULL for GK. A new-data result can have residual.source = "oob" and oob = FALSE: the training residuals are OOB and the new predictions use all trees.

For multivariate and mixed-outcome forests, o$quantreg is a named list of these components, one per returned continuous response. Classification responses have no quantile component.

Training data and residuals

o$xvar contains predictors for the observations being predicted; o$yvar contains their responses, if supplied. The saved forest's xvar and yvar contain the training data. The CDF grid and null CRPS reference use training responses; performance scores use responses for the observations being evaluated.

For a univariate training object, o$residual = o$yvar - o$predicted and o$residual.oob = o$yvar - o$predicted.oob, when the respective predictions are available. For multivariate and mixed outcomes, these vectors are in o$regrOutput[[response.name]]. If predicted or predicted.oob is NULL, its residual vector is also NULL. Missing OOB predictions give missing residuals.

Both sets of training residuals are saved in o$forest$residual and o$forest$residual.oob. They are vectors for univariate regression and matrices for multivariate or mixed outcomes, with training observations in rows and named continuous responses in columns. If residuals are available for only some responses, the remaining columns contain NA. If none are available, the residual matrix is NULL. o$forest$quantreg$residual.source records which set was selected for each response. Prediction objects retain these training residuals under forest.

Nonfinite residuals are excluded together with their weights, and the remaining forest weights are normalized for each observation. Quantiles and summaries are NA if no finite residuals remain, the forest mean prediction is nonfinite, or the remaining weights include nonfinite values or their sum is not positive. Negative weights and incompatible weight-matrix dimensions cause an error.

Helper functions

extract.quantile() always returns a response-named list of quantile components. 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 the 2001 ACM SIGMOD International Conference on Management of Data, 58–66. https://sanjeevkhanna.org/papers/sigmod01_quantiles.pdf.

Meinshausen N. (2006). Quantile regression forests. Journal of Machine Learning Research, 7:983–999. https://jmlr.org/papers/v7/meinshausen06a.html.

Zhang H., Zimmerman J., Nettleton D. and Nordman D.J. (2020). Random forest prediction intervals. The American Statistician, 74(4):392–406. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/00031305.2019.1585288")}.

See Also

rfsrc, predict.rfsrc, plot.quantreg

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)

## Check the stored levels and probability covered by the CDF grid.
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 prediction for new observations
## ------------------------------------------------------------
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 reuses the saved OOB residuals.
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 observations
print(c(evaluated.rows = o.exchangeable$n,
        training.rows = length(o.exchangeable$forest$residual.oob)))

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

## ------------------------------------------------------------
## Changing maxn leaves residual quantiles and summaries unchanged
## ------------------------------------------------------------
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"])))

## Recalculate all responses; select Temp 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)))

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

## Impute missing values using random splitting.
iowa <- housing
iowa$PID <- NULL
iowa$SalePrice <- log(iowa$SalePrice)
iowa <- impute(SalePrice ~. , iowa, splitrule = "random", nimpute = 1)

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

randomForestSRC documentation built on Nov. 1, 2026, 5:07 p.m.