Quantile Regression in R: quantreg Package, Model Beyond the Mean

Quantile regression, fit in R with the quantreg package's rq() function, models the conditional median (or any other quantile) of an outcome instead of its mean. It reveals how predictors act on low, typical, and high responses, captures heteroscedasticity, and resists outliers.

When should you reach for quantile regression instead of lm()?

Decisions about rent caps, hospital stay limits, or insurance premiums care about the tails, not just the average. Ordinary least squares collapses every bit of variability onto a single line through the mean, hiding the fact that one predictor may barely shift the median yet stretch the upper quartile sharply upward. Quantile regression, introduced by Koenker and Bassett (1978), traces that behaviour quantile by quantile.

The classic example ships inside quantreg itself: Engel's 1857 study of how food expenditure depends on household income. Fit the median regression with rq() and the mean regression with lm(), then compare.

The median slope (0.56) is steeper than the mean slope (0.49), and the median intercept is lower. That gap is not noise. Richer households spread their food spending further above the line than poorer households pull it below, so the mean gets tugged toward the high earners. The median tracks the typical household, which is what many policy decisions actually need.

How does rq() fit a quantile instead of a mean?

OLS finds the line that minimises the sum of squared residuals. That loss treats positive and negative errors symmetrically, so the solution lands on the conditional mean. Quantile regression swaps squared loss for an asymmetric absolute loss called the check function. By tilting the penalty, it pushes the fit toward any quantile you name.

Figure 1: The quantreg workflow, fit with rq(), inspect coefficients, bootstrap standard errors, then plot the process.

The check function for quantile $\tau$ is:

$$\rho_\tau(u) = u \cdot (\tau - \mathbb{1}(u < 0))$$

Where:

• $u$ = the residual ($y - \hat{y}$)
• $\tau$ = the target quantile, between 0 and 1
• $\mathbb{1}(u < 0)$ = 1 when the residual is negative, 0 otherwise

At $\tau = 0.5$, both positive and negative residuals get weight 0.5, which recovers least absolute deviations (median regression). At $\tau = 0.9$, positive residuals get weight 0.9 and negative residuals get weight 0.1, so the fit pushes up until only 10% of observations sit above it. The optimiser solves a linear programme under the hood, which is why rq() handles any tau without Newton iterations.

Feeding a vector of taus to rq() fits all of them at once:

The income slope rises from 0.40 at the 10th percentile to 0.69 at the 90th, a 70% increase. That tells you high-spending households scale their food bill with income much more aggressively than frugal households do. OLS returns a single number (0.49) and forfeits this story.

How do you read the coefficients at different quantiles?

A coefficient at tau = 0.9 says: "if x increases by 1, the 90th percentile of y increases by this much." It does not say anything about the average household, only the one sitting at that specific quantile of the conditional distribution.

When slopes differ across taus, the effect of x is heterogeneous: it acts on low, middle, and high responders differently. When slopes are close to constant across taus, the effect is homogeneous and OLS is telling the whole story.

Scanning the income row, the slope climbs monotonically from 0.40 to 0.69. That monotonic rise is the signature of location-scale heterogeneity: both the mean and the spread of foodexp grow with income. If you cared about policy for the poorest 10%, 0.40 is your elasticity, not the OLS 0.49.

How do you get reliable standard errors and confidence intervals?

summary(fit) returns confidence intervals by default, but the method matters because heteroscedasticity breaks simple formulas. The se argument of summary.rq() lets you choose:

When in doubt, bootstrap. It makes no distributional assumptions and handles heteroscedasticity naturally.

The bootstrap standard error for income (0.033) gives a 95% confidence interval of roughly $[0.495, 0.625]$, so the median elasticity is estimated to about two decimal places. The rank-based default would sometimes return narrower intervals, but those can undercover under heteroscedasticity, so the bootstrap is the safer default on real data.

How do you visualize effects across the entire distribution?

Reading coefficients at three taus is fine, but the richer picture is a quantile process plot, one panel per predictor showing the coefficient as a continuous function of tau. The dashed horizontal line marks the OLS estimate, and the shaded band is the confidence interval.

