#  A Comparison Between Quantile Regression and Linear Regression on Empirical
Quantiles for Phenological Analysis in Migratory Response to Climate Change
Running title: Quantile regression for phenological analyses

Måns Karlsson
Ola Hössjer Department of Mathematics, Stockholm University, 106 91 Stockholm,
Sweden<EMAIL_ADDRESS>of Mathematics, Stockholm University, 106
91 Stockholm, Sweden

###### Abstract

It is well established that migratory birds in general have advanced their
arrival times in spring, and in this paper we investigate potential ways of
enhancing the level of detail in future phenological analyses. We perform
single as well as multiple species analyses, using linear models on empirical
quantiles, non-parametric quantile regression and likelihood-based parametric
quantile regression with asymmetric Laplace distributed error terms. We
conclude that non-parametric quantile regression appears most suited for
single as well as multiple species analyses.

Keywords: Phenology, quantile regression, mixed effects, arrival times, linear
regression, bird observatory.

##  1 Introduction

Analysis of “the timing of seasonal activities of animals and plants” is
termed phenology (Walther et al., 2002), and it is a well studied topic. For
instance, in the meta analysis of Usui et al. (2017), 73 different studies
analysing bird phenology by means of linear regression were used to
investigate the phylogenetic signal in the response to climate change. The
problem can be approached in various ways, such as investigating if there is a
connection between migration timing change and variations of some climate
index (Jonzén et al., 2006), or if there is a trend in the location shift of
the arrival distribution (Lehikoinen et al., 2019). These two approaches were
used in the metastudy of Usui et al. (2017).

Phenological studies use data collected at one or several locations and may
cover one or multiple species. When birds are studied, the data used is often
collected in a systematic manner at bird observatories, since these
organisations tend to have the long running series of data needed for
phenological analyses (see Knudsen et al. (2007) for an in depth study). In
more recent years, citizen science data has been used as well (Mayer, 2010).
The sample unit of the data is usually individual birds, and these might be
trapped and ringed or just observed and registered without trapping, cf.
Lehikoinen et al. (2019) for a study where both of these data collection
methods were combined.

Associated with each bird is a varying number of covariates. Every observation
is associated with an arrival or departure date and almost always a species
information, which can be used as a block factor in multiple species analyses.
Ringed birds are also often registered to the hour of the day, and, if
possible, the age and sex of a bird is determined as well. Other biometrics
are sometimes collected, at least if the bird is captured.

The main interest in the studies of Jonzén et al. (2006) and Lehikoinen et al.
(2019), as well as in this study, concerns changes over time in the
distribution of arrival times. We will use year as the covariate of primary
interest, as the interpretation of the associated coefficient will tell us
what change has occured over time. Other covariates may also be used, such as
the age and sex of birds or other biometrics (Aharon-Rotman et al., 2021). The
empirical arrival distribution of any particular year can be visualized as a
step function with variable step sizes. The step sizes depend on the regional
and local weather at the study site, and a study of the weather effects for
fall migration of European Robins has been conducted by Karlsson (2014). In
this study we do not include weather covariates though.

Knudsen et al. (2007) attempted two ways of smoothing these yearly arrival
distributions, before fitting a model to data, and they concluded that
methodological advances were needed. In this paper we will not perform any
smoothing of this kind before fitting models to data.

The purpose of this article is to compare linear models for empirical
quantiles with linear quantile models without (Koenker, 2005) or with (Geraci
and Bottai, 2014) random effects. The latter model, which essentially is a
mixed effects model for quantile regression, has recently been employed for
phenological analyses (Aharon-Rotman et al., 2021). The remainder of our paper
is structured as follows. In Section 2, fixed effect models are first
introduced for single species data and then fixed and mixed effects models are
presented for multiple species data. Section 3 contains a presentation of the
data we will analyse, Section 4 provides a summary of the results of the
analysis, and finally in Section 5 we discuss potential improvements and
extensions. Some mathematical details are gathered in the supplementary
material.

##  2 Statistical approaches

This section contains an overview of the statistical methods used in this
paper for phenological analysis. We will cover the empirical quantiles linear
models used in Jonzén et al. (2006) and Lehikoinen et al. (2019),
nonparametric quantile regression (Koenker, 2005) and a likelihood-based
quantile regression approach using the asymmetric Laplace distribution (Geraci
and Bottai, 2014). These methods are described for singular species as well as
for multiple species analyses.

###  2.1 Single species models

Throughout this paper, the response variable will be julian day. For single
species models, $y_{j}$ will be the julian arrival day of individual
$j\in\\{1,\ldots,n\\}$. With each observation comes a covariate vector
$x_{j}=\left(1,t_{j},x_{j1},\ldots,x_{jp}\right)$, consisting of an intercept
term $1$, the year $t_{j}$ that individual $j$ was recorded, and $p$
additional covariates. Typical covariates are the binary age (juvenile or
adult) and sex (female or male), but continuous covariates such as wing length
may also be used (Aharon-Rotman et al., 2021). We assume that $t_{j}$ takes
values in ${\cal T}=\\{1,\ldots,T\\}$ for all $j$ and that
$x_{j1},\ldots,x_{jp}$ take values in ${\cal X}_{1},\ldots,{\cal X}_{p}$
respectively.

####  2.1.1 Empirical quantile linear models

The empirical quantile (eq) model does not take the raw julian day $y_{j}$ as
response, but rather an empirical quantile of the julian days of a subset of
observations. Let ${\cal X}={\cal X}_{1}\times\ldots\times{\cal X}_{p}$. For
each pair $(x,t)\in{\cal X}\times{\cal T}$ we extract the set of observations

${\cal
Y}_{(x,t)}=\left\\{y_{j}:1\left(t_{j}=t\land(x_{j1},\ldots,x_{jp})=x\right)\right\\},$
(1)

where $1(\cdot)$ is an indicator function that equals 1 if the condition
within the brackets is true, and 0 otherwise. We also create the $T\cdot|{\cal
X}|\times(2+p)$ matrix $X^{\text{eq}}$ by stacking all vectors $(1,t,x)$ on
top of each other. Let $\tau\in(0,1)$ be a quantile. For each of our sets
${\cal Y}_{(x,t)}$ we let $\hat{F}_{(x,t)}$ be the empirical distribution
function formed by the elements of this set, and define the corresponding
empirical quantile

$\hat{Q}_{(x,t)}(\tau)=\inf_{y}\left\\{y\in{\cal
Y}_{(x,t)}:\hat{F}_{(x,t)}(y)\geq\tau\right\\},$ (2)

conditional on ${\cal Y}_{(x,t)}\neq\emptyset$, and stack these quantiles, for
all $(x,t)\in{\cal T}\times{\cal X}$, into the vector $Y^{\text{eq}}(\tau)$.
We then formulate the linear model

$\displaystyle
Y^{\text{eq}}(\tau)=X^{\text{eq}}\beta^{\text{eq}}(\tau)+\varepsilon^{\text{eq}}(\tau)$
(3)

where (omitting $\tau$) the column vector
$\beta^{\text{eq}}=(\beta^{\text{eq}}_{0},\beta^{\text{eq}}_{t},\beta^{\text{eq}}_{1},\ldots,\beta^{\text{eq}}_{p})^{\top}$
contains the intercept $\beta^{\text{eq}}_{0}$, slope on year
$\beta^{\text{eq}}_{t}$ and the $p$ covariate effects
$\beta^{\text{eq}}_{1},\ldots,\beta^{\text{eq}}_{p}$, whereas $\top$ refers to
matrix transposition. The vector of error terms in
$\varepsilon^{\text{eq}}(\tau)$ is assumed to have a multivariate normal
distribution $N\left(0,\sigma_{\text{eq}}^{2}(\tau)I_{T\cdot|{\cal
X}|}\right)$, where $I_{T\cdot|{\cal X}|}$ is the identity matrix of rank
$T\cdot|{\cal X}|$. This approach gives each combination of covariate values
and year the same weight. However, if we want to give each recorded bird the
same weight, some adjustements are needed. Omitting any dependence on $\tau$
in the notation, the log-likelihood for the model in (3) is

