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))
