Theory understanding behind quantile regression

Question

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:

observation	year	julian_day	decade	time
1	y(1)	j(1)	d(1)	t(1)
...	...	...	...	...
n	y(n)	j(n)	d(n)	t(n)

2. 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= .)

3. 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...

4. 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.

I hope that this answer to your question.

Answer 1

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