$l\left(\beta^{\text{eq}},\sigma_{\text{eq}}^{2}\mid
Y^{\text{eq}},X^{\text{eq}}\right)=\sum_{(x,t)\in{\cal T}\times{\cal X}}\log
f\left(\hat{Q}_{(x,t)}(\tau)\mid(X^{\text{eq}}\beta^{\text{eq}})_{(x,t)},\sigma_{\text{eq}}^{2}\right)$
(4)

where $f(\cdot\mid m,\sigma^{2})$ is the Gaussian density with mean $m$ and
variance $\sigma^{2}$. Let $w_{(x,t)}$ denote a weight given to each empirical
quantile $\hat{Q}_{(x,t)}(\tau)$. If we let $w_{(x,t)}=|{\cal Y}_{(x,t)}|$, we
will achieve equal weighting of each bird by reweighting the log-likelihood
(4) as

$l_{w}\left(\beta^{\text{eq}},\sigma_{\text{eq}}^{2}\mid
Y^{\text{eq}},X^{\text{eq}}\right)=\sum_{(x,t)\in{\cal T}\times{\cal
X}}w_{(x,t)}\log
f\left(\hat{Q}_{(x,t)}(\tau)\mid(X^{\text{eq}}\beta^{\text{eq}})_{(x,t)},\sigma_{\text{eq}}^{2}\right).$
(5)

As is well known, the objective function (5) is continous, concave and twice
differentiable, and thus it is easily optimized in order to find the MLE of
$\beta^{\text{eq}}$ and $\sigma_{\text{eq}}^{2}$. This model fitting can then
be repeated for any $\tau$ of interest.

####  2.1.2 Non-parametric quantile regression

For a deeper treatment of quantile regression (qr), we refer the reader to the
standard litterature of Koenker (2005), which is the primary source for the
material that is presented in this section.

Instead of a linear predictor for the mean of the response distribution, a
quantile regression model has a linear predictor for a conditional quantile of
the response distribution. As before, let $y_{j}$ be the observed response for
individual $j$, and stack all the covariate vectors
$x_{j}=(1,t_{j},x_{j1},\ldots,x_{jp})$ into the $n\times(p+2)$-matrix $X$. Let
the conditional distribution function of the random variable $Y_{j}$
corresponding to its observed value $y_{j}$ be

$F_{Y_{j}\mid x_{j}}(y)=\mathbb{P}(Y_{j}\leq y\mid
x_{j}),\quad-\infty<y<\infty.$ (6)

For each $0<\tau<1$ the inverse

$Q(\tau\mid x_{j})=\inf\\{y;F_{Y_{j}\mid x_{j}}(y)\geq\tau\\}$ (7)

of (6) is the conditional quantile function, and we may construct models of
the form

$Q\left(\tau\mid
X\right)=X\beta^{\text{qr}}(\tau)+\varepsilon^{\text{qr}}(\tau)$ (8)

where (omitting $\tau$)
$\beta^{\text{qr}}=(\beta^{\text{qr}}_{0},\beta^{\text{qr}}_{t},\beta^{\text{qr}}_{1},\ldots,\beta^{\text{qr}}_{p})^{\top}$
is analogous to $\beta^{\text{eq}}$. No parametric assumptions are made about
the error terms $\varepsilon^{\text{qr}}(\tau)$.

The objective function in (8) is the loss function

$\sum_{j=1}^{n}\rho_{\tau}(y_{j}-x_{j}b)$ (9)

where

$\rho_{\tau}(\nu)=\nu(\tau-1(\nu<0))$ (10)

(Koenker, 2005, p. 5). We thus choose $\tau$ and estimate the regression
parameters by finding

$\beta^{\text{qr}}(\tau)=\arg\min_{b\in{\mathbb{R}}^{2+p}}\sum_{j=1}^{n}\rho_{\tau}(y_{j}-x_{j}b).$
(11)

As opposed to (5), the objective function in (11) is not differentiable
everywhere. In more detail, the derivative of the objective function has
discontinuities along hyperplanes $\\{b;\,y_{j}=x_{j}b\\}$ for $j=1,\ldots,n$.
Non-gradient based numerical estimation procedures are however available, and
we use the Frisch-Newton interior point method (Portnoy and Koenker, 1997)
provided in the R package quantreg (Koenker, 2021).

As with the empirical quantile models, we may fit quantile regression models
to any number of quantiles $\tau$. Since (11) applies directly to raw data
$\\{(x_{j},y_{j});\,j=1,\ldots,n\\}$ on the sample unit level, no reweighting
is needed to give each bird the same weight.

####  2.1.3 Linear quantile models

The third approach to modelling a single species is the linear quantile model
(lqm). It is based on the assumption that the error terms in (8) are
independent and identically distributed, and follow an asymmetric Laplace (AL)
distribution. A random variable $\Gamma\sim AL(\mu,\sigma,\tau)$ if it has
probability density function

$f_{\gamma}(\mu,\sigma,\tau)=\frac{\tau(1-\tau)}{\sigma}\exp\left\\{-\frac{1}{\sigma}\rho_{\tau}(\gamma-\mu)\right\\},$
(12)

where $\mu\in(-\infty,\infty)$ is the location parameter, $\sigma>0$ is the
scale parameter and $\tau\in(0,1)$ is the asymmetry parameter (Hinkley and
Revankar, 1977). Note that

$\displaystyle\mathbb{P}\left(\Gamma\leq\mu\right)$
$\displaystyle=\int_{-\infty}^{\mu}\frac{\tau(1-\tau)}{\sigma}\exp\left\\{-\frac{1}{\sigma}\rho_{\tau}(\gamma-\mu)\right\\}\,\mathrm{d}\gamma$
(13)
$\displaystyle=\frac{\tau(1-\tau)}{\sigma}\int_{-\infty}^{\mu}\exp\left\\{(\gamma-\mu)\frac{1-\tau}{\sigma}\right\\}\,\mathrm{d}\gamma$
(14)
$\displaystyle=\tau\left[\exp\left\\{(\gamma-\mu)\frac{1-\tau}{\sigma}\right\\}\right]_{\gamma=-\infty}^{\gamma=\mu}$
(15) $\displaystyle=\tau$ (16)

meaning that $\mu$ is the $\tau$:th quantile of the response variable
distribution, when it is assumed to be asymmetric Laplacian. Therefore,
inference about the location parameter $\mu$ will be equivalent to inference
about the $\tau$:th quantile of the response variable distribution $F_{Y}$.
For further details, see e.g. Koenker and Machado (1999).

Assume for simplicity that the error terms in (8) have a scale parameter
$\sigma_{\text{lqm}}=1$. Then the log-likelihood to maximize under the
asymmetric Laplace assumption is

$l(\beta^{\text{lqm}},\sigma_{\text{lqm}}\mid\tau,Y,X)=n\log(\tau(1-\tau))-\sum_{j=1}^{n}\rho_{\tau}(y_{j}-x_{j}\beta^{\text{lqm}}(\tau)).$
(17)

A gradient of the objective function (17) is defined by Bottai et al. (2015),
and it is used in a gradient-based search algorithm in order to find an
estimate of $\beta^{\text{lqm}}$. In a large sample simulation study of the
quantreg-package, Bottai et al. (2015) found computational advantages of the
gradient search method when compared with the abovementioned Frisch-Newton
approach. They also noted that the bias (approximated as the difference
between the estimated and true parameter values, average over many samples)
was practically zero for both the gradient search based method and the Frisch-
Newton method. They did not comment on the largest observed difference between
true parameter value and parameter estimate, but the largest observed relative
absolute difference in the parameter estimates of the two methods was 0.12.

