A New Family of Error Distributions for Bayesian Quantile Regression
Abstract
We propose a new family of error distributions for modelbased quantile regression, which is constructed through a structured mixture of normal distributions. The construction enables fixing specific percentiles of the distribution while, at the same time, allowing for varying mode, skewness and tail behavior. It thus overcomes the severe limitation of the asymmetric Laplace distribution – the most commonly used error model for parametric quantile regression – for which the skewness of the error density is fully specified when a particular percentile is fixed. We develop a Bayesian formulation for the proposed quantile regression model, including conditional lasso regularized quantile regression based on a hierarchical Laplace prior for the regression coefficients, and a Tobit quantile regression model. Posterior inference is implemented via Markov Chain Monte Carlo methods. The flexibility of the new model relative to the asymmetric Laplace distribution is studied through relevant model properties, and through a simulation experiment to compare the two error distributions in regularized quantile regression. Moreover, model performance in linear quantile regression, regularized quantile regression, and Tobit quantile regression is illustrated with data examples that have been previously considered in the literature.
Keywords.
Asymmetric Laplace distribution; Markov chain Monte Carlo;
Regularized quantile
regression; Skew normal distribution; Tobit quantile regression.
1 Introduction
Quantile regression offers a practically important alternative to traditional mean regression, and forms an area with a rapidly increasing literature. Parametric quantile regression models are almost exclusively built from the asymmetric Laplace (AL) distribution the density of which is
(1) 
where , with denoting the indicator function. Here, is a scale parameter, , and corresponds to the th percentile, . Hence, a model for th quantile regression can be developed by expressing as a function of available covariates , for instance, yields a linear quantile regression structure. Note that maximizing the likelihood with respect to under an AL response distribution corresponds to minimizing for the check loss function, , used for classical semiparametric estimation in linear quantile regression (Koenker, 2005).
The AL distribution is receiving increasing attention in the Bayesian literature, originating from work on inference for linear quantile regression (Yu and Moyeed, 2001; Tsionas, 2003). Particularly relevant to the Bayesian framework are the different mixture representations of the distribution (Kotz et al., 2001), which have been exploited to construct posterior simulation algorithms (Kozumi and Kobayashi, 2011), as well as to explore different modeling scenarios; see, for instance, Lum and Gelfand (2012) and Waldmann et al. (2013).
However, the AL distribution has substantial limitations as an error model for quantile regression. Most striking is that the skewness of the error density is fully determined when a specific percentile is chosen, that is, when is fixed. In particular, the error density is symmetric in the case of median regression, since for , the AL reduces to the Laplace distribution. Moreover, the mode of the error distribution is at zero, for any , which results in rigid error density tails for extreme percentiles.
The literature includes Bayesian nonparametric models for the error distribution in the special case of median regression (Walker and Mallick, 1999; Kottas and Gelfand, 2001; Hanson and Johnson, 2002) and in general quantile regression (Kottas and Krnjajić, 2009; Reich et al., 2010). The Bayes nonparametrics literature has also explored inference methods for simultaneous quantile regression (Taddy and Kottas, 2010; Tokdar and Kadane, 2012; Reich and Smith, 2013). However, work on parametric alternatives to AL quantile regression errors is limited, and the existing models do not overcome all the limitations discussed above. For instance, although the class of skew distributions studied in Wichitaksorn et al. (2014) includes the AL as a special case, it shares the same restriction with the AL as a quantile regression error model in that it has a single parameter that controls both skewness and percentiles. Zhu and ZindeWalsh (2009) and Zhu and Galbraith (2011) explored the family of asymmetric exponential power distributions, which does not include the AL distribution. For a fixed probability , the density function has four free parameters and allows for different decay rates in the left and the right tails. However, similar to the AL, the mode of the distribution is fixed at the quantile by construction.
More flexible parametric quantile regression error models are arguably useful both to expand the inferential scope of the asymmetric Laplace in the standard quantile regression setting, as well as to provide building blocks for model development under more complex data structures. The limited scope of results in this direction may be attributed to the challenge of defining sufficiently flexible distributions that are parameterized by percentiles and, at the same time, allow for practicable modeling and inference methods.
Seeking to fill this gap, we propose a new family of distributions that is parameterized in terms of percentiles, and overcomes the restrictive aspects of the AL distribution. The distribution is developed constructively through an extension of an AL mixture representation. In particular, we introduce a shape parameter to obtain a distribution that has more flexible skewness and tail behaviour than the AL distribution, while retaining it as a special case of the new model. The latter enables connections with the check loss function which are useful in studying the utility of the new model in the context of regularized quantile regression. Owing to its hierarchical mixture representation, the proposed distribution preserves the important feature of ready to implement posterior inference for Bayesian quantile regression.
In Section 2, we develop the new distribution and discuss its properties relative to the AL distribution. In Section 3, we formulate the Bayesian quantile regression model, including a prior specification for the regression coefficients that encourages shrinkage resulting in regularized quantile regression, and a Tobit quantile regression formulation. In Section 4, we present results from a simulation study to compare the performance of the AL and the proposed distribution in regularized quantile regression. The methodology is illustrated with three data examples in Section 5, focusing again on comparison with the AL quantile regression model. Finally, Section 6 concludes with a summary and discussion of possible extensions.
2 The generalized asymmetric Laplace distribution
The construction of the new distribution is motivated by the most commonly used mixture representation of the AL density. In particular,
(2) 
where and . Moreover, denotes the normal distribution with mean and variance , and denotes the exponential distribution with mean . We use such notation throughout to indicate either the distribution or its density, depending on the context.
The mixture formulation in (2) enables exploration of extensions to the AL distribution. Extending the mixing distribution is not a fruitful direction in terms of evaluation of the intergal, and, more importantly, with respect to fixing percentiles of the resulting distribution. However, both goals are accomplished by replacing the normal kernel in (2) with a skew normal kernel (Azzalini, 1985). In its original parameterization, the skew normal density is given by , where and denote the density and distribution function, respectively, of the standard normal distribution. Here, is a location parameter, a scale parameter, and the skewness parameter. Key to our construction is the fact that the skew normal density can be written as a location normal mixture with mixing distribution given by a standard normal truncated on (Henze, 1986). More specifically, reparameterize to , where and , such that and . Then, , where denotes the standard normal distribution truncated over .
The proposed model, referred to as generalized asymmetric Laplace (GAL) distribution, is built by adding a shape parameter, , to the mean of the normal kernel in (2) and mixing with respect to a variable. More specifically, the full mixture representation for the density function, , of the new distribution is as follows
(3) 
Note that, integrating over in (3), the GAL density can be expressed in the form of (2) with the kernel replaced with a skew normal kernel, which, in its original parameterization, has location parameter , scale parameter , and skewness parameter . Evidently, when , reduces to the AL density.
To obtain the GAL density, we integrate out first and then in (3). The integrand of can be recognized as the kernel of a generalized inverseGaussian density. Therefore, integrating out , we obtain . This integral involves a normal density kernel, but care is needed with the limits of integration which depend on the sign of and of . Combining the resulting expressions from all possible cases, we obtain that for , the GAL density is given by
(4)  
where , , , with . The relatively complex form of the density in (4) is not an obstacle from a practical perspective, since its hierarchical mixture representation facilitates study of model properties and Markov chain Monte Carlo posterior simulation.
There is a direct link between the GAL distribution and the th quantile for any ; note that parameter no longer corresponds to the cumulative probability at the quantile for . When , the distribution function of (4) at is given by . Hence, letting , the distribution function becomes,
We use above, since this is the general form of that applies also in the case.
Note that, for , , where . The function is monotonically increasing in , since . Moreover, , and . Therefore, for , and thus is monotonically increasing in . Since is an even function, it also obtains that it is monotonically decreasing in .
Consider now setting . Then, the fact that is decreasing in combined with , imply that for each in the domain that respects the condition of and , there is a unique solution of that ensures , and subsequently a unique based on . For , setting and letting leads to the same argument.
The above connection between and suggests that by reparameterization with desired and , we can derive a new family of distributions with the percentile for fixed given by , and with an additional shape parameter . For , the density, , of such quantilefixed GAL distribution is
(5) 
where , , , and . Parameter has bounded support over interval , where is the negative root of and is the positive root of . For instance, takes values in , and when , and , respectively. When , the density reduces to the AL density, which is also a limiting case of (5). The density function is continuous for all possible values.
The quantilefixed GAL distribution has three parameters, , and . Note that has density if and only if has density . Hence, similarly to the AL distribution, is a location parameter and is a scale parameter. The new shape parameter enables the extension relative to the quantilefixed AL distribution. As demonstrated in Figure 1, controls skewness and tail behaviour, allowing for both left and right skewness when the median is fixed, as well as for both heavier and lighter tails than the asymmetric Laplace, the difference being particularly emphatic for extreme percentiles. Moreover, as varies, the mode is no longer held fixed at ; it is less than when and greater than when . The above attributes render the proposed distribution substantially more flexible than the AL distribution.
Finally, we note that parameter satisfies likelihood identifiability. Consider the locationscale standardized density, , which is effectively the model for the errors in quantile regression. Then, assume , for all . Given that parameter controls the mode of the density, this implies that and must have the same sign. Working with either of the two cases (that is, and or and ) in expression (5), we arrive at , which, based on the monotonicity of function , implies .
3 Bayesian quantile regression with GAL errors
3.1 Inference for linear quantile regression
Consider continuous responses and the associated covariate vectors , for . The linear quantile regression model is set up as , where the arise independently from a quantilefixed GAL distribution with . Owing to the mixture representation of the new distribution, the model for the data can be expressed hierarchically as follows
(6) 
where , and and are the functions of given in (2). Since is a function of and , , and are all functions of parameter . The Bayesian model is completed with priors for , and . Here, we assume a normal prior for and an inversegamma prior for , with mean provided . For any specified , is defined over an interval with fixed finite endpoints, and thus a natural prior for is given by a rescaled Beta distribution, with the uniform distribution available as a default choice.
The augmented posterior distribution, which includes the and the , can be explored via a Markov chain Monte Carlo algorithm based on Gibbs sampling updates for all parameters other than . As in Kozumi and Kobayashi (2011), we set , . Then, the posterior simulation method is based on the following updates.

