# Stepwise quantile regression: What's the reason behind these strange results?

So I am attempting to build a model using quantile regression & am using stepwise regression for initial data exploration. I'm well aware that stepwise methods are widely frowned upon & am using this approach purely for exploratory purposes, as I have a very large amount of variables, most of which could theoretically be relevant to our phenomenon.

I am seeing some strange behavior that makes me wonder if something in my data is inappropriate for this type of analysis & possibly points to gaps in my conceptual understanding of quantile regression and/or stepwise selection. I'm working in R, though I suspect it's a conceptual issue & not a problem with R.

My dataset has 59 potential predictor variables & n is close to 20,000.

The code is essentially just:

stepfit <- step(rq(response ~ ., data = data, tau = 0.99), direction = "both")


Here are the main red flags I'm seeing (though they raise somewhat separate questions, I'm presenting them all, since I wonder if the issues are in fact related):

• Running stepwise at tau = 0.9 produces a final model with 7 variables and AIC in the neighborhood of 16,000. The model looks pretty reasonable from the perspective of the phenomenon we're studying. However, running the same at tau = 0.99 produces a monster model that includes almost 90% of our variables, a lot of them with bizarre giant coefficients, & AIC at something like 45,000. Why would a change in tau lead to such a dramatic change in results?
• In a previous version of the data, which had about 70 variables, tau = 0.9 would run (yielding the SAME parsimonious model as in the bullet point above), whereas tau = 0.99 would fail due to singularity in the matrix. We "resolved" this problem by paring down our potential predictor variables -- though we did this on the basis of theory, not because anything about them was obviously liable to lead to matrix singularity (they were not involved in collinearity or multiples of one another or anything like that). Is there a general reason why a higher tau and/or higher number of predictor variables would lead to increased likelihood of singularity?
• For all cases, R will report some warning message like "973 non-positive fis". The number is higher for higher tau & more variables.

Some features of our data that I suspect are causing these (& possibly other less visible?) issues:

• The response variable is very heavily zero-dominated (like at least 90%). However, we can't just toss this out, because these zeroes are meaningful and to be expected.
• The response variable is not really continuous -- it can only take values of 0 (the dominant value), 0.33333, 1.66667, 2, 2.3333, 2.6667, and 3. It is a composite index of about 4 other variables that ARE continuous, but were binned (into "0", "1", "2" etc) & averaged to get the index. With 20,000 data points, there are many, many identical response values.
• 5 of our predictor variables are binary & mutually exclusive (if one has "1", the others will have "0").
• Because we suspected that one or more of the above 3 bullet points would lead to matrix singularity issues, we added a tiny amount of random noise to everything (all predictors & the response), as recommended elsewhere. Doing so + paring down the number of predictors from 70+ to 59 seemed to "resolve" singularity issues, but at that point, the output began to follow the pattern described in the very first bullet point above.

I suspect one or several of the characteristics of our data is leading to the strange outputs (radically different at tau = 0.9 vs. 0.99). Is there a conceptual reason for this & how might we be able to address it?? Do any of the data characteristics described clearly invalidate the approach here?

• Quantile regression at the tails is difficult, often you don't have a lot of data and the algorithm has a hard time . This could explain the strange results. You might want to consider extreme value quantile regression, but I am not familiar enough with this method to recommend anything. – Repmat Sep 27 '17 at 17:29

By setting $\tau=0.99$ you are essentially fitting to the noise at the extreme high-end of your sampling distribution. If you imagine the only-slightly-more-extreme case with $\tau=1.0$, you would be fitting exclusively to the points comprising the upper-envelope of your data, in other words sparsely-distributed points which are almost guaranteed to be unrepresentative outliers. See the plot below for reference. 