In conclusion, since minimizing the objective function in (11) with respect to
the regression parameter is equivalent to maximizing the log likelihood in
(17), the linear quantile model provides the same parameter estimates as the
quantile regression model of the previous section, if the same numerical
estimation algorithm is used for parameter optimization. In Section 4.1 we
will investigate how the choice of optimization method affects parameteter
estimates.

###  2.2 Multiple species models

To incorporate several species in a phenological analysis, we can either
include species as a categorical predictor, and add interaction effects
between this and other covariates accordingly, or we can consider species a
random block factor and make use of mixed effects models (Lehikoinen et al.,
2019; Jonzén et al., 2006). We will consider the categorical predictor
approach for nonparameteric quantile regression (qr), and the block factor
approach for empirical quantile mixed models (meq) as well as for the Linear
Quantile Mixed Model (lqmm) approach (Geraci and Bottai, 2014).

In the case of the presence of a block factor (i.e. species) in data, the
response is $y_{ij}$ where $i=1,\ldots,M$ denotes which species the
observation belongs to, and $j=1,\ldots,n_{i}$ index the observations within
species $i$. In total, we have $\sum_{i=1}^{M}n_{i}=N$ observations.
Analogously the covariate vector associated with each observation is now
$x_{ij}=\left(1,t_{ij},x_{ij1},\ldots,x_{ijp}\right)$. We introduce the
vectors $z_{ij}=\left(1,t_{ij},z_{ij1},\ldots,z_{ijq}\right)$, for
$i=1,\ldots,M$ and $j=1,\ldots,n_{i}$. These vectors $z_{ij}$ will be used
below to introduce a vector $u$ with $q+2$ random effects, for the meq and the
lqmm models. The first two components of this vector correspond to a random
intercept and slope for each species.

####  2.2.1 Empirical quantile linear mixed effects models

With the introduction of the block factor species, we repeat the process of
computing empirical quantiles in Section 2.1.1 for each species $i\in{\cal
M}=\\{1,\ldots,M\\}$ and set of years ${\cal T}_{i}$. In order to simplify our
exposition it is assumed in this section that $q=0$ and $z_{ij}=(1,t_{ij})$
for all birds $(i,j)$. The multispecies empirical quantile (meq) model of
species $i$ then takes the form

$Y^{\text{meq}}_{i}(\tau)=X^{\text{meq}}_{i}\beta^{\text{meq}}+Z^{\text{meq}}_{i}u^{\text{meq}}_{i}+\varepsilon^{\text{meq}}_{i}(\tau)$
(18)

where the $\lvert{\cal T}_{i}\rvert\cdot\lvert{\cal X}_{i}\rvert\times(2+p)$
matrix $X^{\text{meq}}_{i}$ is created by stacking the vectors
$(1,t_{i},x_{i})$, with $x_{i}\in{\cal X}_{i}$, where ${\cal X}_{i}$ is the
set of values $(x_{ij1},\ldots,x_{ijp})$ can take, whereas
$Z^{\text{meq}}_{i}$ consists of the first two columns of
$X^{\text{meq}}_{i}$. The random effects vector $u^{\text{meq}}_{i}$ has a
two-dimensional Gaussian distribution with expected value $(0,0)^{\top}$ and
covariance matrix $\Psi$. The likelihood contribution of species $i$ (cf.
Section 2.1.1) is thus obtained by integrating over the unobserved random
effects, i.e.

$\begin{split}L_{i}\left(\Psi,\beta,\sigma\mid y\right)&=\\\
\prod_{(x,t)\in{\cal X}_{i}\times{\cal
T}_{i}}&\int_{\mathbb{R}^{q}}f\left(\hat{Q}_{(x,t)}(\tau)\mid\left(X^{\text{meq}}_{i}\beta^{\text{meq}}+Z^{\text{meq}}_{i}u_{i}\right)_{(x,t)},\sigma\right)f(u_{i}\mid\Psi)\,\mathrm{d}u_{i}.\end{split}$
(19)

To weight each bird equally, we again follow Section 2.1.1 and create weights
$w_{(i,t)}$ based on the distribution of observations across year and other
covariates. The reweighted log likelihood contribution of species $i$ is then

$l_{i}\left(\Psi,\beta,\sigma\mid y\right)=\sum_{(x,t)\in{\cal
X}_{i}\times{\cal
T}_{i}}w_{(x,t)}\log\left(\int_{\mathbb{R}^{q}}f\left(\hat{Q}_{(x,t)}(\tau)\mid\ldots\right)f(u_{i}\mid\Psi)\,\mathrm{d}u_{i}\right),$
(20)

with the same integration of random effects as in (19). By summing the
reweighted log likelihood over all species we obtain the log likelihood

$l\left(\Psi,\beta,\sigma\mid
y\right)=\sum_{i=1}^{M}l_{i}\left(\Psi,\beta,\sigma\mid y\right)$ (21)

of the full dataset. We will fit the linear mixed effects models using the R
(R Core Team, 2021) package lme4 package (Bates et al., 2015). Details of the
model specification and computational approaches are given in the lme4-package
vignette of Bates (2014).

####  2.2.2 Linear quantile mixed models

The framework called linear quantile mixed models (lqmm) was proposed by
Geraci and Bottai (2014) as an extension of the random intercept quantile
regression model of Geraci and Bottai (2007), and it is implemented in the
R-package lqmm (Geraci et al., 2014). The aim of lqmm is to provide an
analogue of linear mixed effects models in the setting of quantile regression,
exploiting the connection between the asymmetric Laplace likelihood and the
objective function (9) of a nonparametric quantile regression model (Koenker
and Machado, 1999).

In an lqmm for quantile $\tau$, the linear predictor is

$\mu_{ij}(\tau)=x_{ij}\beta(\tau)+z_{ij}u_{i}$ (22)

and it is assumed that

$y_{ij}\mid u_{i}\sim
AL\left(\mu_{ij}(\tau),\sigma_{\varepsilon}(\tau),\tau\right)$ (23)

independently for all $i,j$, with
$\varepsilon_{ij}(\tau)=y_{ij}-\mu_{ij}(\tau)$ for all $i$ and $j$. The random
effects vector $u_{i}$ is independent of
$\\{\varepsilon_{ij}(\tau)\\}_{j=1}^{n_{i}}$ for any $i\in\cal M$, and, for
our models, it is assumed to be multivariate Gaussian with covariance matrix
$\Psi(\tau)$. Omitting in the notation the dependence of parameters on $\tau$,
the marginal likelihood

$L\left(\beta,\sigma_{\varepsilon},\Psi\mid\left\\{y_{ij}\right\\}\right)=\prod_{i=1}^{M}\int_{\mathbb{R}^{q+2}}\prod_{j=1}^{n_{i}}f(y_{ij}\mid\mu_{ij},\sigma_{\varepsilon},\tau)p\left(u_{i}\mid\Psi\right)\,\mathrm{d}u_{i}$
(24)

is obtained by integrating out the random effects of each species. In Geraci
and Bottai (2014) the integral of each term of (24) is approximated
numerically with a Gaussian quadrature procedure.

Geraci and Bottai (2014) prove that the likelihood (24) is log-concave. They
also make use of the iterative gradient-based optimization algorithm of Bottai
et al. (2015) in order to find the argmax of (24). In this context they also
point out that derivative free methods are also viable for maximizing (24),
namely the Nelder-Mead derivative free optimization (Nelder and Mead, 1965) in
the optim function in R. Convergence is, in theory, guaranteed by the log-
concavity of the likelihood, but when maximizing a numerical approximation of
(2.2.2), log-concavity is not guaranteed and parameter estimates may depend on
starting values, as mentioned in Geraci et al. (2014, p.14).

####  2.2.3 Non-parameteric quantile regression with interactions

Since nonparametric quantile regression (qr) does not include random effects,
we will instead include species as a fixed effect and let it interact with
year. The purpose is to compare the quantile regression estimates with
estimates obtained from the meq and lqmm models of Sections 2.2.1 and 2.2.2.
In particular, for the latter two approaches, we will study inference
conditional on predictions of the random effects $u_{i}$, $i=1,\ldots,M$.