Sample from , with covariance matrix and mean vector .

For each , sample from a generalized inverseGaussian distribution, , where and , with density given by .

For each , sample from a normal distribution truncated on , where and .

Sample from a distribution, where , , and .

Update with a MetropolisHasting step, using a normal proposal distribution on the logit scale over .
Based on the hierarchical model structure, the posterior predictive error density can be expressed as , and thus estimated through Monte Carlo integration, using the posterior samples of .
3.2 Quantile regression with regularization
Since the GAL distribution is constructed through modifying the mixture representation of the AL distribution, it retains some of the interesting properties of the AL distribution. In particular, working with the hierarchical representation of the GAL distribution, we are able to retrieve an extended version of the check loss function which corresponds to asymmetric Laplace errors.
Consider the collapsed posterior distribution, , that arises from (3.1) by marginalizing over the . Then, the corresponding posterior full conditional for can be expressed as
where is the prior density for , , and , with the probability associated with the specified quantile modeled through . Hence, ignoring the prior contribution, finding the mode of the posterior full conditional for is equivalent to minimizing with respect to the adjusted loss function ; note that in the special case with asymmetric Laplace errors, that is, for , this reduces to the check loss function with .
Based on the above structure, the positivevalued latent variables can be viewed as responsespecific weights that are adjusted by realvalued coefficient , which is fully specified through the shape parameter . The result is the realvalued, responsespecific terms , which reflect on the estimation of the effect of outlying observations relative to the AL distribution. A promising direction to further explore this structure is in the context of variable selection. For instance, Li et al. (2010) study connections between different versions of regularized quantile regression and different priors for , working with asymmetric Laplace errors. The main example is lasso regularized quantile regression, which can be connected to the Bayesian asymmetric Laplace error model through a hierarchical Laplace prior for . We consider this prior below extending the AL error distribution to the proposed GAL distribution. The perspective we offer may be useful, since it can be used to explore regularization adjusting the loss function, through the response distribution, in addition to the penalty term, through the prior for the regression coefficients.
Here, we denote by the dimensional vector of regression coefficients excluding the intercept . Then, the Laplace conditional prior structure for is given by
The second expression above utilizes the normal scale mixture representation for the Laplace distribution, which has been exploited for posterior simulation in the context of lasso mean regression (Park and Casella, 2008). Moreover, to facilitate Markov chain Monte Carlo sampling, we reparameterize in terms of and place a gamma prior on . The lasso regularized version of model (3.1) is completed with a normal prior for , and with the priors for the other parameters as given in Section 3.1. The posterior simulation algorithm is the same with the one described in Section 3.1 with the exception of the updates for the , , and for . Using the mixture representation of the Laplace prior, each can be sampled from a normal distribution, whereas has a gamma posterior full conditional distribution.
3.3 Tobit quantile regression
Tobit regression offers a modeling strategy for problems involving range constraints on the response variable (Amemiya, 1984). The standard Tobit regression model can be viewed in the context of censored regression where the responses are left censored at a threshold ; without loss of generality, we take . The responses can be written as , where are the observed values and are latent if . In the context of quantile regression, Yu and Stander (2007) and Kozumi and Kobayashi (2011) applied the ALbased model to the latent responses . Here, we consider the Tobit quantile regression setting with GAL errors.
Consider a data set of observations on covariates and associated responses , where consists of positivevalued observed responses with the remaining responses censored from below at . Assuming the GAL distribution for the latent responses, the likelihood can be expressed as . Using data augmentation (Chib, 1992), let be the unobserved (latent) responses corresponding to the data points that are left censored at . Then, using again the hierarchical representation of the GAL distribution, the joint posterior distribution that includes can be written as
where denotes the prior for the model parameters, and . Here, denotes a truncated normal on , and an exponential distribution with mean .
Regarding posterior inference, the posterior full conditional for each auxiliary variable is given by a truncated normal distribution. And, given the augmented data , the model parameters and the latent variables can be sampled as before.
Although results are not reported here, we have tested the posterior simulation algorithm on simulated data sets based on GAL errors, with observations and a censoring rate that ranged from to . Under this scenario, the posterior distributions successfully captured the true values of all parameters in their 95% credible intervals.
4 Simulation study
Here, we present results from a simulation study designed to compare the lasso regularized quantile regression models with AL and GAL errors. We follow a standard simulation setting from the literature regarding the linear regression component (Tibshirani, 1996; Zou and Yuan, 2008; Li et al., 2010), varying the extent of sparsity in the true vector. For the underlying datagenerating error distributions, we consider four scenarios with different types of skewness and tail behavior. For model comparison, we evaluate the accuracy in variable selection, inference for the regression function, and posterior predictive performance, using relevant assessment criteria. Overall, the GALbased quantile regression model performs better in variable selection and prediction accuracy and it is more robust to nonstandard error distributions, particularly for extreme quantiles. The two models yield comparable results in the case of median regression.
4.1 Simulation settings
We consider synthetic data generated from linear quantile regression settings, with , and to study model performance for both extreme and more central percentiles. The rows of the design matrix were generated independently from an 8dimensional normal distribution with zero mean vector and covariance matrix with elements , for . We present detailed results from a relatively sparse case for the vector of regression coefficients, . In Section 4.3, we briefly discuss results form two other scenarios for corresponding to a dense and a very sparse case.
Data were simulated under four different error distributions:

