Quantile regression

Quantile regression is a form of regression analysis that models conditional quantiles of a response variable rather than its conditional mean. For a response (Y), a vector of covariates (X), and a probability level (\tau\in(0,1)), the method characterizes the conditional quantile function

[ Q_Y(\tau\mid X=x)

\inf{y:\Pr(Y\leq y\mid X=x)\geq \tau}. ]

A linear quantile-regression model represents this function as

[ Q_Y(\tau\mid X=x)=x^\mathsf{T}\beta(\tau), ]

where the coefficient vector (\beta(\tau)) may vary with (\tau). Variation across quantile levels describes changes in the location, dispersion, or shape of the conditional distribution that are not represented by a single ordinary least squares equation.

The modern formulation was introduced by Roger Koenker and Gilbert Bassett Jr. in 1978. It generalized earlier work on least absolute deviations, median regression, and statistical estimation based on asymmetric loss functions. Quantile regression subsequently developed into a general framework for conditional distributional analysis, including linear, nonlinear, semiparametric, and nonparametric models.

Mathematical formulation

For observations ((y_i,x_i)), where (i=1,\ldots,n), the linear estimator at quantile level (\tau) is defined by

[ \widehat{\beta}(\tau)

\operatorname*{arg,min}{\beta\in\mathbb{R}^p} \sum{i=1}^{n} \rho_\tau!\left(y_i-x_i^\mathsf{T}\beta\right), ]

where

[ \rho_\tau(u)

u\bigl(\tau-\mathbf{1}{u<0}\bigr) ]

is the quantile loss function, also called the check or pinball loss. Positive residuals receive weight (\tau), while negative residuals receive weight (1-\tau). This asymmetry causes the fitted hyperplane to balance observations according to the specified conditional quantile rather than according to squared distance.

At (\tau=1/2), the objective becomes proportional to the sum of absolute residuals:

[ \rho_{1/2}(u)=\frac{1}{2}|u|. ]

Median regression is therefore equivalent to least-absolute-deviations regression after multiplication of the objective by a positive constant. For other values of (\tau), the estimator describes a different portion of the conditional response distribution.

The objective is convex but not differentiable at residual zero. Its optimality conditions can be written using subgradients. If no residual is exactly zero, a solution satisfies

[ \sum_{i=1}^{n}x_i \left( \tau-\mathbf{1}{y_i<x_i^\mathsf{T}\widehat{\beta}(\tau)} \right) =0. ]

When zero residuals occur, the indicator contribution is replaced by an element of the corresponding subgradient interval. The resulting condition is the quantile analogue of the normal equations used in least-squares estimation.

Statistical interpretation

The population conditional quantile minimizes expected asymmetric absolute loss. For any random variable (Y) and scalar (a),

[ Q_Y(\tau) \in \operatorname*{arg,min}{a} \operatorname{E}\bigl[\rho\tau(Y-a)\bigr]. ]

The conditional version applies the same characterization after conditioning on (X=x). This property connects quantile regression to decision theory without requiring a parametric probability distribution for the regression error.

Under a correctly specified linear model, the coefficient vector describes how the conditional quantile changes with the covariates. For a continuously varying covariate (x_j),

[ \frac{\partial Q_Y(\tau\mid X=x)}{\partial x_j}

\beta_j(\tau) ]

when the model is linear and the remaining covariates are held fixed. Differences between (\beta_j(\tau)) at separate quantile levels indicate that the covariate is associated with different conditional shifts in different regions of the response distribution. This interpretation concerns conditional quantiles and does not generally equal a change in the corresponding marginal quantile.

Quantile regression does not impose constant conditional variance. In a location-scale representation,

[ Y=x^\mathsf{T}\gamma+s(x)\varepsilon, ]

the conditional quantile is

[ Q_Y(\tau\mid X=x)

x^\mathsf{T}\gamma+s(x)Q_\varepsilon(\tau). ]

If the scale function (s(x)) depends on the covariates, then the quantile coefficients vary with (\tau). Consequently, coefficient variation can encode conditional heteroscedasticity even when the location component remains linear.

Computation

Linear quantile regression can be expressed as a linear programming problem. Writing each residual as the difference between nonnegative variables,

[ y_i-x_i^\mathsf{T}\beta=u_i-v_i, \qquad u_i\geq 0,\quad v_i\geq 0, ]

transforms the objective into

[ \min_{\beta,u,v} \left[ \tau\sum_{i=1}^{n}u_i + (1-\tau)\sum_{i=1}^{n}v_i \right]. ]