Following Section 2.1.2, we let each species $i$ contribute with two
covariates, a species specific intercept parameter $\beta_{0i}$ and a species
specific slope parameter $\beta_{ti}$, which is interpreted as an interaction
parameter between species and year. In order not to overparametrize, species
$i=1$ will be included in the intercept $\beta_{0}$, and the overall slope on
year $\beta_{t}$. The model, in the special case $p=0$ and
$x_{ij}=(1,t_{ij})$, is thus

$Q(\tau\mid
x_{ij})=\beta^{\text{qr}}_{0}+\beta^{\text{qr}}_{t}t_{ij}+1(i>1)\left(\beta^{\text{qr}}_{0i}+\beta^{\text{qr}}_{ti}t_{kj}\right)+\varepsilon^{\text{qr}}_{ij},$
(25)

for $i=1,\ldots,n$, keeping dependence on $\tau$ only on the left-hand side of
(25). Notice that the interpretation of $\beta^{\text{qr}}_{0i}$ and
$\beta^{\text{qr}}_{ti}$ is the difference in intercept and slope of species
$i>1$ compared to species $1$.

###  2.3 Predicting $u$

Inference in mixed models can be marginal, i.e. looking at the fixed effects
estimates $\hat{\beta}$ only, or conditional, i.e. looking at the species
specific estimates of quantities obtained by conditioning on the random
effects $u=\left(u_{1},\ldots,u_{M}\right)$. Conditional inference requires
predictions $\hat{u}$ of the random effects, which we now present for the
linear mixed model of Section 2.2.1 and the lqmm of Section 2.2.2.

####  2.3.1 Linear mixed models

We use the lmer function of the lme4-package (Bates et al., 2015) to fit
linear mixed effects models. The predicted random effects $\hat{u}$ are
extracted with the function ranef, and the value obtained is the mode of the
conditional density of the random effects, conditionally on the response and
the estimates of the other model parameters (Bates, 2014, sec. 5.1.4).

####  2.3.2 Linear quantile mixed models

An adaption of the ranef function is provided in the lqmm-package, named
ranef.lqmm. We will briefly describe the random effects output of this
function and how we adapted the prediction of $u$ to be more computationally
efficient for our models. The full details of our adaption are provided in a
supplementary document.

Given the lqmm of Section 2.2.2, we let $Z$ be the random effects design
matrix, $X$ the fixed effects design matrix and $y$ the response vector. The
function ranef.lqmm gives the estimated best linear predictor

$\hat{u}(\tau)=\hat{\Psi}(\tau)Z^{\top}\hat{\Sigma}^{-1}\left\\{y-X\hat{\beta}(\tau)-\hat{\mathbb{E}}\left[\varepsilon(\tau)\right]\right\\},$
(26)

of the random effects in an lqmm (Geraci and Bottai, 2014, eq. 12), where
$\hat{u}(\tau)$ is an $M\times(q+2)$ matrix of the estimated random effects
and $\hat{\Sigma}$ is an estimate of

$\Sigma=\mbox{Cov}(Y)=Z\Psi(\tau)Z^{\top}+\mbox{Cov}(\varepsilon(\tau)),$ (27)

the $N\times N$ block-diagonal covariance matrix of the response vector $Y$,
of which $y$ is an observation. Note that
$\mbox{Cov}(\varepsilon(\tau))=\sigma_{\varepsilon}(\tau)I_{N}$, where $I_{N}$
is the identity matrix of rank $N$ and $\sigma_{\varepsilon}(\tau)$ is the
scale parameter of the asymmetric Laplace distribution for quantile $\tau$. We
thus need to invert $\hat{\Sigma}$ to compute $\hat{u}(\tau)$, and since
ranef.lqmm tries to invert $\Sigma$ head-on by a call to solve, this is
sometimes unfeasible for large data sets on consumer grade computers.
Therefore, we adapted (26)-(27) by inverting $\Sigma$ in blocks, and thereby
made prediction of $u$ a lot faster. The full details of the adapted
computational approach are given in a supplementary document. We tested our
adapted matrix inversion against the ranef.lqmm output for a subset of data,
and the two approaches produced the same results.

##  3 Data

The phenological data used in this paper has been collected and kindly
provided for use by the Falsterbo Bird observatory. This section gives an
overview of the data collection method and how data was filtered.

###  3.1 Data collection

The Falsterbo Bird Observatory (55∘23’01.6”N, 12∘49’00.9”E) uses a method for
systematic trapping and ringing of birds described in Roos and Karlsson
(1981). In short, 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. The assumption is not fully verified. The correlation
between nocturnal migratory activity and the number of ringed birds during
autumn migation in Falsterbo has been studied in Zehnder and Karlsson (2001).
They found a positive correlation under the assumption of weather affecting
the number of trapped birds in the same way as it effects the migratory
activity. Since we work with spring migration, patterns and behavior might be
different. The general perception from field work is that most days with
weather such that the welfare of caught birds is risked, correlate negatively
with the number of birds. This means that a smaller number of birds is
expected to be captured these days. However, some days such weather conditions
rather cause a high number of birds to congregate in Falsterbo.

###  3.2 Data selection and overview

In total, $151\,992$ birds of 89 species have been recorded during the spring
seasons of 1980-2019. We filter these data to remove birds recorded after the
first four hours of effort, and species with less than 150 recorded birds in
total, which yields a data set of $137\,870$ birds distributed over 38
species. From these, we select 25 species that can be assumed to mostly
consist of migratory individuals but with few locally breeding birds mixed in.
This data set is used for all multiple species models of Section 2.2.

The data is very unbalanced with respect to species. In Figure 1 illustrates
the distribution of number of observations of each species, each year. There
are two species that are very numerous, and the other species usually have
much fewer birds each year. This could be an argument to rebalance data so
that each species has the same weight in a mixed effects model. However, the
fixed effects would not be interpretable then as the overall effects of birds
in Falsterbo, but rather the effects when weighting up scarcely recorded
species.

Figure 1: A color-coded barplot of the number of birds of each species each
year. The plot is intended to illustrate the numerical dominance of European
robin and Willow warbler in particular, and that year totals are quite
consistent over time. Specific numbers are available in the supplementary
material.

For the single species models, we chose the species Eurasian Blackcap (Sylvia
atricapilla), since it is easily classified in terms of species, age and sex.
Thus, we have highly reliable data, and it is possible to incorporate the
covariates age and sex in the single species models. The Eurasian Blackcap is
quite a common choice of species for phenological analysis, and an interesting
modeling approach can be found in Aharon-Rotman et al. (2021) for this
species, using lqmm with year as block factor and fixed effect. The
distribution of birds across age and sex levels are shown in Table 1, whereas
the number of birds recorded each year is illustrated in Figure 2. Another
single species phenological analysis can be found in Hössjer and Karlsson
(2021), with an analysis of Common Redstart Phoenicurus phoenicurus. In that
article we investigate not only changes in location of the arrival
distribution, but also changes in scale, skewness and kurtosis, and how
covariates affect these.

| Female | Male
---|---|---
Juvenile | 567 | 786
Adult | 406 | 444
Table 1: Age-sex distribution of the Eurasian blackcap data. There is a larger
number of juveniles and males caught, but the imbalance is not particularly
extreme. Figure 2: The distribution of Eurasian blackcaps across years.

##  4 Results

First, we look at the single species Eurasian blackcap and compare two ways of
estimating a quantile regression model with an empirical quantiles model, in
terms of paramater estimates and pointwise confidence interval coverage.
Thereafter, we fit models for multiple species, with species as a block
factor, and a random intercept and random slope on year for each species, or
as a categorical fixed effect, interacting with year. These models are
scrutinized with special regard to quantile crossing. For all fitted models,
year was centered around 2001.

###  4.1 Species specific modelling

####  4.1.1 Model setup