, with chosen such that the th quantile is 0.

, with chosen such that the th quantile is 0.

, with chosen such that the th quantile is 0.

, with and chosen such that the th quantile is 0. To generate the errors, we first sample from a generalized Pareto distribution, then take the logarithm. Based on the parameterization in Embrechts et al. (1997), the density function of the errors is given by , for .
The normal and Laplace error distributions are symmetric about zero under median regression. The parameters of the twocomponent normal mixture are selected such that the resulting error distribution is skewed. Finally, the logtransformed generalized Pareto distribution is included to study model performance under an error density which is both skewed and does not have exponential tails.
For each setting of the simulation study, we generated 100 data sets, each with observations for training the models and another for testing predictions.
4.2 Criteria for comparison
We consider a number of criteria to assess different aspects of model performance. Since Bayesian lasso regression only shrinks the covariate effects, we consider a threshold on the effect size for the purpose of variable selection. Following Hoti and Sillanpää (2006), we calculate the standardized effects as , , where is the standard deviation of predictor and is the standard deviation of the response. For each posterior sample, if the standardized effect is greater than 0.1 in absolute value, we consider the predictor as included. We count the number of correct inclusion and exclusions (CIE) in the posterior sample and divide it by to normalize it to a number between 0 and 1. By averaging over all the posterior samples, we obtain the mean standardized CIE for each simulated data set.
To assess predictive performance for the regression function, we calculate the mean check loss on the test data points, defined as: , where is the posterior mean estimate from the training data. The mean check loss resembles the standard mean squared error criterion, which is commonly used for evaluating prediction with crossvalidation.
Finally to assess model fitting taking into account predictive uncertainty, we apply the posterior predictive loss criterion from Gelfand and Ghosh (1998). This criterion favors the model that minimizes , where is a goodnessoffit term, and is a penalty term for model complexity. Here, , and and are the mean and variance under model of the posterior predictive distribution for replicated response with corresponding covariate . We also consider the generalized version of the criterion based on the check loss function, under which . For this generalized criterion, the goodnessoffit term can be defined by and the penalty term by , since the check loss function is convex in , and thus ; see Gelfand and Ghosh (1998) for details on defining the model comparison criterion under loss functions different from quadratic loss.
4.3 Results
We used the same hierarchical Laplace prior for under the AL and GAL models, with a gamma prior for with prior mean and variance . Such prior specification is relatively noninformative in the sense that it does not favor shrinkage for the regression coefficients, resulting in marginal prior densities for each that place substantial probability mass away from . The shape parameter of the GAL error distribution was assigned a uniform prior. Results under both models and for each simulated data set are based on 5,000 posterior samples, obtained after discarding the first 50,000 iterations of the Markov chain Monte Carlo sampler and then retaining one every 20 iterations.
Within each simulation scenario, we summarize results from the 100 data sets using the median and standard deviation (SD) of the values for the performance assessment criteria discussed in Section 4.2. Results are reported in Table 1 through Table 4, where we use boldface to indicate the model supported by the particular criterion under each setting.
Overall, the lasso regularized Bayesian quantile regression model performs better under the GAL error distribution. The GALbased model includes/excludes correct regression coefficient values more often than the AL model for almost all combinations of and error distributions (Table 1). It also results in a lower median mean check loss for the test data in most cases, demonstrating better performance in the prediction of the regression function (Table 2). Note that, for both types of assessment in Tables 1 and 2, the GALbased model produces better results across all error distributions for , and, with the exception of one case, when . Results are generally more balanced in the median regression setting, although the GAL model fares better in all cases for which the underlying error distribution is skewed.
Error distribution  
logtransformed  
Model  Normal  Laplace  Normal mixture  generalized Pareto  
0.05  GAL  0.848 (0.063)  0.633 (0.083)  0.911 (0.042)  0.893 (0.052) 
AL  0.746 (0.099)  0.534 (0.087)  0.817 (0.075)  0.840 (0.081)  
0.25  GAL  0.851 (0.049)  0.728 (0.060)  0.918 (0.048)  0.896 (0.050) 
AL  0.843 (0.069)  0.700 (0.068)  0.913 (0.060)  0.900 (0.051)  
0.50  GAL  0.848 (0.052)  0.738 (0.065)  0.909 (0.049)  0.897 (0.055) 
AL  0.850 (0.056)  0.737 (0.065)  0.905 (0.050)  0.870 (0.061) 
Error distribution  
logtransformed  
Model  Normal  Laplace  Normal mixture  generalized Pareto  
0.05  GAL  0.340 (0.083)  1.073 (0.391)  0.224 (0.060)  0.268 (0.081) 
AL  0.523 (0.130)  1.709 (0.485)  0.375 (0.101)  0.388 (0.114)  
0.25  GAL  0.325 (0.086)  0.676 (0.199)  0.225 (0.071)  0.265 (0.080) 
AL  0.360 (0.096)  0.778 (0.215)  0.257 (0.076)  0.274 (0.082)  
0.50  GAL  0.323 (0.092)  0.642 (0.208)  0.235 (0.064)  0.262 (0.081) 
AL  0.322 (0.095)  0.624 (0.207)  0.237 (0.063)  0.294 (0.089) 
For each simulation setting, Table 3 includes the values for the posterior predictive loss criterion with quadratic loss (under , such that ), and Table 4 shows the generalized criterion under check loss. Both versions of the posterior predictive loss criterion support the GAL model when , with differences in values between the two models that are substantially larger than for the other two values of . This reinforces the earlier findings on the potential benefits of the GAL error distribution for extreme percentiles. With the exception of one case under the check loss version of the criterion, the GALbased model is also favored when , whereas results are more mixed in the median regression case.
Finally, although detailed results are not reported here, the simulation study included two more settings for , a dense case with all regression coefficients equal to , and a very sparse case with . The conclusions were overall similar, in particular, the GAL model outperformed the AL model for essentially all combinations of underlying error distribution and value of or . Again, in the median regression case, the distinction between the two models was less clear for the normal, Laplace and normal mixture datagenerating distributions, although the GAL model performed better under all criteria for the setting corresponding to the logtransformed generalized Pareto distribution.
Error distribution  
logtransformed  
Model  Score  Normal  Laplace  Normal mixture  generalized Pareto  
0.05  GAL  1231 (193)  9799 (2483)  653 (112)  1273 (270)  
832 (126)  7046 (1546)  429 (71)  1053 (267)  
2092 (312)  16860 (3839)  1085 (181)  2319 (531)  
AL  3359 (799)  30308 (10763)  1839 (405)  2782 (664)  
952 (165)  8659 (2304)  534 (93)  1168 (279)  
4357 (933)  38766 (12676)  2398 (487)  4020 (873)  
0.25  GAL  1085 (206)  6977 (1607)  608 (95)  1445 (273)  
830 (146)  6897 (1606)  444 (66)  1105 (264)  
1882 (343)  13884 (3115)  1055 (154)  2552 (511)  
AL  1630 (303)  11503 (2727)  884 (148)  1516 (260)  
865 (154)  7395 (1742)  464 (71)  1113 (263)  
2499 (448)  18916 (4349)  1352 (215)  2600 (487)  
0.50  GAL  1283 (205)  7600 (1676)  694 (97)  1189 (217)  
813 (132)  6459 (1509)  424 (60)  1089 (245)  
2111 (328)  14076 (3101)  1121 (152)  2283 (415)  
AL  1177 (191)  7256 (1572)  634 (87)  1318 (247)  
818 (134)  6431 (1509)  426 (60)  1107 (255)  
2008 (318)  13667 (3019)  1058 (143)  2415 (483) 
Error distribution  
logtransformed  
Model  Normal  Laplace  Normal mixture  generalized Pareto  
0.05  GAL  174.2 (13.6)  507.0 (67.8)  122.8 (11.3)  178.5 (17.4) 
AL  209.3 (21.7)  605.4 (70.8)  148.6 (17.3)  200.2 (20.5)  
0.25  GAL  169.5 (15.9)  443.9 (47.1)  126.2 (9.5)  188.0 (17.5) 
AL  178.0 (15.7)  451.4 (45.8)  129.0 (9.8)  185.5 (17.3)  
0.50  GAL  175.7 (13.4)  444.7 (48.0)  127.4 (8.9)  178.6 (16.1) 
AL  172.6 (13.4)  438.5 (47.5)  125.2 (8.7)  183.6 (18.1) 
5 Data examples
In this section, we consider three data examples to illustrate the Bayesian quantile regression models developed in Sections 3.1, 3.2, and 3.3. The main emphasis is on comparison of inference results between models based on the GAL distribution and those assuming an AL distribution for the errors.
We have implemented both models with priors for their parameters that result in essentially the same prior predictive error densities. The two models were applied with the same prior distributions for and . More specifically, for the data sets of Sections 5.1 and 5.3, we used a prior for the vector of regression coefficients, and an prior for the scale parameter . For the data example of Section 5.2, we used a prior for the intercept, and the same conditional Laplace prior for the remaining regression coefficients with the simulation study (see Section 4.3). Finally, a uniform prior was placed on the shape parameter of the GAL error distribution. For all data examples, the posterior densities for model parameters were fairly concentrated relative to the corresponding prior densities.
5.1 ImmunoglobulinG data
We illustrate the proposed model, referred to as model , with a data set commonly used in additive quantile regression; see, for instance, Yu and Moyeed (2001). The analysis focuses on comparison with the simpler model based on asymmetric Laplace errors, referred to as model . The data set contains the immunoglobulinG concentration in grams per litre for children aged between 6 months and 6 years. As in earlier applications of quantile regression for these data, we use a quadratic regression function to model five quantiles, corresponding to , of immunoglobulinG concentration against covariate age ().
The two models result in different posterior predictive error densities, especially for extreme percentiles; see Figure 2. At , under the AL model, both the shape and the skewness of the error distribution are predetermined by and the mode is forced to be 0, resulting in a rigid heavy left tail. The effect of this overly dispersed tail can be observed in the inference for the quantile regression function (Figure 3). The GAL model, on the contrary, yields an error density that has a much thinner left tail, concentrating more of its probability mass around the mode, which is not at 0. Figure 2 shows also the posterior densities for shape parameter , under a uniform prior in all cases. For all three quantile regressions, the 95% posterior credible interval for does not include the value of 0, which corresponds to asymmetric Laplace errors. Median regression is the only case where 0 is within the effective range of the posterior distribution for .
Bayesian information criterion  

