# Theory understanding behind quantile regression

As part of my studies in ecology, I am trying to reproduce a quantile regression method described in Karlsson et. al (2022) (source: https://arxiv.org/pdf/2202.02206.pdf)

1. My starting point is this:

Throughout this paper, the response variable will be julian day. For single species models, yj will be the julian arrival day of individual j ∈ {1, . . . , n}. With each observation comes a covariate vector xj = (1, tj, xj1, . . . , xjp), consisting of an intercept term 1, the year tj that individual j was recorded, and p additional covariates.

So, I assume that their dataset looks like this (+ some covariates columns that I will not take into account):

individual year julian_day
1 y(1) j(1)
... ... ...
n y(n) j(n)

Basically, 1 row = 1 bird, the year it was recorded and the day it was recorded (julian day format).

I have had to create my own dataset, and I think that I have managed to do a similar one, using citizen science data (ebird) on a random species. However they are some differences:

• Firstly, and it might be of importance regarding my question: instead of "individuals" I used "observation" which means that I do not take into account the abundance.

For example, 1 observation could be translated by: 1) I have seen 10 birds at the same time on julian day 35, year 1990 = 1 observation; 2) I have seen 1 bird on julian day 70, year 2010 = 1 observation; 3) I have seen 340 birds on julian day 60, year 1984 = 1 observation.

• Secondly, I divided the time period in decades classes and then attributed a continuous time frame (e.g. 1950-1959 = 1; 1960-1969 = 2 and so on).

Which gives me a dataset like this one:

1 y(1) j(1) d(1) t(1)
... ... ... ... ...
n y(n) j(n) d(n) t(n)
1. Then, they set up quantile regression models, using R, 'quantreg' package:

The quantile is specified by τ , ε denotes the error term and there are n = 2 203 observations. Models were fitted for τ = 0.01, 0.02, . . . , 0.99

And importantly:

we use the Frisch-Newton interior point method (Portnoy and Koenker, 1997) provided in the R package quantreg (Koenker, 2021).

So did I:

model = rq(julian_day~time, tau = 2:98/100, method = "fn", data= .)

1. Parameters estimates:

we checked if the intercept estimates were monotonically increasing in τ .

from Karlsson et al. 2022

NB: the model that I am trying to implement is corresponding to the red one on their representation.

I have built like this 12 models (grouped by a latitude condition) and checked if the intercepts estimates were monotonically increasing in τ, but here comes my difficulty. I have nothing that looks like monotonous. Let me give you two examples (but I can give more if needed):

As their paper is about comparing different mathematical approaches, they do not interpret their results on a biological point of view and what is interesting for me stops here.

I feel that it is not right if my intercet is not monotonous, so my questions are: what does it mean if my intercept is not monotonous? Can I still interpet my quantile regression models? Have I done something wrong from a statistical point of view when I built my dataset, that I am unaware of and that leads to these results? If I want to give sense to my analysis, towards which direction should I aim?

Any advice or thought on my methodology is welcome, I tried my best to explain with as much details as I could a problem that is not easy for me, I really hope that it is clear, but I might have forgotten to give important details: if so, please excuse me and I can still update my post :)

Thank you very much for the help.

UPDATE:

I have kept working and exploring, therefore I can add some more insights.

1. My first attempt was to try to rebuild my dataset as closer as I could to the paper's dataset, so I built some models using this time julian_day~year calculated over abundance (instead of observations). The result is not different, I still do not obtain a monotonous intercept.

2. As my previous attempt just increased a lot the size of my dataset, I went back to my julian_day~time models based on observations to reduce the calculation time for trying different things.

3. That said, I noticed this warning message (I do not know why I did not spot it before...) while computing the models:

 Warning message:
In rq.fit.fnb(x, y, tau = tau, ...) :
Error info =  2 in stepy: possibly singular design


This message lead me to this topic: Cause of singularity in matrix for quantile regression where @Justina Pinch wrote:

I got around it by 'jittering' the data, adding a very small amount of random noise to the response values using rnorm()

I am currently trying this idea, but one question that I have here is: how is it possible that I get this issue, trying to keep closer to the paper method/dataset structure as much as I can, while the authors do not talk at all about having this issue? Is it possible that they had to get through this without mentionning it in the article? Because it should lead to adress some questions like: what amount of jittering do I had and why? Does it have any impact on the coefficient estimation? And probably some more...

1. Finally, I did a plot of my data with the following quantile lines: 0, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 1. Here is the result, where we can see as described by @Justina Pinch:

Duplicating observations (for a single x value, multiple responses) increases chances of singularity.

Still, this feels strange to me as I do not see how it is possible to have any other type of structure for this kind of analysis, for me as for the authors of my reference paper.

@Sextus Empiricus, I add these two plots, because I am not sure to understand what your sidenote remark means.

This first picture is the amount of birds (by abundance) seen for a julian day throughout the whole dataset. As you can see, we clearly note the spring and the autumn migration signals. There is also a good amount of data around christmas because of the christmas bird count, that I get rid off before analysis because it would be a bias.

This second picture is the amount of observations per julian day throughout the dataset. You can see the same migration signals.

I feel it is interesting to study both (I can be wrong here), even if I tend to agree with you when you say that the abundance is the most reliable. But it is not likely to have an impact on the quantile regression itself, does it? I am using citizen science data, based on birdwatching. The alternative would be to use data from a bird ringing station (as done in Karlsson et. al.), where the data collection is described as followed:

an effort is made most days within a fixed number of dates each year (called a season) to capture birds using mist nets at fixed locations. The number of mist nets vary depending on the weather, to ensure bird welfare, as the birds may be harmed under certain weather conditions. The nets are active for at least four hours a day, but some days they might be active for longer, and some days for a shorter time, depending on the weather. The varying effort to capture birds poses some challenges, that we will treat in a simple manner in this paper, since it is not of primary interest. If the effort is extended for more than four hours any given day, a larger portion of all the birds that could possibly be captured that day is expected to be captured. This may cause the resulting arrival distribution to be biased towards days with a larger effort. As a simple corrective measure, we only include birds captured in the first four hours of effort any given day. Ideally, one would however include data on all captured birds in phenological analyses and correct for the varying sampling effort. We will work under the assumption that the number of ringed birds a given day is proportional to the number of birds migrating past the study site the preceeding night.

• Sidenote (not the issue of the question): In the final plot you have time versus julian day. It looks like you have only single points per day and you model julian_day~time. Wouldn't it be interesting to use multiple datapoints for each day according to the number of birds seen? Currently it seems like the number of birds is not relevant. If you have observed hundreds of birds in winter but only on a few days and a few birds in spring but on many days, then this is being seen in your graph as a lot of points in spring and only a few points in winter. Apr 21 at 15:48
• You are trying to say that I should use the abundance of birds instead of the number of observations, right? As explained in my post update, I did it, but it increased a lot the size of my dataset which was not really easy to manipulate. I definitely think that using the abundance is the right way to do (for the reasons that you mentioned), but I thought that I would have an interest to try both, because what I am trying to capture is the signal of migration to study the shifts. Apr 21 at 15:57
• And I first thought that using observations on a long time period was good enough to capture this signal (when I plot it, I have a similar distribution in time compared to abundance). But when I implemented the abundance for my models, I had the same monotonous issue in the methodology. In Karlsson et. al. they also use abundance of birds. Apr 21 at 15:58
• Also, I do not have necessarily a single point per day; for example, I go out, see 10 birds in my garden on the 5th february 1980, and you do your own count the same day in your own garden, you see 3 birds, you will also provide an observation point to the dataset. It will just be recorded as 2 observation points instead of 10+3=13 birds (in the observation-based dataset, it will appear as 2 rows on this day whereas in the abundance-based dataset, you would be able to see 13 rows at this date (1 per bird). Apr 21 at 16:02
• Does it happen that you go out and do not see a bird? Otherwise your plot is just an expression how often you went out to see birds and record it. It does not reflect how many birds are around (or only indirectly, possibly you go bird-watching and record it more likely if you saw many birds around). Apr 21 at 17:38

I feel that it is not right if my intercet is not monotonous, so my questions are: what does it mean if my intercept is not monotonous?

The lines created by quantile regression should not cross each other. But if you extrapolate the curves then this might happen.

So what it means when your intercept is not monotonous is that your intercept falls in a range where extrapolation makes little sense.

In your last image you could plot on a scale where time extends to $$time = 0$$ (the intercept with the y-axis) and then you get that some of those lines cross.

The example below shows that this may occur due to random variations in the estimates of the lines. When the slope varies a little bit then this will have a large influence on the intercept.

The same principle may also occur, not due to random variations, but when you have a relationship that doesn't extrapolate well. For instance if the true relationship is not a linear function and the probability distribution is not having a constant variance/structure.

library(quantreg)

set.seed(1)
x = seq(10,20,0.05)
l = length(x)
y = 1+1*x+rnorm(l,0,3)
plot(x,y, pch = 21, col = 1, bg = 1, cex = 0.7, xlim = c(0,22), ylim = c(0,25))

xs = 0:20
for (tau_i in seq(0.1,0.9,0.1)) {
mod <- rq(y ~ x, tau = tau_i)
lines(xs, mod$$coefficients[1] + mod$$coefficients[2]*xs)
text(20, mod$$coefficients[1] + mod$$coefficients[2]*20, tau_i, pos =4, cex = 0.7)
}

set.seed(1)
x = seq(1,4,0.01)
l = length(x)
y = 1+1*x+rnorm(l,0,x-0.5)
plot(x,y, pch = 21, col = 1, bg = 1, cex = 0.7, xlim = c(0,4.2), ylim = c(-5,15))

xs = 0:4
for (tau_i in seq(0.1,0.9,0.1)) {
mod <- rq(y ~ x, tau = tau_i)
lines(xs, mod$$coefficients[1] + mod$$coefficients[2]*xs)
text(4, mod$$coefficients[1] + mod$$coefficients[2]*4, tau_i, pos =4, cex = 0.7)
}

• Thank you very much for your help. I think that I get what you are explaining, and I believe that it is called the Crossing Problem, which is why we want the intercept to be monotonous, if else a higher number of points may be captured by a lower quantile line as compared to the higher quantile line. I also read in Solution to the Non-Monotonicity and Crossing Problems in Quantile Regression (Saleh & Saleh (2021), source: arxiv.org/pdf/2111.04805.pdf) that "two lines may cross but the estimated model can be still be viewed as acceptable if the crossing occurs outside the convex hull". Apr 21 at 15:15
• To sum up, my thoughts right now are: 1) Try to jitter a little bit to see how the behaviour of my intercept evolves; 2) Try to implement the solution provided in Solution to the Non-Monotonicity and Crossing Problems in Quantile Regression and build their function in R; 3) Question why am I having this issue while it is not mentionned at all in my first reference, Quantile regression for phenological analyses (Karlsson et al.). This last point is what worries me the most... Apr 21 at 15:21
• @Recology to be honest, I had not yet viewed that article. I will look at it and see if I can add more. For the moment I am a bit surprised that this lack of monotonicity is considered a problem. The value of the intercept does not always have a strong physical meaning and it depends on the arbitrary location of the y-axis. In every situation where the slope is not a constant function of time (like in the figures that you have in the quotes) it is possible to choose a point where the intercept will have this crossing "problem". Related: stats.stackexchange.com/a/425547/164061 Apr 21 at 15:35