All models fitted to the Eurasian blackcap data set include the same
covariates, and (as covered in Section 2.1) differ in what is used as response
variable, and how the objective function is constructed. The general
specifications of the empirical quantiles (eq) and quantile regression (qr)
models are provided in equations (3) and (8) respectively. For all models
$p=2$, with $\beta_{0}$ the intercept, $\beta_{t}$ the slope on year $t$,
$\beta_{1}$ the effect of age $x_{1}$ and $\beta_{2}$ the effect of sex
$x_{2}$. The quantile is specified by $\tau$, $\varepsilon$ denotes the error
term and there are $n=2\,203$ observations. Models were fitted for
$\tau=0.01,0.02,\ldots,0.99$. Pointwise bootstrap based approximate confidence
intervals were computed for all regression parameters $\hat{\beta}$, for each
$\tau$. This was done by resampling pairs $(y,x)$ with replacement until 1000
resampled data sets of the same size as the original had been generated. Each
model was fitted to the same resampled data set, to get the most comparable
confidence intervals, which were constructed by extracting the $2.5\%$ and
$97.5\%$ quantiles from the resulting sample of parameter estimates.

####  4.1.2 Estimation methods

For all eq-models we used the maximum likelihood estimates. For the qr-models
we used the Frisch-Newton interior point method for optimizing (Portnoy and
Koenker, 1997) the objective function (11) in order to get our parameter
estimates. For the lqms, we used both the gradient based method of Bottai et
al. (2015), and Nelder-Mead interior point estimation (Nelder and Mead, 1965).
These two optimization approches yielded considerably different estimates. The
estimates obtained through Nelder-Mead optimization were very close to the
Frisch-Newton estimates for qr, which is to be expected due to the similarity
of the two objective functions. However, the gradient based estimates were
different from those obtained with the other two optimization methods, for
$\hat{\beta}_{1}(\tau)$ as well as for $\hat{\beta}_{2}(\tau)$. For this
reason, we chose to include the gradient based estimates in Figure 3, to
highlight the differences occuring with this optimization method, since the
estimates using Nelder-Mead were very close to the qr-estimates.

####  4.1.3 Parameter estimates

For each method we checked if the intercept estimates were monotonically
increasing in $\tau$. It can clearly be seen in Figure 3 that
$\hat{\beta}^{\text{lqm}}_{0}(\tau)$ is not monotone. Moreover, there is one
instance where $\hat{\beta}^{\text{qr}}_{0}$ is decreasing with $\tau$, but
since the decrease is of magnitude $10^{-14}$, we cannot rule out numerical
errors. For eq, the intercept estimates are always monotone in $\tau$. Since
year is centered around 2001, the interpretation of the intercept is the
estimated arrival distribution year 2001.

For the eq model, notice the lacking coverage of the confidence intervals in
the tails of the distribution when it comes to the intercept estimate
$\hat{\beta}^{\text{eq}}_{0}$. The estimates of the effect of age
$\hat{\beta}^{\text{eq}}_{1}$ also approaches the limits of the confidence
interval. Although some missing coverage is expected since the confidence band
consists of interpolated pointwise approximate 95% confidence intervals, the
large discrepancy between estimates and confidence intervals is not ideal.

Further inspection of Figure 3 reveals that the intercept estimates are
relatively close for all three methods, with $\hat{\beta}^{\text{qr}}_{0}$
having a bit heavier tails. The estimates
$(\hat{\beta}^{\text{qr}}_{t},\hat{\beta}^{\text{lqm}}_{t})$ of the effect of
year are also close across $\tau$-values, whereas
$\hat{\beta}^{\text{eq}}_{t}$ does not decrease as much in the two middle
quartiles. Conversely, the effects of age
$(\hat{\beta}^{\text{qr}}_{1},\hat{\beta}^{\text{eq}}_{1})$ and sex
$(\hat{\beta}^{\text{qr}}_{2},\hat{\beta}^{\text{eq}}_{2})$ are similar,
whereas $\hat{\beta}^{\text{lqm}}_{1}$ and $\hat{\beta}^{\text{lqm}}_{2}$
differ markedly from the former two methods in the lower and upper tertials of
the distribution. One can also notice near identical estimates of
$\hat{\beta}^{\text{lqm}}_{1}$ and $\hat{\beta}^{\text{lqm}}_{2}$ in the
lowest quartile. All covariates have a diminishing effect as we approach the
right tail of the distribution. This might be due to data containing some
local breeding birds, especially towards the end of the sampling window. Since
these do not originate from the migration process, they might showcase
different patterns.

Method | intercept | year | age | sex
---|---|---|---|---
eq | 3.358 | 0.164 | 3.900 | 3.832
lqm | 3.646 | 0.158 | 3.207 | 3.514
qr | 2.939 | 0.151 | 3.390 | 3.286
Table 2: The table presents the mean CI width of the single species model, for
each combination of method and parameter, when averaged over all values of
$\tau$ for which estimation was performed. The shortest mean CI is highlighted
with bold font.

####  4.1.4 Conclusions

Our overall assessment is that qr produces the most consistent and detailed
results, as well as having the overall tightest confidence intervals. As a
caution, we do not recommend gradient based optimization when including binary
covariates in the quantile regression models. The large deviations in the
parameter estimates compared to the bootstrapped estimates for the eq model
make us cautions using this method, although it otherwise seems to detect the
more or less same signal as the qr model.

Figure 3: The legend at the top follows the abbreviations already introduced
in the article. The solid line consists of the actual parameter estimates for
all quantiles, the dashed line is a cubic spline smoothing of all the
bootstrap estimates and the shaded band is the interpolated 95% pointwise
confidence interval based on quantiles of the bootstrap distribution. More
detailed comments are provided in Section 4.1.3.

###  4.2 Multiple species modelling

For the multiple species models we used the three approaches described in
Section 2.2. Since some of the species included lack age and sex data, no
other fixed effects than year were included. For the qr-model, species was a
fixed effect interacting with year. In order to ease the computational burden
we only fitted these models to 21 quantiles, namely
$(0.01,0.05,\ldots,0.95,0.99)$. For each quantile we predicted random effects
$\hat{u}_{i}^{\text{meq}}$ and $\hat{u}_{i}^{\text{lqmm}}$ of the meq and lqmm
respectively, for each species $i$.

For the lqmm we were unable to get the gradient based estimation algorithm
(Bottai et al., 2015) to work on anything but a small subset of the data.
Instead we used Nelder-Mead optimization (Nelder and Mead, 1965), implemented
in the optim-function in R, which can be called directly from the lqmm-
function. As for the single species case, the qr-model was estimated using the
Frisch-Newton method (Portnoy and Koenker, 1997). The meq-model was estimated
using REML (Laird and Ware, 1982), described in Section 3.4 of Bates et al.
(2015).

We computed

$\begin{split}\hat{\beta}_{i}^{\text{meq}}(\tau)&=\hat{\beta}^{\text{meq}}(\tau)+\hat{u}^{\text{meq}}_{i}(\tau)\\\
\hat{\beta}_{i}^{\text{lqmm}}(\tau)&=\hat{\beta}^{\text{lqmm}}(\tau)+\hat{u}^{\text{lqmm}}_{i}(\tau)\end{split}$
(28)

for $i=1,\ldots,M$ to get species specific intercepts $\hat{\beta}_{i0}(\tau)$
and slopes $\hat{\beta}_{it}(\tau)$ on year for each model. For the qr-model
there are only species specific effects, since species is a fixed effect
interacting with year. We thus compute $\hat{\beta}_{i}^{\text{qr}}$ by
summing the two coefficients on the right-hand side of (25), for the intercept
and effect of year respectively. In Figures 4 and 5 we have linearly
interpolated and plotted the intercepts and effects of year across all 21
quantiles.

####  4.2.1 Intercept

