Quantile Regression

NoteSelf-study topic

This reading extends the regression framework to conditional quantiles. It is supplementary to the main lecture sequence. Read Chapter 2 — Least Squares Estimation, Chapter 5 — Multiple Regression and Categorical Predictors, and Chapter 8 — Diagnostics and Model Adequacy first: least squares, coefficient interpretation, and diagnostics are the useful prerequisites. The R examples use the quantreg package; install it once with install.packages("quantreg") if needed.

When the mean is not the whole question

OLS describes a conditional mean. In some studies, the scientific question concerns a typical response, a low outcome, or an unusually high outcome. For example, average delivery time and the 90th percentile of delivery time answer different service-planning questions. A covariate may change the spread of the response even when its relationship with the mean is simple.

Definition 1 (Conditional Quantile) For \(0<\tau<1\), define

\[ Q_\tau(Y\mid \mathbf{X}=\mathbf{x})=\inf\{q:P(Y\le q\mid \mathbf{X}=\mathbf{x})\ge\tau\}. \]

Here \(Y\) is a scalar response, \(\mathbf{X}\) is a random predictor vector, and \(\mathbf{x}\) is its observed value. A linear quantile model specifies \(Q_\tau(Y\mid \mathbf{X}=\mathbf{x})=\mathbf{x}^\top\boldsymbol{\beta}_\tau\). The coefficients may vary with \(\tau\). At \(\tau=0.5\), the target is a conditional median.

Interpretation. Without an interaction involving \(x_j\), \(\beta_{\tau,j}\) is the change in the modeled conditional \(\tau\)-quantile for a one-unit change in \(x_j\), holding the other predictors fixed. This compares conditional distributions. It does not track a particular person who remains at the same percentile, and it is not automatically a causal effect.

ImportantUse all observations to estimate a conditional quantile

Fitting the 90th percentile does not mean selecting the largest 10% of responses and fitting OLS to them. Quantile regression uses the full dataset with a loss function that gives different penalties to positive and negative residuals. A marginal percentile of all responses is also different from a percentile conditional on the predictors.

Why a different loss? Least squares penalizes a residual \(u=y-\mathbf{x}^\top\mathbf{b}\) by \(u^2\). Quantile regression instead uses the check loss, also called pinball loss:

\[ \begin{aligned} \rho_\tau(u)&=u\{\tau-\mathbb{1}\{u<0\}\} =\begin{cases}\tau u,&u\ge0,\\(1-\tau)(-u),&u<0,\end{cases}\\ \hat{\boldsymbol{\beta}}_\tau&\in\arg\min_{\mathbf{b}}\sum_{i=1}^n\rho_\tau(y_i-\mathbf{x}_i^\top\mathbf{b}). \end{aligned} \]

At \(\tau=0.5\), minimizing check loss is equivalent to minimizing the sum of absolute residuals. At \(\tau=0.9\), underpredicting the response costs nine times as much per unit as overpredicting it. The fitted target consequently moves toward the upper tail.

Proposition 1 (Why Check Loss Targets a Quantile) If \(\mathbb{E}(|Y|)<\infty\), the minimizers of \(\mathbb{E}\{\rho_\tau(Y-a)\}\) are characterized by \(F(a-)\le\tau\le F(a)\); the conditional-quantile definition above selects the lower endpoint when the minimizer is not unique. Where \(F\) is continuous, the derivative with respect to \(a\) is \(F(a)-\tau\). It is nonpositive below a minimizing quantile and nonnegative above it. The same argument applies to the conditional distribution at a fixed \(\mathbf{x}\).

Scope of the linear model. The population argument allows a different quantile at every \(\mathbf{x}\). A linear quantile regression restricts these values to one linear function. If the true conditional quantile is nonlinear, the fitted line is a best approximation under the chosen loss, rather than an exact description at every predictor value.