Quantile  Model  loglikelihood  BIC 
For formal model comparison, we compute the Bayesian information criterion (BIC), the posterior predictive loss criterion with quadratic loss, and the generalized posterior predictive loss criterion under the check loss. The Bayesian information criterion favors the new model at all five quantiles; see Table 5. Under the posterior predictive loss criterion (Table 6), the two models are comparable in the case of median regression, with model preferred. In all other cases, model is favored by both versions of the model comparison criterion. The improvement in performance over the AL model is particularly conspicuous at the two extreme percentiles. This is in agreement with the difference in the posterior predictive error densities for , reported in Figure 2.
Posterior predictive loss criterion  
Quadratic loss  Check loss  
Quantile  Model  P  G  D  P  G  D 
5.2 Boston housing data
We apply the lasso regularized quantile regression model to the realty price data from the Boston Standard Metropolitan Statistical Area (SMSA) in 1970 (Harrison and Rubinfeld, 1978). The data set contains 506 observations. We take the logtransformed corrected median value of owneroccupied housing in USD 1000 (LCMEDV) as the response, and consider the following predictors: point longitudes in decimal degrees (LON), point latitudes in decimal degrees (LAT), per capita crime (CRIM), proportions of residential land zoned for lots over 25000 square feet per town (ZN), proportions of nonretail business acres per town (INDUS), a factor indicating whether tract borders Charles River (CHAS), nitric oxides concentration (parts per 10 million) per town (NOX), average numbers of rooms per dwelling (RM), proportions of owneroccupied units built prior to 1940 (AGE), weighted distances to five Boston employment centers (DIS), index of accessibility to radial highways per town (RAD), fullvalue propertytax rate per USD 10,000 per town (TAX), pupilteacher ratios per town (PTRATIO), transformed African American population proportion (B), and percentage values of lower status population (LSTAT).
We consider quantiles of and and compare the maximum a posteriori estimates (MAP) of regression coefficients, along with 95% credible intervals, for standardized covariates under the lasso regularized quantile regression models with AL and GAL errors (Figure 4 and 5). For both quantiles, the widths of the 95% credible intervals for the regression coefficients are overall comparable between the two models, but the posterior point estimates can be quite different. For instance, under the 10th quantile regression, the GAL model shrinks the effects of per capita crime (CRIM) and proportytax rate (TAX) to a greater extent compared to the AL model. Similar patterns can be observed for index of accessibility to radial highways (RAD) for the 90th quantile. Moreover, the two models reach different conclusions on the effect of latitude (LAT) for the 10th percentile. Although the posterior point estimates suggest a higher housing price as latitude increases adjusting for all other covariates, the 95% credible interval under the GAL model includes 0, whereas the one under the AL model does not.
Focusing on inference under the GAL error distribution, we note that, although the model selected some common variables for the two quantiles, there is also some discrepancy. For instance, each of higher proportions of residential land zoned for lots over 25000 square feet per town (ZN) and having tracts bordering Charles river (CHAS) increase the price at the 90% percentile, while higher nitrogen oxide value (NOX) has a negative influence on the 90% percentile price. However, none of these covariates have a significant effect on the realty value at the 10% percentile.
Finally, we notice that for both the 10th and 90th quantile regression, 0 is far away from the endpoints of the 95% credible interval for the GAL model shape parameter . This suggests that asymmetric Laplace errors are not suitable for this particular application. This is further supported by the results for the posterior predictive loss criterion reported in Table 7.
Posterior predictive loss criterion  
Quadratic loss  Check loss  
Quantile  Model  P  G  D  P  G  D 
5.3 Labor supply data
We illustrate the Tobit quantile regression model with the female labor supply data from Mroz (1987), which was taken from the University of Michigan Panel Study of Income Dynamics for year 1975. The data set includes records on the work hours and other relevant information of married white women aged between 30 and 60 years old. Of the women, worked at some time during 1975, with the corresponding fully observed responses given by the wife’s work hours (in 100 hours). For the remaining women, the observed zero work hours correspond to negative values for the latent “labor supply” response. We use the quantile regression function considered in Kozumi and Kobayashi (2011), where an ALbased Tobit quantile regression model was applied to the same data set. The linear predictor includes an intercept, income which is not due to the wife (nwifeinc), education of the wife in years (educ), actual labor market experience in years (exper) and its quadratic term (expersq), age of the wife (age), number of children less than 6 years old in household (kidslt6), and number of children between ages 6 and 18 in household (kidsge6). We compare the results from the Bayesian Tobit quantile regression model assuming AL errors (model ) and GAL errors (model ).
Table 8 summarizes the posterior distribution of under the GAL model, and presents results from criterionbased comparison of the two models for , and . Since there is censoring in the data, we use the revised BIC from Volinsky and Raftery (2000). In all three cases, the 95% credible interval for excludes , and the GALbased model is associated with lower BIC values. The results support the GALbased model more emphatically for the extreme percentiles than for median regression.
Quantile  Model  Mean (95% CrI) for  likelihood  BIC 