Figure 4: Intercept plots for the multiple species models. The grey lines
illustrate the interpolated species specific intercept estimates
$\hat{\beta}_{i0}^{\text{qr}}(\tau)$, $\hat{\beta}_{i0}^{\text{meq}}(\tau)$
and $\hat{\beta}_{i0}^{\text{lqmm}}(\tau)$ for $i=1,\ldots,N$. The black line
is the overall intercept estimate for the meq and lqmm methods, and for the
qr, it is the weighted mean
$N^{-1}\sum_{i=1}^{M}n_{i}\hat{\beta}_{i0}^{\text{lqmm}}(\tau)$ of the species
specific estimates.

Of the intercept estimates, $\hat{\beta}_{i0}^{\text{qr}}(\tau)$ was
monotonically increasing in $\tau$ for all $i$, which implies that the
weighted mean also is monotonically increasing. For the lqmm-models
$\hat{\beta}_{i0}^{\text{lqmm}}(\tau)$ increased monotonically for all $i$,
but $\hat{\beta}^{\text{lqmm}}(\tau)$ did not (visible in Figure 4 for several
quantiles). In view of (28), this indicates that the fixed and random
intercept effects of the lqmm-model are difficult to separate for some
quantiles $\tau$. For the meq model, there was one species for which
$\hat{\beta}_{i0}^{\text{meq}}(0.01)>\hat{\beta}_{i0}^{\text{meq}}(0.05)$,
otherwise $\hat{\beta}_{i0}^{\text{meq}}(\tau)$ increased monotonically in
$\tau$ for all $i$ and so did $\hat{\beta}^{\text{meq}}(\tau)$.

The species specific intercept estimates
$\hat{\beta}_{i0}^{\text{lqmm}}(\tau)$ of the lqmm stand out in that most of
them are close to parallel (since the predicted random intercept effect of
each species is virtually independent of $\tau$). This would imply that the
only difference between species’ arrival distributions is a location shift
(i.e. a temporal shift of the distribution). The estimates further stand out
in that they exceed the sampling limits at both ends of the sampling window.
Some species might indeed be migrating outside the sampling window, but
further examination is needed to rule out the possibility that the sampling
window limits are exceeded as a consequence of the difficulty of obtaining a
flexible, species-specific estimates of the arrival distribution.

Although $\hat{\beta}_{i0}^{\text{qr}}(\tau)$ and
$\hat{\beta}_{i0}^{\text{meq}}(\tau)$ are similar, there is a general
difference in that $\hat{\beta}_{i0}^{\text{qr}}(\tau)$ has more steeply
sloping tails, especially in the lower quantiles. The overall conclusion of
the crossing lines is that the arrival distributions of different species
differ in further aspects than location (Oja, 1981).

####  4.2.2 Year

Figure 5: Plots of $\hat{\beta}_{t}^{\text{lqmm}}(\tau)$,
$\hat{\beta}_{t}^{\text{meq}}(\tau)$ and $\hat{\beta}_{t}^{\text{qr}}(\tau)$
as a function of $\tau$. Notice the sharp jumps in
$\hat{\beta}_{t}^{\text{lqmm}}(\tau)$, whereas
$\hat{\beta}_{t}^{\text{meq}}(\tau)$ and $\hat{\beta}_{t}^{\text{qr}}(\tau)$
follow each other closely.

In Figure 5 the sharply oscillating behavior of
$\hat{\beta}_{t}^{\text{lqmm}}(\tau)$, as a function of $\tau$, is clearly
visible. Combined with the flat lines of
$\hat{\beta}_{it}^{\text{lqmm}}(\tau)$ in Figure 6, the conclusion is that
$\beta_{it}^{\text{lqmm}}(\tau)$ is possible to estimate, but the separation
of $\beta^{\text{lqmm}}$ and $u_{i}$ is problematic for each quantile $\tau$.
Since $\hat{\beta}_{t}^{\text{lqmm}}(\tau)$ changes value radically with
$\tau$, whereas $\beta_{it}^{\text{lqmm}}(\tau)$ is a smooth function of
$\tau$ for each species $i$, this implies that $\hat{u}^{\text{lqmm}}_{i}$ is
negatively correlated with $\hat{\beta}^{\text{lqmm}}$, indicating that the
fixed and random effects are difficult to separate.

If $\beta_{it}(\tau)$ is more or less constant in $\tau$, there is not much
need for a quantile regression model. If the location change is the same
across all quantiles, an ordinary linear mixed effects model (not using
empirical quantiles) would be much more suitable. However, inference based on
the eq- and the meq-models indicate that the effect of year varies with $\tau$
for several species.

Figure 6 shows that for all three models and the majority of species and
$\tau$-values, $\hat{\beta}_{it}(\tau)$ are negative. This means that overall,
species have advanced their spring arrival. This is consistent with many
previously published studies of bird phenology (Jonzén et al., 2006; Tøttrup
et al., 2006; Lehikoinen et al., 2019; Usui et al., 2017; Aharon-Rotman et
al., 2021), here with added details of how this effect changes across the
arrival time distribution, both within species and generally across species.
As a closing observations we notice that for most species $i$,
$\hat{\beta}_{it}^{\text{meq}}(\tau)$ and $\hat{\beta}_{it}^{\text{qr}}(\tau)$
follow each other closely.

Figure 6: Estimates of $\beta_{it}(\tau)$ for each species of the multiple
species model. For a few quantiles and species, $\hat{\beta}_{it}(\tau)$ is
located outside the range of the plot. This vertical truncation makes it
easier to inspect where $\hat{\beta}_{it}(\tau)$ is not constant in $\tau$.
Notice the overall positive correlation between the qr and meq estimates,
wheras the lqmm estimates seem more independent of the other methods, and
overall are very flat.

####  4.2.3 Conclusion

The above results indicate that the lqmm framework is unsuitable for a
multiple species analysis of bird phenology, at least for our dataset. The
overall similarities between the qr and meq-models indicates a common signal
being detected in data, and the results seem much more plausible than the
lqmm-results. Considering the overall purpose of the phenological analysis at
hand, we therefore propose to incorporate species as a block factor of fixed
effect.

##  5 Discussion

In this article we fitted empirical quantile models and quantile regression
models to phenological data of bird migration, in the context of single and
multispecies models. Overall the ordinary quantile regression model, without
random effects, provided the most stable and reliable results.

For the multispecies analysis based on a quantile regression model with random
effects (lqmm), the oscillating behavior of the fixed effects estimates
$\hat{\beta}^{\text{lqmm}}(\tau)$, as a function of the quantile $\tau$, is
not ideal (as seen in Figures 4 and 5). Further analyses are needed in order
to find out whether this stems from an objective function with many local
optima (so that the optimization methods could be trapped in one of these), or
if the global maximum varies rapidly with $\tau$. Although it is known that
the log likelihood objective function of the lqmm is concave (Geraci and
Bottai, 2014), the marginalization of the random effects $u$ is done
numerically, using Gaussian quadrature. We cannot rule out that the resulting
approximation of the objective function, unless it is very accurate, has local
optima. The number of knots used in the Gaussian quadrature can be tuned, and
we found that the oscillating behavior of the fixed effect estimates were less
pronounced when we set it to 25.

It would be natural to complement the results of this paper by a simulation
study. A good starting point for such a study could be the Negative Binomial
approach of Lindén (2011, p51), which can be employed repeatedly to simulate a
multiple species data set. Another approach to simulating phenological data
for smoothing can be found in Knudsen et al. (2007), with potentially
influential factors collected in Table 2.

We chose to not include confidence intervals for the parameters of the
multiple species models for two reasons. First, the main focus was to
investigate the difference in estimates between the three approaches (qr, meq
and lqmm). In this context a goal was to investigate whether lqmm is suitable
for phenological analysis of multiple species. Our results indicate that this
is not the case. Therefore, and secondly, we did not make the large
computational effort to perform the built in block bootstrap approach of the
lqmm-package.