Show R code
library(ggplot2)
theme_set(theme_minimal(base_size = 11) +
            theme(legend.position = "bottom",
              panel.grid.minor = element_blank()))

check_loss <- function(u, tau) u * (tau - (u < 0))
taus <- c(0.1, 0.5, 0.9)
loss <- expand.grid(u = seq(-3, 3, length.out = 301), tau = taus)
loss$value <- check_loss(loss$u, loss$tau)
ggplot(loss, aes(u, value, colour = factor(tau))) +
  geom_vline(xintercept = 0, linetype = "dotted", colour = "grey70") +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = c("steelblue4", "grey35", "darkorange3")) +
  labs(x = "Residual: observed - predicted", y = "Check loss",
    colour = "Quantile")
Figure 1: Check loss for three quantiles. A positive residual means the prediction is too low. The upper-quantile loss penalizes that side more strongly.

A mean model can miss changing spread

Consider the simulated model \(Y=5+2x+(0.5+0.8x)Z\), with \(0\le x\le4\) and \(Z\sim N(0,1)\) independent of \(x\). Its conditional mean and median are both \(5+2x\), while

\[ Q_\tau(Y\mid x)=5+2x+(0.5+0.8x)\Phi^{-1}(\tau). \]

The conditional-quantile slopes are \(2+0.8\Phi^{-1}(\tau)\). Different slopes here reflect changing spread, even though every observation follows the same data-generating model. Quantile regression need not be reserved for skewed distributions.

Show R code
set.seed(8561)
n <- 350
x <- runif(n, 0, 4)
dat_qr <- data.frame(x, y = 5 + 2 * x +
                      (0.5 + 0.8 * x) * rnorm(n))
fit_mean <- lm(y ~ x, data = dat_qr)
fit_quantiles <- quantreg::rq(y ~ x, tau = taus, data = dat_qr)
round(coef(fit_quantiles), 3)
            tau= 0.1 tau= 0.5 tau= 0.9
(Intercept)    4.509    5.078    5.850
x              0.768    1.816    2.965
Show R code
data.frame(tau = taus, true_slope = 2 + 0.8 * qnorm(taus),
           estimated_slope = unname(coef(fit_quantiles)["x", ]))
  tau true_slope estimated_slope
1 0.1  0.9747587       0.7676385
2 0.5  2.0000000       1.8156789
3 0.9  3.0252413       2.9651924

The matrix returned by coef() contains one column per requested quantile. R’s rq() documentation describes its formula interface and the tau argument.

Show R code
new_x <- data.frame(x = seq(0, 4, length.out = 200))
qhat <- predict(fit_quantiles, newdata = new_x)
quantile_lines <- data.frame(x = rep(new_x$x, times = length(taus)),
                              prediction = as.vector(qhat),
                              target = rep(paste("Quantile", taus),
                                each = nrow(new_x)))
prediction_lines <- rbind(quantile_lines,
  data.frame(x = new_x$x, prediction = predict(fit_mean, newdata = new_x),
              target = "OLS mean"))
prediction_lines$target <- factor(prediction_lines$target,
                                  levels = c(paste("Quantile", taus),
                                    "OLS mean"))
ggplot(dat_qr, aes(x, y)) +
  geom_point(alpha = 0.35, colour = "grey40", size = 1.5) +
  geom_line(data = prediction_lines,
             aes(y = prediction, colour = target, linetype = target),
               linewidth = 0.8) +
  scale_colour_manual(values = c("steelblue4", "grey35", "darkorange3",
    "black")) +
  scale_linetype_manual(values = c("solid", "solid", "solid", "dashed")) +
  labs(x = "x", y = "Response", colour = NULL, linetype = NULL)
Figure 2: Fitted conditional quantiles and the OLS mean in simulated data with increasing spread. The growing distance between the lower and upper quantiles reveals variation that a mean line alone does not describe.