Fit rq() at a fine grid of taus, wrap the fit in summary() (bootstrap SE recommended), then pipe to plot().

Two panels appear, one for the intercept and one for the income slope. The income slope curve rises steadily from about 0.40 to 0.70 as tau moves from 0.05 to 0.95, well outside the OLS line for most of its range. That visible upward slope is the signature of conditional scale heterogeneity. A flat band hugging the OLS line would say "OLS is enough"; a tilted band, like this one, says "quantile regression is worth the extra work."

How does quantile regression expose heteroscedasticity?

Classical regression detects non-constant variance with a residual-vs-fitted plot: a visible funnel means the spread grows with x. Quantile regression offers a quantitative alternative. If the slope at tau=0.9 exceeds the slope at tau=0.1, the upper tail grows faster than the lower tail, a direct fingerprint of heteroscedasticity.

The 10th-percentile slope is 1.02 and the 90th-percentile slope is 1.98, almost double. The true generating process made residual SD proportional to $0.5 + 0.4x$, so the tails pull apart as x grows. Quantile regression reads this off directly: coef(fit)[2, "tau= 0.9"] - coef(fit)[2, "tau= 0.1"] estimates how much faster the upper tail climbs.

Practice Exercises

Exercise 1: Weight slope across mtcars quantiles

Fit rq() on mtcars for mpg ~ wt + cyl at tau = c(0.25, 0.5, 0.75). Extract the wt slope at each tau into a named numeric vector called my_wt_slopes.

Exercise 2: Find the heterogeneous predictor on airquality

Fit rq() on na.omit(airquality) for Ozone ~ Solar.R + Wind + Temp at tau = c(0.1, 0.5, 0.9) with bootstrap SE (R = 300). Identify which predictor's slope changes the most between the 10th and 90th percentiles, assign the predictor name (string) to my_hetero.

Exercise 3: Engel income elasticity at each decile

Fit rq() on Engel for foodexp ~ income at tau = seq(0.1, 0.9, by = 0.1). Return a data frame my_engel_decile with columns tau and income_slope.

Complete Example

A full airquality workflow, from data cleaning to interpretation:

Reading the four panels (intercept plus three predictors): Temp climbs from 0.75 at tau=0.1 to 1.92 at tau=0.9, meaning each extra degree lifts the upper-tail ozone more than twice as much as the lower-tail ozone. Wind tilts steeply negative in the upper tail, so windy days clip the highest readings sharply. Solar.R shows modest heterogeneity, mostly a location shift. A single OLS model would have reported one average effect per predictor and hidden all three patterns.

Summary

Key takeaways from fitting quantile regression with quantreg:

Figure 2: Quantile regression essentials, fit, inference, visualization, and core use cases.

Quantile regression earns its keep whenever the tails of the response matter, whenever variance changes with predictors, or whenever outliers corrupt OLS. For problems where the conditional mean is all you need and residual variance is constant, lm() remains simpler and sufficient. The quantreg package also extends beyond rq(): rqss() adds nonparametric smoothing, crq() handles censored responses, and nlrq() fits nonlinear quantile models.

References

1. Koenker, R. & Bassett, G. (1978). Regression Quantiles. Econometrica, 46(1), 33-50. Link
2. Koenker, R. & Hallock, K. (2001). Quantile Regression. Journal of Economic Perspectives, 15(4), 143-156. Link
3. Koenker, R. (2005). Quantile Regression. Cambridge University Press.
4. Koenker, R. quantreg package CRAN vignette. Link
5. Cade, B. S. & Noon, B. R. (2003). A gentle introduction to quantile regression for ecologists. Frontiers in Ecology and the Environment, 1(8), 412-420. Link001%5B0412:AGITQR%5D2.0.CO%3B2)
6. CRAN, Package 'quantreg' reference manual. Link

Continue Learning

• Polynomial and Spline Regression in R, the parent tutorial on modelling curvature without leaving the conditional mean.
• Linear Regression in R, the OLS baseline quantile regression generalises.
• Beta Regression in R, for bounded responses where neither OLS nor quantile regression is appropriate.