We filtered the data provided by the Falsterbo Bird Observatory so that
species with a relatively large number of local breeding birds were omitted.
Some birds of the included species might still be local breeding individuals,
rather than migratory individuals, and it is likely that this contamination is
largest at the end points of the sampling window. In particular, the fact that
the estimate $\hat{\beta}_{t}(\tau)$ of the fixed effect of year approches 0
as the quantile $\tau$ approaches 1, might be caused by non-migratory
individuals. Ideally, one would filter data for birds that probably are not
migrating, or at least assign weights to birds depending on how likely it is
that they are migrating. Such a weighting scheme could employ covariates such
as fat score (Rogers, 1991) or data on repeated captures of birds within a
season.

It would also be ideal to include birds that are recaptured between seasons.
These recaptures generally constitute a very small fraction of the number of
first captures, and probably they have a propensity for being local breeders.
In spite of this, birds recaptured between seasons might give deeper insight
into the variation in arrival date within individual birds, as as an instance
of repeated measurements.

At the beginning of the sampling window some species might already have begun
their migration, and at the end of the sampling window some species might
still be arriving. This results in incomplete measurements for the arrival
distribution of these species. It is important to account for such incomplete
data when presenting measures of location, scale, skewness and kurtosis for
the arrival distributions. Such $L$-functionals measures have been used by
Hössjer and Karlsson (2021), for quantile regression models, in the context of
single species phenological analysis. These methods are also applicable to
censored and truncated data.

The equal weighting of birds in this paper aims at creating a complete picture
of the passage of migratory birds. A natural alternative would be to weight
species equally, and thus focus the study on patterns in the behaviors of
various species. In some studies it might be prefable to focus on the general
behavior of species rather than birds. For the multiple species models, such a
shift is easily accomodated by adjusting the species specific weights.

## Acknowledgements

The authors would like to thank the Falsterbo Bird Observatory for providing
the data.

##  Appendix A Efficient matrix inversion for $\hat{u}$

For all fitted models in this paper, $q=0$. Thus, this appendix will cover
random effects estimation for the specific case that each group (species) has
a random intercept and a random slope. Omitting the random slope will make the
computations less involved, and adding further random slopes will complicate
it. Several grouping factors are not supported by lqmm.

Recall that species $i$ has $n_{i}$ observations, and now denote the number of
observed birds of a particular species year $t$ as $n_{t}$, $t=1,\ldots,T$. We
will consider one species at a time, and therefore omit the species index $i$
in the notation, unless it appears in $n_{i}=\sum_{t=1}^{T}n_{t}$. The arrival
time of bird $j$ year $t$ is then

$y_{tj}=x_{tj}\beta(\tau)+u_{1}+u_{2}t+\varepsilon_{tj}(\tau),$ (29)

where $x_{tj}=(1,t,x_{j1},\ldots,x_{jp})$ includes the covariates of bird $j$
from year $t$. The arrival time of all birds of the particular species can be
written in matrix form as

$Y=X\beta(\tau)+Zu+\varepsilon(\tau)$ (30)

where $Y$ and $\varepsilon(\tau)$ are column vectors that contain the arrival
times $y_{tj}$ and the error terms $\varepsilon_{tj}(\tau)$ of all birds,
stacked yearwise on top of each other,

$\displaystyle Z=\begin{pmatrix}&1_{n_{1}}\\\ 1_{n_{i}}&\vdots\\\
&1_{n_{T}}T\end{pmatrix},\quad u_{i}=\begin{pmatrix}u_{i1}\\\
u_{i2}\end{pmatrix},\quad\Psi(\tau)=\begin{pmatrix}\psi_{11}&\psi_{12}\\\
\psi_{21}&\psi_{22}\end{pmatrix}$ (31)

and $1_{n_{\cdot}}$ denotes a column vector with ones of length $n_{\cdot}$.
We introduce

$B=Z\Psi(\tau)Z^{\top}=\begin{pmatrix}B_{11}&\dots&B_{1T}\\\
\vdots&\ddots&\vdots\\\ B_{T1}&\dots&B_{TT}\end{pmatrix}$ (32)

where

$\displaystyle B_{st}$
$\displaystyle=\left[\psi_{11}+\psi_{12}\left(s+t\right)+\psi_{22}st\right]\left(1_{n_{s}}\otimes
1_{n_{t}}\right)$ (33) $\displaystyle:=b_{st}\left(1_{n_{s}}\otimes
1_{n_{t}}\right).$ (34)

This means that each block $B_{st}$ in $B$ is a matrix of dimension
$n_{s}\times n_{t}$ with value $b_{st}$ everywhere. The repeated occurence of
the same values within blocks is the structure we will make use of to make the
matrix inversion more efficient, by finding the inverse to a reduced version
of $\Sigma$ and then rebuilding it using the block structure. We first
postulate that $\Sigma^{-1}$ should have the same type of block structure as
$\Sigma$, i.e.

$\Sigma^{-1}=C+\delta_{\varepsilon}I_{n_{i}}.$ (35)

Here, $C$ is a matrix with the same structure as $B$, but with $b_{st}$
replaced with $c_{st}$, and $\delta_{\varepsilon}>0$ a scalar. If we first
find $\delta_{\varepsilon}$ and all $c_{st}$, we can then make use of (35) and
compute $\Sigma^{-1}$ from these values and
$\left\\{n_{t}\right\\}_{t=1}^{T}$. From the equality

$\left(B+\sigma_{\varepsilon}I_{n_{i}}\right)\left(C+\delta_{\varepsilon}I_{n_{i}}\right)=BC+\delta_{\varepsilon}B+\sigma_{\varepsilon}C+\sigma_{\varepsilon}\delta_{\varepsilon}I_{n_{i}}=I_{n_{i}}$
(36)

we introduce

$BC:=A=\begin{pmatrix}A_{11}&\dots&A_{1T}\\\ \vdots&\ddots&\vdots\\\
A_{T1}&\dots&A_{TT}\end{pmatrix}$ (37)

and consequently, $A_{st}=a_{st}\left(1_{n_{s}}\otimes 1_{n_{t}}\right)$,
where

$a_{st}=\sum_{r=1}^{T}n_{r}b_{sr}c_{rt}$ (38)

which can be seen by performing the matrix multiplication in (37). Our goal is
to find a matrix of reduced dimension to invert, so to this end we introduce
the matrices

$\check{A}=\begin{pmatrix}a_{11}&\dots&a_{1T}\\\ \vdots&\ddots&\vdots\\\
a_{T1}&\dots&a_{TT}\end{pmatrix}\quad\check{B}=\begin{pmatrix}b_{11}&\dots&b_{1T}\\\
\vdots&\ddots&\vdots\\\
b_{T1}&\dots&b_{TT}\end{pmatrix}\quad\check{C}=\begin{pmatrix}c_{11}&\dots&c_{1T}\\\
\vdots&\ddots&\vdots\\\ c_{T1}&\dots&c_{TT}\end{pmatrix}$ (39)

and the diagonal matrix $\check{N}$ with diagonal elements contained in the
vector $\text{tr}(\check{N})=\left(n_{1},\ldots,n_{T}\right)$. In conjunction
with (38) this gives the equality

$\check{A}=\check{B}\check{N}\check{C}.$ (40)

Applying the dimension reduction ${\cal M}^{n_{i}\times n_{i}}\mapsto{\cal
M}^{T\times T}$ to the RHS of (36) we get

$\check{B}\check{N}\check{C}+\delta_{\varepsilon}\check{B}+\sigma_{\varepsilon}\check{C}+\sigma_{\varepsilon}\delta_{\varepsilon}I_{T}=I_{T}.$
(41)

A solution to (41) is

$\begin{cases}\sigma_{\varepsilon}\delta_{\varepsilon}=1,\\\
\check{B}\check{N}\check{C}+\delta_{\varepsilon}\check{B}+\sigma_{\varepsilon}\check{C}=0_{T}\otimes
0_{T},\end{cases}$ (42)