Read the graph. Compare the predicted 10th, 50th, and 90th percentiles at a fixed \(x\). The gap between the 10th and 90th percentiles describes an interval containing the central 80% of a continuous conditional distribution when the quantile model is correct. The fitted gap is not a confidence interval for the mean or for a regression coefficient, and its predictive coverage is not guaranteed in a finite sample.

Uncertainty and practical checks

Quantile regression avoids a normal-error or constant-variance requirement for defining its target, but it is not assumption-free. Reliable estimation and inference still depend on the quantile specification, the design, the sampling structure, and sufficient information near the quantile of interest. Very extreme quantiles are difficult to estimate from small datasets.

A pairs-bootstrap illustration. These simulated observations are independent draws, so resampling whole \((x_i,y_i)\) pairs is appropriate. The code requests a bootstrap standard error for the median-regression slope and forms an approximate normal interval. The replicate count is kept moderate for the notes; increase it and check Monte Carlo stability in a substantive analysis.

Show R code
set.seed(8561)
fit_median <- quantreg::rq(y ~ x, tau = 0.5, data = dat_qr)
median_summary <- summary(fit_median, se = "boot", R = 499,
                          bsmethod = "xy")
b <- coef(fit_median)["x"]
se_b <- median_summary$coefficients["x", "Std. Error"]
c(slope = unname(b), bootstrap_SE = se_b,
  lower = unname(b - qnorm(0.975) * se_b),
  upper = unname(b + qnorm(0.975) * se_b))
       slope bootstrap_SE        lower        upper
    1.815679     0.136845     1.547468     2.083890 

This is not the exact finite-sample OLS t interval from Chapter 3 — Distribution Theory of OLS and Inference. For dependent or clustered observations, an independent-pairs bootstrap generally fails to preserve the sampling structure. The package’s summary.rq() documentation describes alternative inference methods.

ImportantTwo fitted quantiles need a joint comparison

If one slope is significant at \(\tau=0.1\) and another is not significant at \(\tau=0.9\), that does not establish that the slopes differ. They are estimated from the same observations and are dependent. A difference such as \(\beta_{0.9,j}-\beta_{0.1,j}\) requires joint inference that retains this dependence, for example refitting both quantiles in each common bootstrap sample.

Before interpreting a fit: check the functional form, influential predictor values, the amount of information in the tails, and whether separately fitted quantiles cross. Linear loss is less sensitive to very large response residuals than squared loss, but it does not eliminate leverage problems. Crossing fitted quantiles cannot describe a valid ordered distribution at that \(x\); investigate misspecification or use methods that enforce ordering rather than silently relabeling the curves.

Evaluate the stated target. For held-out data, compare candidate predictions at the same \(\tau\) using mean check loss. A coverage diagnostic compares the fraction of responses below the predicted quantile with \(\tau\), including checks across meaningful predictor ranges. Good overall coverage can hide poor conditional coverage. Use splits that match future prediction, as in Chapter 10 — Model Assessment and Prediction.

Application and self-study

NoteApplication: household expenditure

The engel dataset bundled with quantreg contains income and food expenditure. Use data("engel", package = "quantreg"), plot the variables, and fit the 0.1, 0.5, and 0.9 conditional quantiles of foodexp given income. Explain how a household-expenditure quantile differs from mean expenditure. Treat the fitted relationships as descriptive associations, and avoid extrapolating beyond the observed incomes.

  1. Use the loss definition to explain why the 90th percentile is penalized more for underprediction.
  2. Derive the true intercept and slope at each quantile in the simulated example. Explain why the median and mean agree here.
  3. Predict the 10th and 90th percentiles at \(x=1\) and \(x=3\). Explain what their widening gap means.
  4. Design a bootstrap for comparing two quantile slopes. State which observations must stay together if the data are clustered.

Suggested reading. Start with the quantreg package documentation, then the vignette available through vignette("rq", package = "quantreg"). Roger Koenker’s Quantile Regression develops the estimation and inference theory in greater depth.