Why are the predictions of a quantile regression model changed by an increasing transformation of the DV?

Question

Take a dataset and suppose we fit two quantile regression models to it, one with the raw dependent variable (DV) and one with the logged DV. Then look at each model's predictions for the training data, reversing the transformation in the case of the logged model. Here's an example in R:

library(quantreg)

set.seed(1)

x1 = rnorm(100)
y = exp(x1 + 3*rnorm(100))

m = rq(y ~ x1, tau = .5)
p1 = predict(m)

m = rq(log(y) ~ x1, tau = .5)
p2 = exp(predict(m))

print(head(p1))
print(head(p2))

The two models give different predictions:

        1         2         3         4         5         6 
0.4631787 0.8619277 0.3602179 1.5567692 0.9337257 0.3676802 
        1         2         3         4         5         6 
0.2776789 0.7859069 0.2122640 4.8162168 0.9478212 0.2164372

But how can this be? The models are fit in terms of the quantiles of the DV, and the natural logarithm is an increasing function, so it preserves quantile ranks.

Answer 1

The equivariance to monotone transformations property that implies $Q_q(y \vert x)=\exp \{Q_q(\ln y|x )\}$ is exact only if the conditional quantile function is correctly specified. This is unlikely to be the case in practice, and is not the case in your simulation, since $\exp \{x+\varepsilon\} \ne x + \varepsilon$. The only case where the linear model will be exact is when all regressors are discrete and we specify a fully saturated model with dummy variables as regressors that exhaust all the possible interactions.

This will give you much better results, since I am adding a constant to avoid undefined logs of zeros and negatives, instead of exponentiating:

library(quantreg)
set.seed(1)
x1 = rnorm(100)
y = 10 + x1 + 3*rnorm(100)
m = rq(y ~1 + x1, tau = .5)
p1 = predict(m)
m2 = rq(log(y) ~1 + x1, tau = .5)
p2 = exp(predict(m2))
print(head(p1))
print(head(p2))

Answer 2

This isn't a complete answer, but it may help.

The covariate in your model may be a red herring, the issue is reproducible with an intercept only model

set.seed(1)

y = runif(100, 0, 1)  # Need to keep it positive for the log.

m = rq(y ~ 1, tau = .5)
p1 = predict(m)

m2 = rq(log(y) ~ 1, tau = .5)
p2 = exp(predict(m2))

print(head(p1))
print(head(p2))

which results in

> print(head(p1))
        1         2         3         4         5         6 
0.4820801 0.4820801 0.4820801 0.4820801 0.4820801 0.4820801 
> print(head(p2))
        1         2         3         4         5         6 
0.4935413 0.4935413 0.4935413 0.4935413 0.4935413 0.4935413 

But if we redo this whole thing with an odd number of data points

set.seed(1)

y = runif(101, 0, 1)  # <- Now it's odd!

m = rq(y ~ 1, tau = .5)
p1 = predict(m)

m2 = rq(log(y) ~ 1, tau = .5)
p2 = exp(predict(m2))

print(head(p1))
print(head(p2))

we get

    1         2         3         4         5         6 
0.4935413 0.4935413 0.4935413 0.4935413 0.4935413 0.4935413 
> print(head(p2))
        1         2         3         4         5         6 
0.4935413 0.4935413 0.4935413 0.4935413 0.4935413 0.4935413

Which is on the nose.

It's easy to reproduce this without the quantile regression call

> x <- c(1, 2, 3, 4, 5, 6)
> median(x)
[1] 3.5
> exp(median(log(x)))
[1] 3.464102
> 
> x <- c(1, 2, 3, 4, 5, 6, 7)
> median(x)
[1] 4
> exp(median(log(x)))
[1] 4

I suspect this has something to do with the behaviour you're seeing, but it's not clear how to incorporate the covariate into the argument.

Answer 3

The fitted values are calculated by substituting in the vector values of the independent variable to the coefficient estimates (in this case of the median).

$y = mx +c$.

There is a direct linear relationship between the original (m) and log (m2) predictions (ie. same x with different slope and intercept).

Here's a code snippet to get to the nuts and bolts of it:

set.seed(1)
x1 = rnorm(100)
y = exp(x1 + 3*rnorm(100))
X <- as.matrix(cbind("(Intercept)"=rep(1,length(x1)),x1))
rq1 <- rq.fit.fnb(X, y, tau = 0.5)
resid <- (y - X %*% rq1$coefficients)
fit <- X %*% rq1$coefficients

rq2 <- rq.fit.fnb(X, log(y), tau = 0.5)
residlog <- (y - X %*% rq2$coefficients)
fitlog <- X %*% rq2$coefficients

model <- coef(lm(fit~fitlog))
plot(fitlog ~ fit)
head(fit)
head((fitlog*model[2])+model[1])
      [,1]
[1,] 0.4769196
[2,] 0.8292802
[3,] 0.3859367
[4,] 1.4432873
[5,] 0.8927255
[6,] 0.3925308

      [,1]
[1,] 0.4769196
[2,] 0.8292802
[3,] 0.3859367
[4,] 1.4432873
[5,] 0.8927255
[6,] 0.3925308

Refs for the mathematics of quantile regression:

https://projecteuclid.org/download/pdf_1/euclid.ss/1030037960

http://www.econ.uiuc.edu/~roger/research/rq/rq.pdf