(, )  
(, )  
(, ) 
Figure 6 shows the posterior distributions of labor supply quantiles corresponding to , and for women with 0, 1, 2 and 3 children less than 6 years old. For all other predictors, we use the median values from the data as input values to represent an average wife. As the number of young children increases, the AL model estimates the 5th quantile and the median of labor supply of an average wife to be closer to each other. Under the GAL model, the distance between the densities of the 5th quantile and median labor supply also decreases with increasing number of young children, albeit at a lower rate. When estimating the 95th quantile, the proposed model is more conservative than the AL model about the labor contribution of an average wife with an increasing number of children less than 6 years old. When there are 3 children less than 6 years old in the household, the center of the posterior distribution for the 95th quantile is below zero under the GAL model, meaning that even at the top 5th percentile of labor supply, an average wife may still produce negative labor supply as she takes care of many young family members. More specifically, the posterior probability of the 95th labor supply quantile being positive is 0.19 under the GAL model, as opposed to 0.97 under the AL model. These results demonstrate that the choice of error distribution in quantile regression can have an effect on practically important conclusions for a particular application.
6 Discussion
We have developed a Bayesian quantile regression framework with a new error distribution that has flexible skewness, mode and tail behavior. The proposed model has better performance compared with the commonly used asymmetric Laplace distribution, particularly for modeling extreme quantiles. Owing to the hierarchical structure of the new distribution, posterior inference and prediction can be readily implemented via Markov chain Monte Carlo methods.
The main motivation for this work was to develop a sufficiently flexible parametric distribution that can be used as a building block for different types of quantile regression models. The extension to quantile regression with ordinal responses is a possible direction. Expanding the model to a spatial quantile regression process, along the lines of Lum and Gelfand (2012), is another direction. Finally, current work is exploring a composite quantile regression modeling framework, built from structured mixtures of generalized asymmetric Laplace distributions, to combine information from multiple quantiles of the response distribution in inference for variable selection.
Acknowledgements
This research was supported in part by the National Science Foundation under award SES1631963.
References
 Amemiya (1984) Amemiya, T. (1984), “Tobit models: A survey,” Journal of econometrics, 24, 3–61.
 Azzalini (1985) Azzalini, A. (1985), “A class of distributions which includes the normal ones,” Scandinavian Journal of Statistics, 12, 171–178.
 Chib (1992) Chib, S. (1992), “Bayes inference in the Tobit censored regression model,” Journal of Econometrics, 51, 79–99.
 Embrechts et al. (1997) Embrechts, P., Klüppelberg, C., and Mikosch, T. (1997), Modelling extremal events, vol. 33, Springer Science & Business Media.
 Gelfand and Ghosh (1998) Gelfand, A. and Ghosh, S. (1998), “Model choice: A minimum posterior predictive loss approach,” Biometrika, 85, 1–11.
 Hanson and Johnson (2002) Hanson, T. and Johnson, W. O. (2002), “Modeling regression error with a mixture of Pólya trees.” Journal of the American Statistical Association, 97, 1020–1033.
 Harrison and Rubinfeld (1978) Harrison, D. and Rubinfeld, D. L. (1978), “Hedonic housing prices and the demand for clean air,” Journal of environmental economics and management, 5, 81–102.
 Henze (1986) Henze, N. (1986), “A probabilistic representation of the skewnormal distribution,” Scandinavian Journal of Statistics, 13, 271–275.
 Hoti and Sillanpää (2006) Hoti, F. and Sillanpää, M. (2006), “Bayesian mapping of genotype expression interactions in quantitative and qualitative traits,” Heredity, 97, 4–18.
 Koenker (2005) Koenker, R. (2005), Quantile Regression, New York: Cambridge University Press.
 Kottas and Gelfand (2001) Kottas, A. and Gelfand, A. E. (2001), “Bayesian semiparametric median regression modeling,” Journal of the American Statistical Association, 96, 1458–1468.
 Kottas and Krnjajić (2009) Kottas, A. and Krnjajić, M. (2009), “Bayesian semiparametric modelling in quantile regression,” Scandinavian Journal of Statistics, 36, 297–319.
 Kotz et al. (2001) Kotz, S., Kozubowski, T., and Podgorski, K. (2001), The Laplace Distribution and Generalizations: A Revisit With Applications to Communications, Exonomics, Engineering, and Finance, Boston: Birkhäuser.
 Kozumi and Kobayashi (2011) Kozumi, H. and Kobayashi, G. (2011), “Gibbs sampling methods for Bayesian quantile regression,” Journal of Statistical Computation and Simulation, 81, 1565–1578.
 Li et al. (2010) Li, Q., Xi, R., and Lin, N. (2010), “Bayesian regularized quantile regression,” Bayesian Analysis, 5, 533–556.
 Lum and Gelfand (2012) Lum, K. and Gelfand, A. E. (2012), ‘‘Spatial quantile multiple regression using the asymmetric Laplace process (with discussion),” Bayesian Analysis, 7, 235–276.
 Mroz (1987) Mroz, T. A. (1987), “The sensitivity of an empirical model of married womenś hours of work to economic and statistical assumptions,” Econometrica: Journal of the Econometric Society, 765–799.
 Park and Casella (2008) Park, T. and Casella, G. (2008), “The Bayesian lasso,” Journal of the American Statistical Association, 103, 681–686.
 Reich et al. (2010) Reich, B. J., Bondell, H. D., and Wang, H. J. (2010), “Flexible Bayesian quantile regression for independent and clustered data,” Biostatistics, 11, 337–352.
 Reich and Smith (2013) Reich, B. J. and Smith, L. B. (2013), “Bayesian quantile regression for censored data,” Biometrics, 69, 651–660.
 Taddy and Kottas (2010) Taddy, M. and Kottas, A. (2010), “A Bayesian nonparametric approach to inference for quantile regression,” Journal of Business and Economic Statistics, 28, 357–369.
 Tibshirani (1996) Tibshirani, R. (1996), “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society. Series B (Methodological), 267–288.
 Tokdar and Kadane (2012) Tokdar, S. T. and Kadane, J. B. (2012), “Simultaneous linear quantile regression: A semiparametric Bayesian approach,” Bayesian Analysis, 7, 51–72.
 Tsionas (2003) Tsionas, E. G. (2003), “Bayesian quantile inference,” Journal of Statistical Computation and Simulation, 73, 659–674.
 Volinsky and Raftery (2000) Volinsky, C. T. and Raftery, A. E. (2000), “Bayesian information criterion for censored survival models,” Biometrics, 56, 256–262.
 Waldmann et al. (2013) Waldmann, E., Kneib, T., Yue, Y. R., Lang, S., and Flexeder, C. (2013), “Bayesian semiparametric additive quantile regression,” Statistical Modelling, 13, 223–252.
 Walker and Mallick (1999) Walker, S. G. and Mallick, B. K. (1999), “A Bayesian semiparametric accelerated failure time model,” Biometrics, 55, 477–483.
 Wichitaksorn et al. (2014) Wichitaksorn, N., Choy, B. S. T., and Gerlach, R. (2014), “A generalized class of skew distributions and associated robust quantile regression models,” The Canadian Journal of Statistics, 42, 579–596.
 Yu and Moyeed (2001) Yu, K. and Moyeed, R. A. (2001), “Bayesian quantile regression,” Statistics & Probability Letters, 54, 437–447.
 Yu and Stander (2007) Yu, K. and Stander, J. (2007), “Bayesian analysis of a Tobit quantile regression model,” Journal of Econometrics, 137, 260–276.
 Zhu and Galbraith (2011) Zhu, D. and Galbraith, J. W. (2011), “Modeling and forecasting expected shortfall with the generalized asymmetric Studentt and asymmetric exponential power distributions,” Journal of Empirical Finance, 18, 765–778.
 Zhu and ZindeWalsh (2009) Zhu, D. and ZindeWalsh, V. (2009), “Properties and estimation of asymmetric exponential power distribution,” Journal of Econometrics, 148, 86–99.
 Zou and Yuan (2008) Zou, H. and Yuan, M. (2008), “Composite quantile regression and the oracle model selection theory,” The Annals of Statistics, 36, 1108–1126.