Quantile regression on a constant: is this different from unconditional quantile?

Question

With a linear model, estimating an OLS regression of y on a constant only will give us the mean of y. I was wondering whether the same would not hold with a quantile regression?

In other words, is running a quantile regression of y on a constant equivalent to estimating the unconditional quantile? My intuition was telling me that yes, but a quick example in R shows me this is not necessarily the case? Or are they supposed to be equivalent, yet there is some small machine imprecision between algorithms?

Maybe this is linked to the post on What is the difference between conditional and unconditional quantile regression? My motivation for wondering about the quantile regression comes from the quantile pseudo $R^2$, which contains precisely a quantile regression on a constant, see R-squared in quantile regression

> set.seed(1)
> y <- rnorm(100)
> 
> ## Numerically equivalent for linear regression
> all.equal(coef(lm(y~1)), mean(y), check.attributes = FALSE)
> #> [1] TRUE
> sapply(1:9, function(i) quantile(y, probs=0.3, type=i))
       30%        30%        30%        30%        30%        30%        30%        30%        30% 
-0.3942900 -0.3807557 -0.3942900 -0.3942900 -0.3807557 -0.3861694 -0.3753420 -0.3825603 -0.3821091 
> coef(quantreg::rq(y~1, method="fn", tau=0.3))
(Intercept) 
 -0.3672215 
> coef(quantreg::rq(y~1, method="br", tau=0.3))
(Intercept) 
 -0.3672215 
Warning: In rq.fit.br(x, y, tau = tau, ...) : Solution may be nonunique

^{Created on 2021-06-08 by the reprex package (v2.0.0)}

Answer 1

You can generally reproduce the desired behaviour

1. if you use the inverse ECDF method to calculate the empirical quantiles and
2. if $np$ is not an integer ($n$ is the sample size, $p$ the probability).

Here the code:

library(quantreg)

#====================================
# n*probs is an integer
#====================================

n <- 100
set.seed(1)
y <- rnorm(n)
p <- 0.3

if (n*p == trunc(n*p)) {
  print("quantile not uniquely defined")
}

# Comparison: FALSE
all.equal(
  coef(rq(y ~ 1, tau = p)),
  quantile(y, probs = p, type = 1),
  check.attributes = FALSE
)

#====================================
# n*probs not an integer
#====================================

n <- 99
set.seed(1)
y <- rnorm(n)
p <- 0.3

if (n*p == trunc(n*p)) {
  print("quantile not uniquely defined")
}

# Comparison: TRUE
all.equal(
  coef(rq(y ~ 1, tau = p)),
  quantile(y, probs = p, type = 1),
  check.attributes = FALSE
)

#====================================
# Factor as covariable
#====================================

set.seed(1)
x <- rep(0:1, each = 99)
y <- rnorm(99 * 2) + x

# Comparison: TRUE
all.equal(
  coef(rq(y ~ 0 + factor(x), tau = 0.3)),
  c(by(y, x, quantile, type = 1, probs = 0.3)),
  check.attributes = FALSE
)

The last snippet shows how to calculate quantiles within subgroups with the same trick.