where the upper equation of (42) gives
$\delta_{\varepsilon}=\sigma_{\varepsilon}^{-1}$, whereas the lower equation
of (42) yields

$\displaystyle\check{C}$
$\displaystyle=-\delta_{\varepsilon}\left(\check{B}\check{N}+\sigma_{\varepsilon}\right)^{-1}\check{B}.$
(43)

With $\delta_{\varepsilon}$ and $\check{C}$ determined, we plug them into (35)
to find $\Sigma^{-1}$, which is then inserted into (26).

However, if $n_{i}$ is large, the memory requirements of (26) might still lead
to a crash on some computers, and thus we adapted (26) to compute
$\hat{u}^{(\tau)}$ sequentially. We divide $\Sigma^{-1}$ into $T$ submatrices
$\left[\Sigma^{-1}_{1}|\ldots|\Sigma^{-1}_{T}\right]$ where each submatrix
consists of adjacent columns of $\Sigma^{-1}$ of dimension $n_{i}\times n_{t}$
for $t=1,\ldots,T$. We then compute for each $t$ the $2\times n_{t}$ matrices

$W_{t}=\hat{\Psi}^{(\tau)}Z^{\top}\hat{\Sigma}^{-1}_{t}$ (44)

and combine them into the $2\times n_{i}$ matrix
$W=\left[W_{1}|\ldots|W_{T}\right]$. As a final step we then compute

$\hat{u}^{(\tau)}=W\left\\{y-X\hat{\beta}^{(\tau)}-\hat{\mathbb{E}}\left[\varepsilon^{(\tau)}\right]\right\\}.$
(45)

The whole procedure is still quite memory intensive, but we manage to run it
on a computer with 32GB of RAM with the most numerous species having 38 738
observations. For species with a few thousand observations, the sequential
procedure is not needed for computing (26), but the dimension reduced approach
to inverting $\Sigma$ is recommended, as it speeds up computation
considerably.

The only added requirement on data for our particular implementation (cf. the
supplementary material), is that birds are ordered by species, and ordered
chronologically within species.

## References

  * Aharon-Rotman et al. (2021) Aharon-Rotman, Y., G. Perlman, Y. Kiat, T. Raz, A. Balaban, and T. Iwamura (2021).  Limited flexibility in departure timing of migratory passerines at the east-mediterranean flyway.  Scientific reports 11(1), 1–13.
  * Bates (2014) Bates, D. (2014).  Computational methods for mixed models.
  * Bates et al. (2015) Bates, D., M. Mächler, B. Bolker, and S. Walker (2015).  Fitting linear mixed-effects models using lme4.  Journal of Statistical Software 67(1), 1–48.
  * Bottai et al. (2015) Bottai, M., N. Orsini, and M. Geraci (2015).  A gradient search maximization algorithm for the asymmetric laplace likelihood.  Journal of Statistical Computation and Simulation 85(10), 1919–1925.
  * Geraci et al. (2014) Geraci, M. et al. (2014).  Linear quantile mixed models: the lqmm package for laplace quantile regression.  Journal of Statistical Software 57(13), 1–29.
  * Geraci and Bottai (2007) Geraci, M. and M. Bottai (2007).  Quantile regression for longitudinal data using the asymmetric laplace distribution.  Biostatistics 8(1), 140–154.
  * Geraci and Bottai (2014) Geraci, M. and M. Bottai (2014).  Linear quantile mixed models.  Statistics and computing 24(3), 461–479.
  * Hinkley and Revankar (1977) Hinkley, D. V. and N. S. Revankar (1977).  Estimation of the pareto law from underreported data: a further analysis.  Journal of Econometrics 5(1), 1–11.
  * Hössjer and Karlsson (2021) Hössjer, O. and M. Karlsson (2021).  On the use of $l$-functionals in regression models.
  * Jonzén et al. (2006) Jonzén, N., A. Lindén, T. Ergon, E. Knudsen, J. O. Vik, D. Rubolini, D. Piacentini, C. Brinch, F. Spina, L. Karlsson, et al. (2006).  Rapid advance of spring arrival dates in long-distance migratory birds.  Science 312(5782), 1959–1961.
  * Karlsson (2014) Karlsson, M. (2014).  Modelling daily numbers of ringed birds with negative binomial generalized linear models.  Stockholm University bachelor thesis in mathematical statistics, 2014:10.
  * Knudsen et al. (2007) Knudsen, E., A. Lindén, T. Ergon, N. Jonzén, J. O. Vik, J. Knape, J. E. Røer, and N. C. Stenseth (2007).  Characterizing bird migration phenology using data from standardized monitoring at bird observatories.  Climate Research 35(1-2), 59–77.
  * Koenker (2005) Koenker (2005).  Quantile Regression (Econometric Society monographs; no. 38).  Cambridge university press.
  * Koenker (2021) Koenker, R. (2021).  quantreg: Quantile Regression.  R package version 5.86.
  * Koenker and Machado (1999) Koenker, R. and J. A. Machado (1999).  Goodness of fit and related inference processes for quantile regression.  Journal of the american statistical association 94(448), 1296–1310.
  * Laird and Ware (1982) Laird, N. M. and J. H. Ware (1982).  Random-effects models for longitudinal data.  Biometrics, 963–974.
  * Lehikoinen et al. (2019) Lehikoinen, A., A. Lindén, M. Karlsson, A. Andersson, T. L. Crewe, E. H. Dunn, G. Gregory, L. Karlsson, V. Kristiansen, S. Mackenzie, et al. (2019).  Phenology of the avian spring migratory passage in europe and north america: Asymmetric advancement in time and increase in duration.  Ecological indicators 101, 985–991.
  * Lindén (2011) Lindén, A. (2011).  Using first arrival dates to infer bird migration phenology.
  * Mayer (2010) Mayer, A. (2010).  Phenology and citizen science: volunteers have documented seasonal events for more than a century, and scientific studies are benefiting from the data.  BioScience 60(3), 172–175.
  * Nelder and Mead (1965) Nelder, J. A. and R. Mead (1965, 01).  A Simplex Method for Function Minimization.  The Computer Journal 7(4), 308–313.
  * Oja (1981) Oja, H. (1981).  On location, scale, skewness and kurtosis of univariate distributions.  Scandinavian Journal of statistics, 154–168.
  * Portnoy and Koenker (1997) Portnoy, S. and R. Koenker (1997).  The gaussian hare and the laplacian tortoise: computability of squared-error versus absolute-error estimators.  Statistical Science 12(4), 279–300.
  * R Core Team (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Rogers (1991) Rogers, C. M. (1991).  An evaluation of the method of estimating body fat in birds by quantifying visible subcutaneous fat (una evaluación del método para estimar la cantidad de grasa en aves, mediante la cuantificación de la grasa subcutanea visible).  Journal of field ornithology, 349–356.
  * Roos and Karlsson (1981) Roos, G. and L. Karlsson (1981).  Ringmarkningsverksamheten vid falsterbo fågelstation 1980.  Anser 20, 99–108.
  * Tøttrup et al. (2006) Tøttrup, A. P., K. Thorup, and C. Rahbek (2006).  Patterns of change in timing of spring migration in north european songbird populations.  Journal of Avian Biology 37(1), 84–92.
  * Usui et al. (2017) Usui, T., S. H. Butchart, and A. B. Phillimore (2017).  Temporal shifts and temperature sensitivity of avian spring migratory phenology: A phylogenetic meta-analysis.  Journal of Animal Ecology 86(2), 250–261.
  * Walther et al. (2002) Walther, G.-R., E. Post, P. Convey, A. Menzel, C. Parmesan, T. J. Beebee, J.-M. Fromentin, O. Hoegh-Guldberg, and F. Bairlein (2002).  Ecological responses to recent climate change.  Nature 416(6879), 389–395.
  * Zehnder and Karlsson (2001) Zehnder, S. and L. Karlsson (2001).  Do ringing numbers reflect true migratory activity of nocturnal migrants?  Journal für Ornithologie 142(2), 173–183.