Parametric Modeling of Quantile Regression Coefficient Functions with Longitudinal Data

PARAMETRIC MODELING OF QUANTILE REGRESSION COEFFICIENT FUNCTIONS WITH LONGITUDINAL DATA

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if11 { \\ Paolo Frumento ([email removed]) is Associate Professor at the Department of Political Sciences, University of Pisa, Italy.\\ Matteo Bottai ([email removed]) is Professor at the Unit of Biostatistics at Karolinska Institutet, Institute of Environmental Medicine, Stockholm, Sweden.\\ Iv\'an Fern\'andez-Val ([email removed]) is Professor at the Department of Economics, Boston University.\\ Address for correspondence: Department of Political Sciences, University of Pisa, Via F. Serafini, 3, 56126 Pisa, Italy. } \fi

\if01 {

} \fi

{\it Keywords:} Longitudinal quantile regression, two-level quantile function, parametric quantile function, penalized fixed-effects, R package qrcm.

\spacingset{1.45}

INTRODUCTION

Quantile regression (e.g., kb, koenker) has become a standard method in many fields, including medicine, epidemiology, economics, and social sciences. Different solutions have been proposed to extend quantile regression to longitudinal data, in which the same individuals or clusters are observed repeatedly.

In conditional models, that include fixed- and random-effects models, the dependence between observations is accounted for by introducing individual-specific parameters, or “individual effects”. In fixed-effects models, the individual effects are treated as parameters, avoiding distributional assumptions and allowing for a simple computation. A penalized fixed-effects estimator for longitudinal quantile regression has been proposed by koenker2004, and similar approaches have been used in lamarche, canay, and kato.\footnote{cfw18 considered an alternative to quantile regression for estimation of quantile effects in longitudinal data based on distribution regression.} In random-effects models, the individual effects are described by a parametric distribution. Different methods have been proposed to combine the parametric likelihood of the random effects with the estimating equation of ordinary quantile regression. geraci (geraci, geraci2) used the log-likelihood of an asymmetric Laplace distribution, and kim described an empirical likelihood method. ad08 adapted the correlated random effect approach of c84 to quantile regression, and arellano marginalized the loss function of quantile regression with respect to the posterior distribution of the individual effects. farcomeni, marino, and alfo used finite mixtures to approximate the probability density function of the individual effects through a discrete distribution.

Marginal models have also been described in the literature. leng defined a set of unbiased estimating equations carrying information on the correlation structure. A similar approach was used by zhao to implement longitudinal single-index quantile regression.

In this paper, we adopt the conditional paradigm and introduce a two-level quantile function, in which both the distribution of the within-subject response (level 1) and that of the individual effects (level 2) are described by quantile regression models. With this approach, the distribution of the individual effects is not subject to strong parametric assumptions and is allowed to depend on level-2 covariates. Following iqr (2016, 2017) and yang, we describe the level-1 and level-2 quantile regression coefficients by (flexible) parametric functions of the order of the quantile. Compared with standard quantile regression, in which quantiles are estimated one at a time, this modeling approach presents numerous advantages, that include a simpler computation and inference (owing to a smooth objective function), increased statistical efficiency, and easy interpretability of the results.

To fit the model, we introduce a new form of penalized fixed-effects estimator in which the penalty term carries information on level-2 parameters. This method presents important advantages over standard $\ell_1$ and $\ell_2$ penalization. In particular, it avoids the problem of selecting a tuning constant, and allows to estimate the level-2 coefficients of the model using fixed-effects techniques.

The paper is structured as follows. We describe a general model in Section (ref), and discuss model building in Section (ref). We introduce an estimator in Section (ref), and in Section (ref) we derive its asymptotic properties. In Section (ref) we present goodness-of-fit measures and tools for model selection, and in Section (ref) we report simulation results. Section (ref) concludes the paper with the analysis of a dataset relating plasma neutrophil gelatinase-associated lipocalin (NGAL) to sepsis status. Appendix A provides a general asymptotic expansion for fixed effects estimators with mixed-rates asymptotics and applies it to derive the asymptotic distribution of the proposed estimator. We discuss computation in Appendix B, and present extended simulation results in Appendix C. The R package qrcm implements the described estimator and provides a variety of auxiliary functions for model building, summary, plotting, prediction, and goodness-of-fit assessment.

THE MODEL

A two-level quantile function

We consider a cluster data structure, in which $N$ individuals or clusters are observed repeatedly. We denote by $i = 1, \ldots, N$ the index of the subject, and by $t = 1, \ldots, T$ the within-subject index, such that the total sample size is $NT$. Designs in which $T$ varies across clusters are also possible, at the cost of a slightly more complicated notation.

We denote by $Y_{it}$ a response variable of interest, and assume that

where $\bm{x}_{it}$ is a $d_x$-dimensional vector of level-1 covariates, with associated parameter $\bm{\beta}(\cdot)$; and $\bm{z}_i$ is a $d_z$-dimensional vector of level-2 covariates, with associated parameter $\bm{\gamma}(\cdot)$.

We assume that (i) $\bm{x}_{it}^{ \mathrm{\scriptscriptstyle T} }\bm{\beta}(\cdot)$ and $\bm{z}_{i}^{ \mathrm{\scriptscriptstyle T} }\bm{\gamma}(\cdot)$ are a.s. non-decreasing functions of their arguments, and (ii) $U_{it}$ and $V_i$ are $U(0,1)$ variables, independent of each other and of the covariates. Based on model ((ref)), $\alpha_i = \bm{z}_{i}^{ \mathrm{\scriptscriptstyle T} }\bm{\gamma}(V_{i})$ is an individual effect with conditional quantile function $\bm{z}_{i}^{ \mathrm{\scriptscriptstyle T} }\bm{\gamma}(\cdot)$, while $\bm{x}_{it}^{ \mathrm{\scriptscriptstyle T} }\bm{\beta}(\cdot)$ is the conditional quantile function of $Y_{it} - \alpha_i$.

The level-1 quantile regression model, $\bm{x}_{it}^{ \mathrm{\scriptscriptstyle T} }\bm{\beta}(\cdot)$, has the standard interpretation (e.g., koenker2004): it characterizes the “within” part of the distribution, purged of the individual effects. The level-2 regression model, $\bm{z}_{i}^{ \mathrm{\scriptscriptstyle T} }\bm{\gamma}(\cdot)$, describes the distribution of the between-subject differences with respect to a reference value which typically corresponds to a “mean” or “median” individual.

Consider, for example, a clinical study in which patients are repeatedly measured their body mass index (BMI) during their lifetime. The level-1 part of the model describes the conditional quantiles of BMI in a “typical” patient, i.e., someone with an individual effect equal to $0$. Level-1 predictors include time-varying characteristics, such as the age of the patient at each observation, as well as constant traits, such as the gender of the patient. The level-2 model accounts for the between-patient heterogeneity, and describes the conditional quantiles of the individual effects. Level-2 covariates can only include time-invariant traits, such as the gender, and summary statistics of level-1 covariates, e.g., the age at the first examination. Note that the dimension of the level-1 covariates, $\bm{x}_{it}$, is $NT$, while that of the level-2 covariates, $\bm{z}_i$, is $N$.

Unlike the “standard” approaches, that do not consider the effect of level-2 covariates, our modeling framework allows to investigate the determinants of the between-subject variability. For example, in the linear random-intercept model, the level-2 response is described by a $N(0, \phi^2)$ distribution, in which $\phi^2 = \text{var}(\alpha_i)$ is interpreted as the “between” variance and is assumed to be unaffected by predictors. This model may fail to capture important features of the data, such as the fact that the variance of the individual effects is different in males and females. Model (ref), instead, allows including gender as level-2 predictor.

Using a quantile regression approach permits avoiding strong parametric assumptions such as normality and homoskedasticity, that are often used in likelihood-based modeling. In the existing literature on longitudinal quantile regression, however, a quantile regression model is usually applied to the level-1 response, but not to the individual effects, that are treated as nuisance parameters. In our paradigm, instead, the two parts of the distribution are considered “equally important”, in the sense that the same modeling structure is used to describe the quantiles of the within-subject response, and those of the between-subject differences. As shown later in the paper, working with model (ref) permits using the same techniques to estimate both the level-1 and the level-2 parameters, and avoids combining level-1 quantile regression methods with likelihood-based level-2 estimators as for example in kim. This leads to rather simple procedures for estimation and inference, in which a fundamental role is played by the two independent uniform random variables ($U_{it}, V_i$) that generate the data.

Parametric coefficient functions

Through the paper, we assume that the quantile regression coefficient functions, $\bm{\beta}(\cdot)$ and $\bm{\gamma}(\cdot)$, can be modeled parametrically:

where $\bm{\theta}$ and $\bm{\phi}$ are unknown model parameters. This modeling approach was used by iqr (2016, 2017), and is exemplified in Figure 1. The broken line in figure represents standard regression coefficients at quantiles $u = (0.01,0.02, \ldots, 0.99)$. The estimated coefficients show a non-smooth, volatile trend and, although consistently positive, are almost never significant. A parametric model can be used to characterize the coefficient function with few parameters and describe it by a simple, closed-form mathematical expression. In Figure 1 we propose a linear fit, $\beta(u \mid \bm{\theta}) = \theta_0 + \theta_1 u$, that is represented by a dashed line. This simple model reveals the underlying trend and permits achieving statistical significance.

Compared with standard quantile regression, which works in a quantile-by-quantile fashion, modeling quantile functions parametrically simplifies estimation and inference and yields important advantages in terms of parsimony, efficiency, and ease of interpretation. Moreover, it allows for model identification in the presence of latent structures or missing information, making it simple to apply quantile regression to censored and truncated data ctiqr.

On the other hand, this approach requires formulating a parametric model for the coefficient functions, $\bm{\beta}(u\mid\bm{\theta})$ and $\bm{\gamma}(v\mid\bm{\phi})$. This task is not straightforward and the existing literature on the subject is lacking. In Section (ref) we describe in details model building, provide guidelines, and suggest a variety of possible parametrizations.

TWO-LEVEL MODELING OF QUANTILE REGRESSION COEFFICIENT FUNCTIONS

We assume model ((ref)) to hold, and parametrize the quantile regression coefficient functions as follows:

where $\bm{b}(u) = \left[b_1(u), \ldots, b_{d_{b}}(u)\right]^{{ \mathrm{\scriptscriptstyle T} }}$ and $\bm{c}(v) = \left[c_1(v), \ldots, c_{d_{c}}(v)\right]^{{ \mathrm{\scriptscriptstyle T} }}$ are $d_{b}$- and $d_{c}$-dimensional sets of known functions. With this notation, $\bm{\theta}$ is a $d_{x}\times d_{b}$ matrix, and $\bm{\phi}$ is a $d_{z}\times d_{c}$ matrix. The data-generating process can be written as

Although other parametrizations are possible (e.g., $\bm{\beta}(u\mid\bm{\theta})$ and $\bm{\gamma}(v\mid\bm{\phi})$ may be allowed to be nonlinear functions of $\bm{\theta}$ and $\bm{\phi}$), model (ref) is very flexible and computationally convenient. We illustrate the potentials of this modeling approach with a number of examples, and provide general guidelines for model building.

A simple model

Consider the following model with a single level-1 covariate $x$, and no level-2 predictors: $$Y_{it} = \beta_0(U_{it}) + \beta_1(U_{it})x_{it} + \gamma_0(V_i).$$ Denote by $\zeta(\cdot)$ the quantile function of a standard normal distribution, and assume that $$\beta_0(u\mid\bm{\theta}) = \theta_{00} + \theta_{01}\zeta(u),$$ $$\beta_1(u\mid\bm{\theta}) = \theta_{10},$$ $$\gamma(v \mid\bm{\phi}) = \phi \zeta(v).$$ This is just a reformulation of the standard linear random-intercept model, in which $Y_{it} = \theta_{00} + \theta_{10}x_{it} + \alpha_i + \epsilon_{it}$ with $\alpha_i \sim N(0,\phi^2)$ and $\epsilon_{it} \sim N(0,\theta_{01}^2)$. In this model, $\theta_{00}$ corresponds to the intercept of the “fixed” part, while $\theta_{01}^2$ and $\phi^2$ are interpreted as the “within” and “between” variance components. In the equivalent quantile regression model, $\theta_{00}$ is the “intercept” of $\beta_0(u\mid\bm{\theta})$ and corresponds to $\beta_0(0.5\mid\bm{\theta})$, while $\theta_{01}$ and $\phi$ are “slopes” associated with $\zeta(\cdot)$ in the level-1 and level-2 part of the quantile function, respectively. The regression coefficient of $x$, $\beta_1(u\mid\bm{\theta})$, is assumed to be constant across quantiles, forcing homoskedasticity.

A more flexible model

The standard linear random-intercept model is rather restrictive and, within the described framework, can be easily generalized by choosing a different specification of $\bm{\beta}(\cdot \mid \bm{\theta})$ and $\bm{\gamma}(\cdot \mid \bm{\phi})$. For example, one may define $$\beta_0(u\mid\bm{\theta}) = \theta_{00} + \theta_{01}u + \theta_{02}u^2 + \theta_{03}u^3 + \theta_{04}\zeta(u),$$ $$\beta_1(u\mid\bm{\theta}) = \theta_{10} + \theta_{11}u,$$ $$\gamma(v \mid\bm{\phi}) = \phi_1 \log{(2v)} + \phi_2\log{(2(1 - v))}.$$ The intercept, $\beta_0(u\mid\bm{\theta})$, is modeled by a linear combination of $\zeta(u)$, the quantile function of a standard normal distribution, and three additional components, $u$, $u^2$ and $u^3$, that allow for a deviation from the normal model. The resulting quantile function can be asymmetric or multimodal and does not correspond to any “standard” family of random variables. The coefficient associated with $x$, $\beta_1(u\mid\bm{\theta})$, is now assumed to be a linear function of $u$, allowing for data heteroskedasticity. In particular, the variance of the level-1 response is an increasing function of $x$, if $\theta_{11} > 0$, and a decreasing function of it, if $\theta_{11} < 0$. Finally, the individual effects are assumed to follow a zero-median asymmetric logistic distribution, which is much more flexible than the commonly used normal model.

As shown in this example, $\bm{\beta}(\cdot)$ and $\bm{\gamma}(\cdot)$ can be constructed as linear combinations of relatively simple functions, $\bm{b}(\cdot)$ and $\bm{c}(\cdot)$, such that $\bm{\beta}(u\mid\bm{\theta}) = \bm{\theta}\bm{b}(u)$, and $\bm{\gamma}(v\mid\bm{\phi}) = \bm{\phi}\bm{c}(v)$. In this framework, the model is entirely determined by the choice of $\bm{b}(\cdot)$ and $\bm{c}(\cdot)$. Useful guidelines for model building are provided in the rest of this section. Various modeling approaches are illustrated in Sections (ref) and (ref) of this paper, while a general discussion on quantile modeling can be found in the book by gil. Finally, the documentation of the qrcm package (in particular the functions iqr and iqrL) includes an extensive tutorial for the practitioners.

Model building: level 1

Modeling $\beta_0(u \mid \bm{\theta})$. Assuming that the support of $\bm{x}$ includes the zero (which can be obtained by centering the covariates), $\beta_0(\cdot \mid \bm{\theta})$ must be a monotonically increasing function. Prior belief or knowledge can be used to identify a meaningful parametric model. For instance, one may use the quantile function of a known distribution. Possible parametrizations of $\beta_0(u \mid \bm{\theta})$ include: $\theta_{00} + \theta_{01}\zeta(u)$, the normal distribution, $\text{N}(\theta_{00}, \theta_{01}^2)$; $-\theta_{01}\log(1 - u)$, the exponential distribution, $\text{Exp}(\theta_{01})$; $\theta_{00} + \theta_{01}\log(u/(1 - u))$, the logistic distribution, $\text{Logis}(\theta_{00}, \theta_{01})$; $\theta_{00} + \theta_{01}\log(u) + \theta_{02}\log(1 - u)$, the asymmetric logistic, $\text{aLogis}(\theta_{00}, \theta_{01}, \theta_{02})$; $\theta_{00} + \theta_{01}u$, the uniform distribution, $\text{U}(\theta_{00}, \theta_{00} + \theta_{01})$. Note that, in this framework, the parameters of well-known distributions may have an unusual interpretation. For example, the value of $\theta_{01}$ in a $\text{U}(\theta_{00}, \theta_{00} + \theta_{01})$ distribution corresponds to its range, but can also be seen as the slope of a linear quantile function, $\theta_{00} + \theta_{01}u$.\\ Modeling $\beta_1(u \mid \bm{\theta}), \beta_2(u \mid \bm{\theta}), \ldots$. There are no general constraints to the parametric form of the regression coefficients associated with the covariates. However, the coefficient functions are usually bounded and exhibit a rather simple behavior. Sometimes, it is possible to assume that covariates only affect the location of the level-1 response, and force homoskedasticity by choosing a constant-slope model in which $\beta_j(u \mid \bm{\theta}) = \theta_{j0}, j = 1, 2, \ldots$. In a more general scenario, a useful approximation is often given by a linear-slope model, $\beta_j(u \mid \bm{\theta}) = \theta_{j0} + \theta_{j1}u$, or a quadratic-slope model, $\beta_j(u \mid \bm{\theta}) = \theta_{j0} + \theta_{j1}u + \theta_{j2}u^2$, which does not impose monotone effect with respect to $u$.\\

Model building: level 2

A similar model strategy can be applied to the level-2 quantile function. There are, however, some important differences.\\ Modeling $\gamma_0(v \mid \bm{\phi})$. The distribution of the individual effects is typically assumed to have zero mean or median, and, for identifiability, $\gamma_0(v \mid \bm{\phi})$ does not usually include a constant term. Meaningful definitions of $\gamma_0(v \mid \bm{\phi})$ include: $\phi_{01}\zeta(v)$, the normal distribution, $\text{N}(0, \phi_{01}^2)$; $-\phi_{01}\log(1 - v)$, the exponential distribution, $\text{Exp}(\phi_{01})$; $\phi_{01}\log(v/(1 - v))$, the logistic distribution, $\text{Logis}(0, \phi_{01})$; $\phi_{01}\log(2v) + \phi_{02}\log(2(1 - v))$, a zero-median asymmetric logistic; $\phi_{01}[\log(v) + 1] + \phi_{02}[\log(1 - v) + 1]$, a zero-mean asymmetric logistic; $\phi_{01}2(v - 0.5)$, a centered uniform distribution, $\text{U}(-\phi_{01}, \phi_{01})$. In most cases, the coefficients can be interpreted as scale parameters, while the centrality parameter is fixed and equal to zero. In the exponential case, the value $0$ is the minimum of the support of the individual effects, and not a measure of central tendency, while both the mean and the standard deviation of the individual effects correspond to $\phi_{01}$.\\ Modeling $\gamma_1(v \mid \bm{\phi}), \gamma_2(v \mid \bm{\phi}), \ldots$. Importantly, the described framework permits investigating how the conditional quantile function of the individual effects depends on level-2 covariates $\bm{z}_i$, which typically include cluster-invariant characteristics (e.g., gender) or summary measures of the level-1 covariates, e.g., the cluster means or medians. The variance of the individual effects is likely to differ across subgroups of the population. Also, as suggested by some authors (e.g., lancaster), agents may select their covariates' values based on prior knowledge about their own individual effect, which induces a correlation between $\alpha_i$ and $\bm{z}_i$.

Modeling the effect of level-2 covariates is not trivial. To make an example, suppose that $\alpha_i = \gamma_0(V_i) + \gamma_1(V_i)z_i$, and consider the following alternative parametrizations:

In model (i), where $\gamma_0(v \mid \bm{\phi})$ and $\gamma_1(v \mid \bm{\phi})$ are symmetric around the zero, the conditional distribution of $\alpha_i$ has zero mean and median at all values of $z_i$. The covariate only affects the scale of the individual effects by introducing heteroskedasticity, while no linear correlation between $z_i$ and $\alpha_i$ is present. Model (i) assumes normality, but allows the variance of the individual effects to be a function of the level-2 covariates, i.e. $\alpha_i \mid z_i \sim \text{N}(0, \phi_{01}^2 + \phi_{11}^2 z_i^2)$. For example, if $z_i$ is binary, the “between” variance is $\phi_{01}^2$ when $z_i = 0$, and $\phi_{01}^2 + \phi_{11}^2$ when $z_i = 1$.

In models (ii) and (iii), $z_i$ and $\alpha_i$ have a non-zero correlation unless $\phi_{10} = 0$. In model (ii), where $\int_0^1\gamma_0(v \mid \bm{\phi}) \mathrm{d}v = 0$, the marginal distribution of the individual effects has zero mean if $z_i$ is centered around its mean or $\phi_{10} + \phi_{11}/2 = 0$. In model (iii) the mean and the median of the individual effects are functions of the parameters and cannot be determined in advance. However, if $z_i \ge 0$, model (iii) generates $\alpha_i \ge 0$ for any positive value of the parameters, implying that the “reference” individual ($\alpha_i = 0$) corresponds to someone with the smallest possible individual effect.

Additional remarks

The problem of formulating a parametric quantile function is equivalent, at least in principle, to that of choosing a parametric form for a probability density function, a hazard function, or a survival function. For example, as shown in Section (ref), standard parametric assumptions such as normality and homoskedasticity can be directly translated into a quantile function with a simple closed-form expression. However, as suggested in Section (ref), the models that can be used to describe a quantile function are often very different from most of the “conventional” parametric distributions, and frequently much more flexible.

An exploratory semiparametric fit can be obtained by letting $\bm{b}(\cdot)$ and $\bm{c}(\cdot)$ be the basis of a linear or polynomial spline. A flexible model can be used as a guide to find more parsimonious and efficient parametrizations. Note that standard quantile regression, in which quantiles are estimated one at a time, can be thought of as a model in which $\bm{b}(\cdot)$ and $\bm{c}(\cdot)$ are allowed to be arbitrarily flexible and the parameters $\bm{\theta}$ and $\bm{\phi}$ are virtually infinite-dimensional.

In absence of prior knowledge, one may define $\bm{b}(\cdot)$ and $\bm{c}(\cdot)$ using polynomials $\left[\text{e.g.,} u, u^2, u^3, \ldots\right]$, roots $\big[\text{e.g.,} u^{1/2}$, $(1 - u)^{1/2}$, $u^{1/3}$, $(1 - u)^{1/3}, \ldots\big]$, trigonometric functions $\left[\text{e.g.,} \cos(2\pi u), \sin(2\pi u)\right]$, splines, and combinations of the above. A possible strategy is to consider a “simple” quantile function (e.g., that of a normal or an exponential distribution, depending on the nature of the outcome) and allow for a departure from it, as suggested in Section (ref).

Importantly, the model specification should reflect assumptions on the shape, support, and boundedness (or unboundedness) of the level-1 and level-2 responses. For example, if the individual effects are believed to be symmetric, $\gamma_0(v \mid \bm{\phi})$ could be formed by the quantile function of a normal or logistic distribution. If the level-1 distribution has a long right tail, $\beta_0(u \mid \bm{\theta})$ may have a positive asymptote in $u = 1$, e.g., $\beta_0(u \mid \bm{\theta}) = \theta_{00} - \theta_{01} \log(1 - u) + \ldots$. On the other hand, if the outcome is strictly positive, building blocks such as $\log(u)$ or $\zeta(u)$, that present a negative asymptote in $u = 0$, may not be appropriate.

Apart from the above important considerations, the choice of $\bm{b}(\cdot)$ and $\bm{c}(\cdot)$ is not as crucial as it appears. For example, the coefficient function defined by $\beta(u) = (u - 0.3)^3$ is almost identical to $\beta(u) = -1.87 + 6.20u + 1.84\cos(u) - 5.92\sin(u)$, the correlation between the two being about $0.99999$. The fact that very different model specifications can be indistinguishable in terms of model fit is unsurprising (for example, it is almost impossible to distinguish a Normal distribution, a Student's t distribution with large degrees of freedoms, and a Gamma distribution with large shape parameter), and suggests that meaningful criteria for model selection should include parsimony and interpretability.

Often, a rather restrictive model may provide a reasonable approximation of the true data distribution, and can be preferred to a more correct model because of its simplicity. Also, parsimonious models are very rewarding in terms of precision, although they may introduce some bias. This explains why strong parametric assumptions, such as homoskedasticity and proportionality of hazards or odds, are used routinely in statistical analysis. In quantile regression, very convenient assumptions are represented by the constant-slope model (e.g., $\beta(u \mid \bm{\theta}) = \theta_0$), in which a certain predictor has the same effect at all quantiles, and the linear-slope model (e.g., $\beta(u \mid \bm{\theta}) = \theta_0 + \theta_1u$), in which a quantile regression coefficient is assumed to be a linear function.

THE ESTIMATOR

iqr considered cross-sectional data $(y_i, \bm{x}_i)$ and defined $\bm{\beta}(u\mid\bm{\theta}) = \bm{\theta}\bm{b}(u)$ as in ((ref)). To estimate $\bm{\theta}$, they suggested minimizing

which is the integral, with respect to the order of the quantile, of the loss function of standard quantile regression, $\rho_u(w) = w(u - I(w \le 0))$ being the “check” function. This estimation method is referred to as integrated loss minimization (ilm) and is currently implemented in the qrcm R package.

To generalize this idea to longitudinal data, assume model ((ref)) holds,

and denote by $y_{it}$ a realization of $Y_{it}$. If the individual effects $\alpha_i = \bm{z}_{i}^{ \mathrm{\scriptscriptstyle T} }\bm{\phi}\bm{c}(V_i)$ were known, one could directly apply the $\textsc{ilm}$ estimator to $y_{it} - \alpha_i$, to compute an estimate of $\bm{\theta}$; and to $\alpha_i$, to compute an estimate of $\bm{\phi}$. This would require solving \[\min_{\bm{\theta}} L_1(\bm{\theta}, \bm{\alpha}_N), \hspace{0.2cm} \min_{\bm{\phi}} L_2(\bm{\phi}, \bm{\alpha}_N)\] where $\bm{\alpha}_N = (\alpha_1, \ldots, \alpha_N)$,\footnote{We index $\bm{\alpha}_N$ by $N$ to emphasize that the dimension grows with the sample size. }

To obtain expressions ((ref))\footnote{ The expression for $L_1(\bm{\theta}, \bm{\alpha}_N)$ bears some similarity to koenker2004's\ (koenker2004) loss function for unpenalized fixed-effects quantile regression, which is defined by $L(\bm{\beta}, \bm{\alpha}_N) = \sum_j \sum_i \sum_t w_j \rho_{u_j}(y_{it} - \alpha_i - \bm{x}_{it}\bm{\beta}(u_j))$ and can be seen as a discretized, non-parametrized, and weighted version of $L_1(\bm{\theta}, \bm{\alpha}_N)$.} and ((ref)), we used equation (9) from iqr, and define

In the formulas, $u_{it}(\bm{\theta}, \alpha_i)$ and $v_i(\bm{\phi}, \alpha_i)$ are such that $y_{it} - \alpha_i = \bm{x}_{it}^{ \mathrm{\scriptscriptstyle T} }\bm{\theta}\bm{b}(u_{it}(\bm{\theta}, \alpha_i))$ and $\alpha_i = \bm{z}_{i}\bm{\phi}\bm{c}(v_i(\bm{\phi}, \alpha_i))$, respectively. This also implies that

is the cumulative distribution of $Y_{it} - \alpha_i$, given $\bm{x}_i$, with parameter $\bm{\theta}$; and

is the cumulative distribution of $\alpha_i$, given $\bm{z}_i$, with parameter $\bm{\phi}$.

In practice, the vector $\bm{\alpha}_N$ of individual effects is not known and must be estimated. We propose estimating $(\bm{\theta}, \bm{\phi}, \bm{\alpha}_N)$ by solving

The proposed loss function is similar to that of a penalized fixed-effects estimator in which $L_2(\bm{\phi}, \bm{\alpha}_N)$ plays the role of a penalty term. Intuitively, $L_2(\bm{\phi}, \bm{\alpha}_N)$ shrinks the estimated fixed effects towards their assumed conditional distribution, introducing some degree of smoothing, improving model identification and efficiency, and avoiding overfitting. At the same time, $L_2(\bm{\phi}, \bm{\alpha}_N)$ carries information on the parameter $\bm{\phi}$ that describes the quantile function of $\bm{\alpha}_N$.

Since both $\bm{\alpha}_N$ and $\bm{\phi}$ are treated as parameters, this approach combines features of fixed-effects estimators, which only estimate $\bm{\theta}$ and $\bm{\alpha}_N$, and random-effects models, which directly estimate $\bm{\theta}$ and $\bm{\phi}$. Computation, however, is much simpler than that of purely random-effects methods (e.g., kim, arellano).

The gradient functions of $L(\bm{\theta}, \bm{\phi}, \bm{\alpha}_N) = L_1(\bm{\theta}, \bm{\alpha}_N) + L_2(\bm{\phi}, \bm{\alpha}_N)$ can be written as

where $\operatorname{vec}$ denotes the vectorization operator and $\otimes$ the kronecker product. The model parameters, $(\bm{\theta}, \bm{\phi}, \bm{\alpha}_N)$, only enter equations (ref)--(ref) through the cumulative distribution functions $u_{it}(\bm{\theta}, \alpha_i)$ and $v_i(\bm{\phi}, \alpha_i)$ defined in (ref) and (ref). Note that $\bm{G}_{\bm{\theta}}(\bm{\theta}, \bm{\alpha}_N)$ does not carry information on $\bm{\phi}$, and $\bm{G}_{\bm{\phi}}(\bm{\phi}, \bm{\alpha}_N)$ does not carry information on $\bm{\theta}$; while $G_{\alpha_i}(\alpha_i, \bm{\theta}, \bm{\phi})$ depends on both $\bm{\theta}$ and $\bm{\phi}$. As shown by iqr, $\bm{G}_{\bm{\theta}}(\bm{\theta}, \bm{\alpha}_N)$ and $\bm{G}_{\bm{\phi}}(\bm{\phi}, \bm{\alpha}_N)$ approach zero when the distributions of $u_{it}(\bm{\theta}, \alpha_i)$ and $v_i(\bm{\phi}, \alpha_i)$ tend to be uniform. This reflects the data-generating process described in ((ref)), which involves the two independent uniform variables $U_{it}$ and $V_i$.

Equation ((ref)) clarifies the role of the “penalty” term $L_2(\bm{\phi}, \bm{\alpha}_N)$:

A desirable property of the proposed penalization is that it only affects the estimates of $\bm{\alpha}_N$ when the clusters are relatively small. As $T \to \infty$, each cluster contains sufficient information to estimate its own individual effect and, consistently, the penalty term $(v_i(\bm{\phi}, \alpha_i) - 0.5)$ in equation ((ref)) becomes irrelevant.

Estimation can be performed by the following iterative process: (i) given $\bm{\alpha}_N$, estimate $\bm{\theta}$ and $\bm{\phi}$ separately by solving $\bm{G}_{\bm{\theta}}(\bm{\theta}, \bm{\alpha}_N)=0$ and $\bm{G}_{\bm{\phi}}(\bm{\phi}, \bm{\alpha}_N)=0$; (ii) given $(\bm{\theta}, \bm{\phi})$, compute a new estimate of $\bm{\alpha}_N$ by solving $G_{\alpha_i}(\alpha_i, \bm{\theta}, \bm{\phi}) = 0$, $i = 1, \ldots, N$. Step (i) can be implemented with standard routines available in the qrcm package, while step (ii) requires finding the zero of $N$ univariate estimating equations. Neither $u_{it}(\bm{\theta}, \alpha_i)$ nor $v_i(\bm{\phi}, \alpha_i)$ are generally available in closed form, and can be evaluated by using a bisection algorithm. Note that the objective function defined by ((ref)) is a smooth function of all parameters, unlike the loss function of standard quantile regression.

The fact that the quantile function may be ill-defined at some value of the parameters can be an issue during estimation. In the implementation of the qrcm package, we use unconstrained optimization from carefully chosen initial values. The algorithm is described in detail in Appendix B.

A new family of penalized fixed-effects estimators

A possible interpretation of the proposed loss function,

is to consider $L_2(\bm{\phi}, \bm{\alpha}_N)$ as a penalty term that shrinks the estimated individual effects towards their conditional median, $\bm{z}_i^{ \mathrm{\scriptscriptstyle T} }\bm{\gamma}(0.5 \mid \bm{\phi})$. Unlike standard penalizations, however, $L_2(\bm{\phi}, \bm{\alpha}_N)$ may depend on level-2 covariates and is a function of estimated parameters.

To clarify this idea, consider a more traditional penalized loss function

where $L_2(\bm{\alpha}_N)$ is a penalty term which does not contain $\bm{\phi}$, and $\lambda$ is a tuning parameter. Common choices of $L_2(\bm{\alpha}_N)$ are the $\ell_1$-penalization, $L_2(\bm{\alpha}_N) = \sum_{i = 1}^N|\alpha_i|$, which was used by koenker2004 to implement longitudinal quantile regression, and the $\ell_2$-penalization, $L_2(\bm{\alpha}_N) = \sum_{i = 1}^N \alpha_i^2$.

Standard $\ell_1$- and $\ell_2$-penalized fixed-effects methods are computationally simple and can substantially improve efficiency of the estimates of the structural parameters. However, besides the fact that they do not allow for estimation of $\bm{\phi}$, they present some important limitations: (i) they do not use prior knowledge on the distribution of the individual effects; (ii) they can introduce bias; (iii) they apply the same penalization to all clusters; and (iv) they require to specify a tuning parameter.

For instance, $\ell_1$-penalized estimators of quantile regression coefficients are asymptotically biased unless $(\alpha_1, \ldots, \alpha_N)$ are independent and identically distributed with zero median lamarche. This is just a consequence of the $\ell_1$-penalty term being a sum of absolute deviations from zero, which does not generally reflect the true distribution of $\bm{\alpha}_N$ and the effect of level-2 covariates on it. Moreover, the same value of $\lambda$ is used for all clusters, ignoring the fact that the variance of the individual effects may differ across subgroups of the population.

The tuning constant $\lambda$ determines the degree of shrinking and, in the standard random-intercept linear model, its optimal value is $\sigma^2_{\epsilon}/\sigma^2_{\alpha}$, i.e., a function of nuisance scale parameters (e.g., koenker2004). Outside the restrictive conditions of linear models, not only the choice of $\lambda$ becomes problematic, but also the use of a single value of $\lambda$ for all clusters is questionable.

The novelty of our approach is that, unlike the $\ell_1$- and $\ell_2$-penalizations, the term $L_2(\bm{\phi}, \bm{\alpha}_N)$ reflects the true (conditional) distribution of $\bm{\alpha}_N$ and carries information about its parameters, $\bm{\phi}$. Our estimator presents the following advantages over standard penalized methods: (i) it enables incorporating parametric assumptions on the distribution of $\bm{\alpha}_N$; (ii) it permits estimating all parameters consistently; (iii) it applies a different degree of shrinking to each cluster, by modeling the effect of level-2 covariates on the distribution of the individual effects; and (iv) it does not require selecting a tuning constant, as no nuisance parameters are present.

To clarify point (iv), consider the loss function of an $\ell_2$-penalized linear regression model: $L_{\lambda}(\bm{\beta}, \bm{\alpha}_N) = L_1(\bm{\beta}, \bm{\alpha}_N) + \lambda L_2(\bm{\alpha}_N) = \sum_{i = 1}^N\sum_{t = 1}^T{(y_{it} - \bm{x}_i^{ \mathrm{\scriptscriptstyle T} }\bm{\beta} - \alpha_i)^2} + \lambda \sum_{i = 1}^N \alpha_i^2$. Here, $L_1(\bm{\beta}, \bm{\alpha}_N)$ and $L_2(\bm{\alpha}_N)$ lack information on the nuisance scale parameters $\sigma^2_{\epsilon} = \text{var}(Y_{it} - \bm{x}_i^{ \mathrm{\scriptscriptstyle T} }\bm{\beta} - \alpha_i)$ and $\sigma^2_{\alpha} = \text{var}(\alpha_i)$. This is adjusted for by the tuning constant $\lambda = \sigma^2_{\epsilon}/\sigma^2_{\alpha}$. In our special type of penalized estimator, instead, $L_1(\bm{\theta},\bm{\alpha}_N)$ and $L_2(\bm{\phi},\bm{\alpha}_N)$ carry information on all model parameters. Intuitively, this means that $L_1(\bm{\theta},\bm{\alpha}_N)$ and $L_2(\bm{\phi},\bm{\alpha}_N)$ are already “properly scaled”. The tuning constant can be thought of as an implicit parameter, a function of $\bm{\theta}$ and $\bm{\phi}$. Although a more general estimator with criterion function $L_1(\bm{\theta},\bm{\alpha}_N) + \lambda L_2(\bm{\phi},\bm{\alpha}_N)$ could in principle be formulated, choosing $\lambda = 1$ appears natural and avoids the problem of selecting the tuning parameter.

INFERENCE

The asymptotic properties of fixed-effects estimators are complicated by the fact that, as $N \to \infty$, the dimension of the parameter $\bm{\alpha}_N$ tends to infinity. Unless $T \to \infty$, the individual effects $\alpha_i$ are estimated using a fixed number of observations. This is often referred to as the “incidental parameter” problem (NS; lancaster), which causes widely used estimators, such as maximum likelihood and M-estimators, to be inconsistent.

To develop the asymptotic theory of our estimator, we follow the recent panel data literature in econometrics and deal with the incidental parameter problem by considering asymptotic sequences where both $N$ and $T$ tend to infinity (e.g., Phillips:1999p733, hn, koenker2004, fv, ArellanoHahn2007, lamarche, and kato). Under this approximation, we show that our estimators are consistent but might have biases in the asymptotic distribution depending on the relative rate of convergence of $N$ and $T$. We apply the theory of M-estimators (e.g., newey), and use well-established results to handle the following non-standard features of our problem: (i) the estimators of $\bm{\theta}$, $\bm{\phi}$ and $\bm{\alpha}_N$ converge at different rates (e.g., radchenko; cs; masuda); and (ii) additional conditions are required on the relative growth rate of $N$ and $T$ (e.g., hn; fv; neweynotes).

Let $\bm{x}_{it} = (\bm{x}_{1i}^{ \mathrm{\scriptscriptstyle T} },\bm{x}_{2it}^{ \mathrm{\scriptscriptstyle T} })^{ \mathrm{\scriptscriptstyle T} }$, where $\bm{x}_{1i}$ contains the time-invariant components including the constant and $\bm{x}_{2it}$ contains the time-varying covariates. We use the following sufficient conditions to establish the identification of the parameters and derive the asymptotic properties of the estimators:

We use Assumptions (ref)(i)-(iv) to establish the identification of all the model parameters. Assumption (ref)(i) imposes that the model is correctly specified. It also requires $\alpha_i$ and $Y_{it} - \alpha_i$ to be conditionally independent across $(i,t)$, and to be independent of each other. We do not impose any sampling condition on the covariate sequences $\{(\bm{x}_{it}, \bm{z}_i) : 1 \leq i \leq N, 1 \leq t \leq T \}$, other than the existence of some limits. Assumptions (ref)(ii)-(iii) apply standard regularity conditions for parameter identification in quantile regression to our longitudinal model (e.g., ACF06). For example, Assumption (ref)(iii) imposes that the conditional quantile and density functions of $Y_{it} - \alpha_i$ and $\alpha_i$ are bounded. These conditions, together with a location normalization on the fixed effects in Assumption (ref)(i), guarantee that $\bm{\beta}(\cdot)$ and $\bm{\gamma}(\cdot)$ in the model (ref) are identified.\footnote{We normalize the mean of the fixed effects. Alternative normalizations on the median or other quantile of the fixed effects are also possible.} Then, Assumption (ref)(iv) pins down $\bm{\theta}^0$ and $\bm{\phi}^0$ from the system of linear equations $\bm{\beta}(\cdot) = \bm{\theta}^0 \bm{b}(\cdot)$ and $\bm{\gamma}(\cdot) = \bm{\phi}^0 \bm{c}(\cdot)$. Assumption (ref)(iv) provides a sufficient condition to guarantee existence and uniqueness of solution to the system from a subset of equations, which is easy to verify in practice. It can be replaced by any other existence and uniqueness condition.

Assumptions (ref)(v)-(vii) impose regularity conditions to derive the distribution of the estimators in large samples. The derivation relies on a general asymptotic expansion for fixed effects M-estimators given in Appendix A, which extend the results of hn and fv to estimators with mixed-rates asymptotics. Assumption (ref)(v) requires sufficient smoothness and bounded moments of the objective functions (ref) and (ref) and their partial derivatives, which are needed to carry out higher-order expansions of these functions. Assumption (ref)(vi) guarantees that all the terms of the expansions are well-defined. Finally, Assumption (ref)(vii) is a standard condition imposing that the limit Hessian matrices of the objective functions are non-singular.

Theorem (ref) shows that the parameters $\bm{\theta}^0$ and $\bm{\phi}^0$ are identified and their estimators $\hat{\bm{\theta}}$ and $\hat{\bm{\phi}}$ have a normal distribution in large samples with different rates of convergence. The large sample distribution of the plugin estimators of $\bm{\beta}^0(\cdot) = \bm{\theta}^0\bm{b}(\cdot)$ and $\bm{\gamma}^0(\cdot) = \bm{\phi}^0 \bm{c}(\cdot)$ can be obtained by the delta method. Let $\hat{\bm{\beta}}(u) = \hat{\bm{\theta}} \bm{b}(u)$ and $\hat{\bm{\gamma}}(v) = \hat{\bm{\phi}} \bm{c}(v)$, for any $u,v \in (0,1)$. Then, if $N = O(T)$, $$ \sqrt{NT}\left( \hat{\bm{\beta}}(u) - \bm{\beta}^0(u) + \frac{\operatorname{vec}_{d_x,d_{b}}^{-1}(\bar{\bm{H}}_{\bm{\theta}}^{-1} \bar{\bm{b}}_{\bm{\theta}}) \bm{b}(u)}{T} \right) \to_d \text{N}(0, \tilde{\bm{b}}(u)^{{ \mathrm{\scriptscriptstyle T} }}\bar{\bm{H}}_{\bm{\theta}}^{-1}\bar{\bm{\Omega}}_{\bm{\theta}}\bar{\bm{H}}_{\bm{\theta}}^{-1} \tilde{\bm{b}}(u)), $$ where $\operatorname{vec}_{d,k}^{-1}$ is the inverse vectorization operator that maps a $dk$-vector to a $d \times k$ matrix, i.e. $\operatorname{vec}_{d,k}^{-1}(\bm{v}) = \{[\operatorname{vec}(\bm{I}_k)]^{{ \mathrm{\scriptscriptstyle T} }} \otimes \bm{I}_d \} (\bm{I}_k \otimes \bm{v})$, $\bm{I}_n$ is the identity matrix of size $n$, and $\tilde{\bm{b}}(u) = \bm{b}(u) \otimes \bm{I}_{d_x}$. Similarly, if $N = O(T^2)$,

where $\tilde{\bm{c}}(v) = \bm{c}(v) \otimes \bm{I}_{d_z}$.

The rates of convergence of all the estimators agree with the square roots of the dimensions of the observations that are informative about the corresponding parameters. Thus, the rate is $\sqrt{NT}$ for $\bm{\theta}^0$ and $\bm{\beta}^0(u)$, and $\sqrt{N}$ for $\bm{\phi}^0$ and $\bm{\gamma}^0(v)$. All the estimators might suffer from bias in short panels due to the estimation of the fixed effects. The order of this bias is the inverse of the number of observations that are informative about each fixed effect, i.e. $T^{-1}$. Comparing the rates of convergence with the order of the bias, we can see that the biases of $\hat{\bm{\theta}}$ and $\hat{\bm{\beta}}(u)$ are negligible in the asymptotic distribution when $N/T \to 0$, whereas the biases of $\hat{\bm{\phi}}$ and $\hat{\bm{\gamma}}(v)$ are negligible when $N/T^2 \to 0$. These biases can be reduced by using analytical or jackknife corrections (e.g., hn, fv, or dj15). We provide consistent analytical estimators of the components of the biases and variances below.

We construct estimators of the components of the asymptotic distribution using sample analogs evaluated at the estimated value of the parameters, e.g., $\hat{u}_{it} = u_{it}(\hat{\bm{\theta}}, \hat \alpha_i)$ and $\hat{v}_{i} = v_{i}(\hat{\bm{\phi}}, \hat \alpha_i)$. Then,

where

The consistency of these estimators follows from the law of large numbers and consistency of $\hat \bm{\theta}$, $\hat \bm{\phi}$ and $\hat \alpha_i$, $1 \leq i \leq N$, together with the continuous mapping theorem, as all the components are continuous functions of the parameters.

GOODNESS-OF-FIT MEASURES AND MODEL SELECTION

Goodness-of-fit

To assess the model fit, we use the fact that, under the true model, $\hat u_{it}$ and $\hat v_i$ consistently estimate the realizations of the two independent uniform variables $U_{it}$ and $V_i$ that generated the data.

A graphical inspection of the joint and marginal distributions of $\hat u_{it}$ and $\hat v_i$ is always recommended. For a formal test, we suggest comparing the empirical distribution of $\hat u_{it}$ and $\hat v_i$, $\hat F_{\hat u \hat v}$, with the distribution of $U_{it}$ and $V_i$, given by \[F_{uv}(u,v) = uv.\] Following iqr (iqr, ctiqr), we compute a p-value for the null hypothesis $H_0:$ \{the model is correct\} using a Monte Carlo procedure:

After repeating steps 1-2 for a sufficient number of times, the p-value is computed as the empirical proportion of cases in which $D^* > D$. In step 1, it is also possible to take a random sample of clusters, and to resample the covariates' value within each cluster. To assess local fit, the test could be repeated within subsets of the original sample identified by specific values of the covariates.

In the implementation of the qrcm package, we chose $D$ to be the Kolmogorov-Smirnov statistic, $\sup_{u,v} |\hat F_{\hat u \hat v}(u, v) - uv|$. This testing procedure is usually reliable, as indicated by the simulation results reported in Section (ref).

Model selection

As suggested in Section (ref) and exemplified in the real-data example presented in Section (ref), it is usually possible to identify numerous alternative models that have a similar fit and are not rejected by a goodness-of-fit test. This can be explained by the fact that the same coefficient functions can be well approximated by different parametric functions.

Important criteria for model selection include parsimony, flexibility, and interpretability. Nested models can be compared by standard Wald test. Let $\hat{\bm{\alpha}}_N = (\hat \alpha_1, \ldots, \hat \alpha_N)$. To compare non-nested models, the value of $L_1(\hat\bm{\theta},\hat{\bm{\alpha}}_N)$ and $L_2(\hat\bm{\phi},\hat{\bm{\alpha}}_N)$ can be used to construct information criteria such as the AIC akaike and the BIC schwarz. These criteria were initially designed for likelihood-based estimators, but can be extended to estimators defined by the minimizer of a loss function. For example, a modification of BIC criterion for M-estimators has been described by machado, while koenker used AIC to compare quantile regression models. Consider the probability density function of the asymmetric Laplace distribution, $$f_p(y \mid \mu,\sigma) = \frac{p(1 - p)}{\sigma(p)}\exp\Bigl\{-\frac{\rho(y - \mu(p))}{2\sigma(p)}\Bigr\},$$ where $\mu(p)$ is a location parameter and corresponds to the $p$-th quantile of the distribution, while $\sigma(p)$ is a scale parameter. Although this distribution is not generally considered a plausible model, its log-likelihood has been used by numerous authors, including km1999 and lee, to obtain measures of goodness-of-fit for quantile regression. Simple algebra permits showing that

i.e., $\hat\bm{\theta}$ and $\hat\bm{\phi}$ minimize an “average” Laplace log-likelihood, in which $u$ and $v$ have been integrated away. After substituting $\hat \sigma_1 = L_1(\hat\bm{\theta}, \hat{\bm{\alpha}}_N)/(2NT)$ and $\hat \sigma_2 = L_2(\hat\bm{\phi}, \hat{\bm{\alpha}}_N)/(2N)$, we obtain the following AIC and BIC: $$ \textsc{AIC}_1 = \log L_1(\hat\bm{\theta}, \hat{\bm{\alpha}}_N) + \frac{q_1}{NT}, \hspace{0.2cm} \textsc{BIC}_1 = \log L_1(\hat\bm{\theta}, \hat{\bm{\alpha}}_N) + \frac{q_1\log{(NT)}}{2NT}, $$ $$ \textsc{AIC}_2 = \log L_2(\hat\bm{\phi}, \hat{\bm{\alpha}}_N) + \frac{q_2}{N}, \hspace{0.2cm} \textsc{BIC}_2 = \log L_2(\hat\bm{\phi}, \hat{\bm{\alpha}}_N) + \frac{q_2\log{(N)}}{2N} $$ where $q_1$ and $q_2$ are the number of non-zero elements of $\bm{\theta}$ and $\bm{\phi}$, respectively. Note that BIC$_1$ can be obtained from equation 2.3 of lee by replacing the loss of standard quantile regression with $L_1(\hat\bm{\theta}, \hat{\bm{\alpha}}_N)$.

The proposed criteria seem to work well in simulation (see Appendix C). However, they often tend to reward parsimony, possibly sacrificing goodness of fit. The testing procedure described in Section (ref) should always be used to perform a preliminary screening of the candidate models.

SIMULATION RESULTS

We analyze the performance of the estimators $\hat{\bm{\beta}}(u)$ and $\hat{\bm{\gamma}}(v)$ in finite samples through numerical simulations. In particular, we report the biases and standard errors of these estimators for different values of the dimensions $T$ and $N$ and the orders of the quantiles $u$ and $v$. We also evaluate the empirical size and power of the goodness of fit test.

We used the following design to generate the data: \[Y_{it} = \beta_0(U_{it}) + \beta_1(U_{it})x_{it} + \gamma_0(V_i) + \gamma_1(V_i)z_i\] where $x_{it} \sim \text{Beta}(2,2)$ and $z_i \sim \text{U}(0,3)$. In simulation 1, we defined:

where $\zeta(v)$ is the quantile function of a standard normal distribution. In simulation 2, we defined:

To fit the true model, we used $\bm{b}(u) = \left[1, -\log(1 - u), (u - 0.5)^3\right]^{{ \mathrm{\scriptscriptstyle T} }}$ and $\bm{c}(v) = \left[\zeta(v)\right]$ in simulation 1, and $\bm{b}(u) = \left[1, 1 - (1 - u)^{1/4}, u\right]^{{ \mathrm{\scriptscriptstyle T} }}$ and $\bm{c}(v) = \left[\log(1 - \log(1 - v))\right]$ in simulation 2. We ran $R = 1000$ Monte Carlo simulations, with $N = \{150, 300\}$ and $T = \{5, 10\}$. In Tables (ref) and (ref), we report the true value of $\bm{\beta}(\cdot)$ and $\bm{\gamma}(\cdot)$ at the quintiles, their average estimates, the empirical standard errors across simulations, and the average estimates of the asymptotic standard errors. Despite the incidental parameters problem, a small bias was found, even with small values of $T$. Also, as $T$ increased, the observed bias decreased rapidly as predicted by the asymptotic theory of Section (ref). The estimated standard errors were, on average, very close to their true values.

To assess the performance of the goodness-of-fit procedure described in Section (ref), we selected two nominal significance levels, $\alpha = 0.05$ and $\alpha = 0.10$, and computed the empirical probability of type I error ($\tilde{\alpha}$) and the power ($1 - \tilde{\beta}$) of the Kolmogorov-Smirnov goodness-of-fit test described in Section (ref). The power was estimated by the empirical probability to reject a misspecified model in which the quantile function was described by an incorrect basis function. In simulation 1, we incorrectly parametrized $\beta_1(u)$ as a linear function, $\beta_1(u) = \theta_{01} + \theta_{11}u$. In simulation 2, we incorrectly assumed that the individual effects have a logistic distribution, defined by $\bm{c}(v) = \left[\log(v/(1 - v))\right]$. Results are shown in the bottom rows of Tables (ref) and (ref). The risk of type I error was very close to its nominal level, and approached it as the value of $T$ increased. With small values of $N$ and $T$, the risk of type II error was relatively large, and the power was often less than $50\%$. However, with $N = 300$ and $T = 10$, and a nominal level of $0.10$ for rejection, the incorrect models were rejected in more than $90\%$ of cases in both scenarios.

Additional simulation results are reported in Appendix C, where we compare our estimator with koenker2004's\ (koenker2004) penalized fixed-effects quantile regression, and discuss the performance of the model selection criteria presented in Section (ref).

ANALYSIS OF NGAL DATA

We analyzed data from app, aiming to investigate the role of plasma neutrophil gelatinase-associated lipocalin (NGAL) as a marker of sepsis and acute kidney disfunction. The dataset included 139 patients admitted to the general intensive care unit at Karolinska University Hospital in Solna, Sweden, between August 2007 and November 2010. Baseline information was collected, and patients were classified daily as having sepsis or not. NGAL (mg/mL), procalcitonin (PCT), C-reactive protein (CRP), and creatinine changes relative to baseline ($\Delta{\text{creat}}$) were measured daily before discharge, for a total of 1317 plasma samples. After removing missing data, individuals with only one observation, and one patient with severe complications, the final sample included 135 patients for a total sample size of $\sum_{i = 1}^{135} T_i = 1263$. The number of observations per patient varied between $2$ and $38$, and more than $80\%$ of patients had $T_i \le 14$.

The goal of our analysis was to estimate conditional quantiles of NGAL, and in particular to measure its association with sepsis. The between-patient variability appeared to be very large, reflecting the presence of important differences in the initial health conditions. We formulated a regression model with the following predictors: a binary indicator of sepsis status, an indicator of $\Delta{\text{creat}} \ge 50$, age (centered at its median, 52 years, and divided by 10), an indicator of female gender, and time since hospitalization (weeks). Age and gender were cluster-invariant and were also included as level-2 predictors.

The response variable was log-transformed, which made it more plausible to define individual effects on the additive scale as in model ((ref)). The regression function was

We formulated a variety of models, in which $\beta_0(\cdot)$ and $\gamma_0(\cdot)$ were unbounded, while the other coefficients were modeled by bounded functions. To facilitate interpretation, we forced $\bm{\gamma}(0.5) = 0$, assigning the individual effects a zero-median distribution in which level-2 covariates only affect the scale parameter.

Selected modeling options are illustrated in Table (ref). Different models appeared to fit the data well, and were not rejected by the goodness-of-fit test described in Section (ref). The following model combined simplicity and flexibility, and was selected for illustrative purposes:

The level-1 and level-2 intercepts were described by different versions of the asymmetric Logistic distribution. The coefficient functions associated with level-1 covariates were a combination of linear and root-4 functions, while those of the level-2 predictors were assumed to be linear.

The p-value of the Kolmogorov-Smirnov test was $0.21$. To assess local fit, the test was repeated in subsamples with different values of the covariates (e.g., the females, those with $\Delta{\text{creat}} > 50$, etc.). No significant evidence of model misspecification was found.

All 27 model parameters are reported in Table (ref), while regression coefficients at selected quantiles are summarized in Table (ref). We represent graphically the quantile regression coefficient functions in Figure 2, where we also report a “nonparametric” fit obtained by modeling all coefficients as piecewise linear functions with knots at the deciles.

Results showed that the distribution of the individual effects was almost symmetric (as suggested by the fact that $\hat\phi_{01} \simeq -\hat\phi_{02}$) and that its variance was not significantly affected by cluster-level predictors. Instead, all predictors apart from gender appeared to be associated with the level-1 response. In particular, the coefficients associated with sepsis, $\Delta{\text{creat}}_{it} > 50$ and age were consistently positive at all quantiles, while the coefficient of time was always negative. The sepsis status was associated with a percentile difference of about 0.12 at quantiles 0.2, 0.4, 0.6, 0.8. As shown by Figure 2, an even larger percentile difference was found at quantiles above 0.8.

CONCLUSIONS

We introduced a general framework for longitudinal quantile regression, extending the work of iqr (iqr, ctiqr) on quantile regression coefficients modeling. We defined a two-level quantile function in which both the “within” and the “between” part of the distribution are described by a quantile regression model. This allows to investigate how covariates affect not only the level-1 response, but also the distribution of the individual effects, which is generally overlooked in the existing literature on longitudinal quantile regression. Identification is achieved by modeling the coefficient functions parametrically, and estimation is carried out by minimizing a smooth objective function.

The proposed method is computationally simple and can be viewed as a special type of penalized fixed-effects estimator that presents important elements of novelty of its own. The penalty term carries information on the parameters that describe the conditional distribution of the individual effects. This permits estimating both level-1 and level-2 parameters, as in random-effects models, but allows carrying out estimation and inference using fixed-effects techniques. Moreover, it avoids the problem of choosing a tuning constant as in standard $\ell_1$- or $\ell_2$-penalization. The described form of penalized fixed-effects method is not limited to a quantile regression framework and could be applied to different estimation problems.

The proposed modeling framework can be generalized in different directions. An interesting possibility is to include multiple individual effects as in random-slope models. In our framework, individual effects are represented by a pure location shift as in koenker2004, geraci (geraci, geraci2), and canay. Using the proposed penalized fixed-effects approach, it is relatively simple to incorporate not only an individual intercept, $\alpha_i$, but also a set of individual slopes, say $\{\delta_{1i}, \delta_{2i}, \ldots\}$. This, however, would typically result in cumbersome computation and, unless $N$ and $T$ are sufficiently large, would probably undermine model identifiability. Using koenker2004's\ (koenker2004) words: “At best we may be able to estimate an individual specific location-shift effect, and even this may strain credulity”.

Another interesting extension is represented by varying-coefficients models (e.g., hastie, fan1, fan2, chiang,kim2) that could be implemented by allowing the level-1 regression coefficients to be functions of time. A possible approach is to describe the coefficients, say $\bm{\beta}(u, t)$, using tensor products of splines. Finally, the proposed method could be used to estimate static and dynamic quantile autoregressive models (e.g., arellano).

An important problem that has not been discussed in the paper is represented by quantile crossing, occurring when either $\bm{x}_{it}^{ \mathrm{\scriptscriptstyle T} }\bm{\beta}^{\prime}(u \mid \hat\bm{\theta}) < 0$ or $\bm{z}_{i}^{ \mathrm{\scriptscriptstyle T} }\bm{\gamma}^{\prime}(v \mid \hat\bm{\phi}) < 0$. One may want to determine in advance which values of the parameters $\bm{\theta}$ and $\bm{\phi}$ would ensure that no crossing occurs, i.e., that $\bm{x}_{it}^{ \mathrm{\scriptscriptstyle T} }\bm{\beta}(\cdot \mid \hat\bm{\theta})$ and $\bm{z}_{i}^{ \mathrm{\scriptscriptstyle T} }\bm{\gamma}(\cdot \mid \hat\bm{\phi})$ are monotonically increasing functions. This is only possible in very simple models with few covariates, or in presence of restrictive assumptions. However, simulation evidence suggests that parametric models are relatively immune to quantile crossing, compared with the “nonparametric” approaches based on ordinary quantile regression. Additionally, the parametric structure makes it particularly simple to verify crossing, taking advantage of the closed-form analytical expression of the quantile function, and admits the application of monotonization methods such as the rearrangement of cfg10 to produce increasing estimates of conditional quantiles.

This paper is accompanied by an R package qrcm, that includes a function named iqrL that performs model fitting, and a variety of auxiliary functions for prediction, plotting, and goodness-of-fit assessment. The documentation contains a rich set of examples, and can serve as tutorial for the practitioners. The package is available upon request to the authors.