What type of data should I generate to observe/amplify a crossing problem in quantile regression?

1. Background

Crossing problem in quantile regression can be observed when we want to estimate several conditional quantiles (e.g. τ = 0.1, 0.2, . . . , 0.9), as two or more estimated conditional quantile functions could cross or overlap. However, a higher percentile must have more points below it than a lower percentile. That is, we should not have more points below the 75th percentile than we do at the 80th percentile. This would violate the basic definition of percentiles, making the quantile regression invalid. (Takeuchi et al., 2005 ; Saleh, 2021)

2. Question

This being said, I would like to generate synthetic data to observe and maximise this issue. To perform the quantile regression, I am using the quantreg framework in R. (https://cran.r-project.org/web/packages/quantreg/quantreg.pdf)

Since starting my research, I have tried to generate various datasets that present mean or SD non-linearly changes and have yet to be successful.

Does anyone know how to do this, or any direction to look for?

1 Answer

If we are to avoid trivialities related to outlying values, evidently crossing is an indication the regression model isn't correct. So, all you need is suitably poor goodness of fit along with a little leverage and not too many data points.

For example, generate conditional distributions, some of whose quantiles vary nonlinearly with the explanatory variables, and fit linear models to the quantiles. Something like $$Y\mid x \sim\mathcal{N}(0,\exp(x))$$ ought to do the trick.

Here is an example, zoomed at the right to show the details of the crossing regression lines. It was created by setting df = 1 at the beginning of the R code below, which generates 601 data points.

The colors of the fitted lines start dark at low quantiles and graduate to lighter colors at higher quantiles. They are in the correct order for $$x gt 2$$ or so, but cross at lower values of $$x.$$

The crossings are severe: ultimately, near $$x=0,$$ the fitted $$5^\text{th}$$ percentile exceeds the fitted $$95^\text{th}$$ percentile!

Even more complex regression models are subject to this problem. Here are the same data fitted with a natural spline with three degrees of freedom (using the code below unchanged).

The nonlinearity is achieved through the strong variation of the conditional response variance, as suggested above. The points at the right are given greater leverage by sampling them disproportionately more often. This causes the trends at the right to overwhelm the fewer amounts of data at the left.

These patterns are consistent with other datasets (comment out the set.seed call to see). In most cases the software does not issue a warning that the fit might not be unique.

library(quantreg)
#
# Generate a dataset.
#
set.seed(17)
n <- 601                  # Size of dataset
lambda <- 5               # Upper limit of x
N <- 31                   # Number of distinct possible values of x, equally spaced
probs <- seq_len(N) ^ 2.5 # Relative sampling probabilities for x
df <- 3                   # Spline degrees of freedom (1 = linear)
X <- data.frame(x = sort(sample(seq(0, lambda, length.out = N),
n, replace = TRUE, prob = probs)))
X\$y <- with(X, rnorm(n, sd = exp(x)))
#
# Fit a range of quantiles, not too extreme.
#
tau <- c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95) # quantiles
fit <- rq(y ~ splines:::ns(x, df), tau = tau, X)
#
# Compute the predicted values.
#
x.0 <-seq(0, lambda, length.out = 501)
y.hat <- predict(fit, newdata = data.frame(x = x.0))
#
# Plot the data and the results.
#
hue <- function(x) log10(1/2 - abs(x - 1/2))
hmax <- max(abs(hue(tau)))
par(mfrow = c(1, 2))
for (mu in c(1, 1/10))
with(X, {
plot(range(x), range(c(y, y.hat)) * mu, type = "n",
main = bquote(Quantiles * phantom(.) * .(paste(tau, collapse = ", "))),
cex.main = 1,
xlab = "x", ylab = "y")
invisible(sapply(seq_len(ncol(y.hat)), function(i)
lines(x.0, y.hat[,i],
col = hsv(-0.8 / hmax * hue(tau[i]), .9, (1 + tau[i])/2), lwd = 2)))
points(x, y, col = "#00000060")
})
par(mfrow = c(1, 1))

• Thank you very much for your detailed answer! I will need some time to study it, and I will write back my feedback. Commented Jan 22, 2023 at 18:37
• To understand your answer, I have tried to manipulate the graphic representation of the crossings. Thus, I tried these lines: rq_model = summary(rq(y ~ splines:::ns(x, df), data=X, tau=1:99/100), se = "boot", bsmethod="xy", R=100, level=0.95) plot(rq_model, mfrow = c(2,3)) If I am correct, this should give me the intercept and the regression curves of your regression. I expected it to be awkward, but it surprisingly looks pretty normal despite all the crossings. How would you interpret the results of these 2 lines? Is it even interpretable? if not, what am I missing to say so? Commented Jan 28, 2023 at 15:55