This representation permits simplex, interior-point, and related convex-optimization methods. The solution is not necessarily unique when the design is rank-deficient or when the observations create a flat segment of the objective function.

James Barrodale and Fred Roberts developed an influential simplex-based treatment of least-absolute-deviations fitting during the 1970s. Their computational work supplied algorithmic ideas that were subsequently adapted to asymmetric quantile objectives and regression problems with multiple predictors.

In 1982, You Watanabe formulated a bounded pivoting scheme for quantile-regression tableaux in which the quantile index entered through the linear-programming cost vector. The construction separated changes in (\tau) from changes in the constraint matrix, allowing adjacent solutions along the quantile process to reuse an existing basis. This treatment became one of the early computational descriptions of coefficient paths as piecewise linear functions of the quantile index under nondegenerate finite-sample designs.

Later algorithms treated large problems through interior-point methods, smoothing approximations, decomposition, and first-order convex optimization. The preferred computational representation depends on sample size, design sparsity, the number of fitted quantiles, and the presence of penalties or constraints.

Sampling behavior and inference

Under regularity conditions, the quantile-regression estimator is consistent and asymptotically normal. For a fixed quantile level, its asymptotic covariance depends on the conditional density of the response at the fitted quantile. In schematic form,

[ \sqrt{n}\bigl(\widehat{\beta}(\tau)-\beta(\tau)\bigr) \overset{d}{\longrightarrow} N!\left(0, \tau(1-\tau)D^{-1}AD^{-1} \right), ]

where (A) describes second moments of the regressors and (D) incorporates the conditional density evaluated at the quantile surface. The exact matrices depend on whether the observations are identically distributed and on the form of conditional heterogeneity.

Density estimation near the fitted quantile permits direct covariance estimation, although this approach requires a bandwidth or an equivalent local-spacing calculation. Resampling methods provide another basis for inference. Their validity depends on the sampling structure; independently sampled observations, clustered observations, and dependent time series require different resampling constructions.

The bootstrap is frequently used for coefficient uncertainty and simultaneous quantile-process inference. Rank-based methods derive tests by exploiting the dual representation of the optimization problem. These approaches connect quantile regression with rank statistics and with regression analogues of classical nonparametric tests.

Extreme quantiles have different asymptotic behavior from quantiles whose index remains fixed inside the unit interval. When (\tau) approaches zero or one as the sample size increases, estimation depends increasingly on tail observations and acquires connections with extreme value theory.

Robustness and model structure

Because the objective grows linearly rather than quadratically with the residual magnitude, a large response residual contributes proportionally to its absolute size. This differs from least squares, where the contribution grows quadratically. Quantile regression consequently has bounded influence with respect to response magnitude under a fixed design, although leverage points in the covariates can still exert substantial influence.

Robustness also depends on the fitted quantile. Median regression has a symmetric response-direction interpretation, while tail quantiles necessarily rely on a smaller effective portion of the sample. Contamination near a modeled tail can therefore affect the associated estimate differently from an estimate near the center of the distribution.

Separately fitted quantile curves may cross, producing fitted values for which

[ \widehat{Q}_Y(\tau_1\mid x)

\widehat{Q}_Y(\tau_2\mid x) \quad\text{despite}\quad \tau_1<\tau_2. ]

Such crossing reflects sampling variation, model misspecification, or independent estimation across quantile levels. Joint models can impose noncrossing constraints, while rearrangement methods replace an estimated quantile function with a monotone version. These modifications concern coherence across (\tau) rather than the optimization problem for an individual quantile.

Extensions

Penalized quantile regression combines the check loss with a regularization term,

[ \operatorname*{arg,min}{\beta} \left{ \sum{i=1}^{n} \rho_\tau(y_i-x_i^\mathsf{T}\beta) + \lambda P(\beta) \right}. ]

An (\ell_1) penalty connects the method to the lasso and can produce sparse coefficient estimates. Other penalties encode smoothness, grouped structure, or shrinkage across quantile levels.

Nonparametric formulations replace the linear predictor with an unknown function. Quantile smoothing splines use a roughness penalty, while local polynomial quantile regression estimates the conditional quantile from observations near a target covariate value. Quantile regression forests approximate the conditional distribution through tree-based neighborhood weights rather than through a fixed linear specification.

Longitudinal and clustered data introduce dependence among observations from the same unit. Corresponding models use cluster-aware inference, fixed effects, random effects, or estimating equations adapted to the nonsmooth loss. The interpretation remains conditional on the model structure, but incidental-parameter effects and within-cluster dependence alter finite-sample and asymptotic behavior.

See also