Bayesian quantile regression analysis for continuous data with a discrete component at zero

Abstract

In this work, we propose a Bayesian quantile regression method to response variables with mixed discrete-continuous distribution with a point mass at zero, where these observations are believed to be left censored or true zeros. We combine the information provided by the quantile regression analysis to present a more complete description of the probability of being censored given that the observed value is equal to zero, while also studying the conditional quantiles of the continuous part. We build up a Markov Chain Monte Carlo method from related models in the literature to obtain samples from the posterior distribution. We demonstrate the suitability of the model to analyse this censoring probability with a simulated example and two applications with real data. The first is a well-known dataset from the econometrics literature about women labour in Britain, and the second considers the statistical analysis of expenditures with durable goods, considering information from Brazil.

Keywords

asymmetric Laplace distribution Bayesian quantile regression Durable goods left censoring two-part model

1 Introduction

In the econometrics literature, a well-known problem is the case when there is a non-negative continuous variable with a point mass at zero. Tobin (1958) considered the problem of expenditures of durable goods, assuming that all zero observations of expenditures of a household were actually censored observations. So in order to have a better understanding of the conditional mean of this response variable given other unknowns a variable Y^★, which is only observable when its value is positive, is added in the modelling scheme. Then considering that the expenditures are normal distributed, one could use the normal cumulative distribution function (cdf) in the likelihood computation. On the other hand, Cragg (1971), still investigating the effects of explanatory variables in durable goods purchases, defined a two-part model, where one part of the model studies whether the individual makes the purchase or not, then another part tries to explain how much is spent on average conditional on some variables. It is notable that these approaches deal with a similar problem using very different assumptions. Here in this article, we try to combine these ideas, without relying solely on the conditional mean to make inference about the continuous distribution, but rather on conditional quantiles of this distribution.

Recently, there has been a surge of methods that give attention to other parts than the mean of the conditional distribution, often without trying to describe this distribution with just one family of distributions [see (Kneib (2013)) for a good discussion of such models]. One of these models is the quantile regression model, which was first proposed by Koenker and Bassett (1978), and is the one we consider in this work. This model was extended from the previous work of Wagner (1959) for the median case, which was based on the results by Harris (1950). Essentially, if one believes that the regression parameters are not fixed for the entire distribution, but rather depend on the quantile of interest, then these models are able to capture this effect. For example, these models are capable of measuring differences between central and tail estimates in the conditional distribution of the response variable, which is not usually the aim of conditional mean models. Some interesting applications of quantile regression models can be seen in Yu et al. (2003), Koenker (2005) and Elsner et al.(2008).

In the frequentist framework, the quantile regression parameters are obtained using linear programming algorithms, since the minimization problem proposed can be written as linear programming problem. For a Bayesian setting, the asymmetric Laplace distribution can be useful in obtaining posterior conditional quantile estimates. Yu and Moyeed (2001) proposed the use of this distribution in order to introduce a Bayesian quantile regression model. The authors verified in simulation studies that this assumption was helpful in approximating conditional quantiles for different probability distributions. Yue and Rue (2011) used this distribution to build additive mixed quantile regression models, where they adopted the integrated nested Laplace approximations (INLA) approach for posterior inference. Later, Sriram et al. (2013) proved that this distribution is capable of providing good estimates, if certain conditions are satisfied. Moreover, Kozumi and Kobayashi (2011) suggested a more efficient Markov Chain Monte Carlo (MCMC) method, making use of location-scale mixture representation of the asymmetric Laplace distribution, and even implementing a Tobit quantile regression when left censoring is present. This new approach facilitated other extensions in the literature on Bayesian quantile regression models. For instance, Lum and Gelfand (2012) developed the asymmetric Laplace process to produce quantile estimates when the data presents spatial correlation. Alhamzawi and Yu (2012); Alhamzawi and Yu (2013) presented variable selection methods also considering this representation. Luo et al. (2012) used this approach for longitudinal data models with random effects.

Considering left censored data, Yu and Stander\rq s (2007) was one of the first proposals to consider the asymmetric Laplace distribution in order to produce a Bayesian Tobit quantile regression model. Kozumi and Kobayashi (2012) developed a Bayesian quantile regression for censored dynamic panel data. Zhao and Lian (2015) introduced a Bayesian Tobit quantile regression for single-index models. Alhamzawi (2016) proposed a Tobit quantile regression with the elastic net penalty in a Bayesian framework. Kobayashi (2016) presented Bayesian Tobit quantile regression models with endogenous variables. Alhamzawi and Yu (2015) considered the estimation of Tobit quantile regression models using a g-prior distribution with a ridge parameter for more flexibility when dealing with censored data. And as an illustration of Tobit quantile regression models, Yue and Hong (2012) explained medical expenditures in panel survey data using this approach.

Furthermore, Santos and Bolfarine (2015), considering proportion data with inflation of zeros or ones, developed a two-part model using a Bayesian quantile regression model to explain the conditional quantiles of the continuous part. They considered the equivariance property of the quantile function to work with a transformed variable in the modelling process to match the support of the asymmetric Laplace distribution. We aim to extend their model in this article, but only concentrating in the zero inflation process, where we assume that a percentage of the zeros are left censored. By examining the continuous distribution with their conditional quantiles, we plan to convey more information about this probability of being censored given that the observation is zero.

This article is organized in the following manner. In Section 2, we show the Bayesian two-part model using quantile regression for the continuous part, with its prior and posterior setting. In Section 3, we extend the two-part model to allow that zero observations are either censored or true zeros, defining the posterior probability of being censored given that it is zero. Advancing, in Section 4, we present the suitability of the method using a simulated example that compares the censoring probabilities for zero observations. Two applications of the model are presented to illustrate the results in Section 5. We complete with our final remarks in Section 6.

2 Two-part model review

If we consider the possibility of modelling the response variable by a mixture of two distributions, a point mass distribution at zero and a continuous distribution for the positive values, we can use the two-part model introduced by Cragg (1971). Then we can write the probability density function of $Y$ as

g (y) = p I (y = 0) + (1 - p) f (y) I (y > 0),

(2.1)

where

I (a)

is the indicator function, which is equal to one when

a

is true and

p = P (Y = 0)

. For the example with durable goods, this model assumes that a person decides whether it is going to make a purchase or not in the first place and later one decides about how much it will spend. For each part of the process, the analysis can verify the importance of predictors to explain the probability of being zero and also in the outcome of the positive part. Cragg (1971) considered a probit link to study the probability

p

and a truncated normal probability to fit the values greater than zero. With these assumptions, or even other distributions to the continuous part, it is possible to make inference about the conditional mean of the positive values. But we believe that it is important to check if other parts of the conditional distribution might present different conclusions about the association between the response variable and its explanatory variables. And in order to help this goal, we select to use quantile regression models.

Santos and Bolfarine (2015) proposed the use of the asymmetric Laplace distribution for the continuous part, for proportion data, when the response variable is defined between zero and one. This distribution allows the posterior inference about the quantiles of the response variable, as proved by Sriram et al. (2013), proposed initially by Yu and Moyeed (2001) and later improved by Kozumi and Kobayashi (2011).

We consider here a location-scale representation mixture of the asymmetric Laplace distribution, which was used by Kozumi and Kobayashi (2011) to present a more efficient Gibbs sampler for quantile regression models. We can say that if $Y$ is distributed according to an asymmetric Laplace distribution with parameters ${μ, σ, τ}$ , then we have

\begin{matrix} Y_{i} | v_{i} & \sim N (μ + θ v_{i}, ψ^{2} σ v_{i}) \\ v_{i} & \sim Exp (σ), \end{matrix}

where

θ = (1 - 2 τ) / (τ (1 - τ))

and

ψ^{2} = 2 / (τ (1 - τ)) .

We can also add covariates to explain the probability of being equal to zero for each observation, $p_{i}$ , by making $η (p_{i}) = z_{i}^{'} γ$ , where $η (.)$ is a link function that could be, for instance, the inverse of a normal cdf—producing the probit model—or the inverse of a logistic cdf—producing the logistic model. The set of variables, $z_{i}$ , in this case can be either different or the same as the one used to describe the continuous density, $x_{i}$ , and $γ$ is a vector parameter used to better describe the probability $p_{i}$ given the set of variables $z_{i}$ .

In order to match the support of the asymmetric Laplace distribution with the positive values in this model, we need to transform the response variable with a non-decreasing function $h : ℝ^{+} \to ℝ$ . As discussed in Santos and Bolfarine (2015), due to the equivariance property, one is allowed to model this transformed variable and still make inference about the conditional quantile in the original scale of the response variable.

If we define the sets $J = {y_{i} : y_{i} = 0}$ and $K = {y_{i} : y_{i} > 0}$ , the augmented likelihood function for the two-part model, considering the location-scale mixture of the asymmetric Laplace distribution, can be defined as

L (β (τ), γ, σ) = \prod_{y_{i} \in J} η^{- 1} (z_{i}^{'} γ) \prod_{y_{i} \in K} (1 - η^{- 1} (z_{i}^{'} γ)) f (h (y_{i}) | v_{i}) f (v_{i}) .

(2.2)

To complete the Bayesian specification, we assume priors distributions for the parameters, with a normal distribution for

γ

and

β (τ)

, and an inverse gamma distribution for

σ

. Given these definitions, we can write the full hierarchical model as

\begin{matrix} h (Y_{i}) | v_{i} & \sim p_{i} I (y_{i} \in J) + (1 - p_{i}) N (x_{i}^{'} β (τ) + θ v_{i}, ψ^{2} σ v_{i}) I (y_{i} \in K), \\ v_{i} & \sim E (σ), \\ p_{i} & = η^{- 1} (z_{i}^{'} γ), \\ β (τ) & \sim N (b_{0}, B_{0}), \\ σ & \sim IG (n_{0}, s_{0}), \\ γ & \sim N (g_{0}, G_{0}) . \end{matrix}

The full conditional distributions for all parameters, after we combine the likelihood with the prior information, are

\begin{matrix} β (τ) ∣ y, v, γ, σ & \sim N (b_{1}, B_{1}), \\ v_{i} ∣ h (y), β (τ), γ, σ & \sim GIG (1 / 2, {\hat{δ}}_{i}, {\hat{ξ}}_{i}), \\ σ ∣ y, v, β (τ), γ & \sim IG (\tilde{n} / 2, \tilde{s} / 2), \\ π (γ ∣ y, v, β (τ), σ) & p \prod_{i \in J} η^{- 1} (z_{i}^{'} γ) \prod_{i \in K} (1 - η^{- 1} (z_{i}^{'} γ)) exp \{- \frac{1}{2} (γ - g_{0})^{'} G_{0}^{- 1} (γ - g_{0})\} . \end{matrix}

(2.3)

With the exception of

γ

, all parameters of the full conditional distributions are similar to the ones in the MCMC proposed by Kozumi and Kobayashi (2011), and for that reason are not provided here.

GIG (ν, δ, ζ)

represents a generalized inverse Gaussian distribution, for which we can use the algorithm by Dagpunar (1989) to generate values from the posterior distribution for each

v_{i}

. We use a random walk Metropolis–Hastings algorithm to sample from the posterior distribution of

γ

, as in Santos and Bolfarine (2015).

Furthermore, in all Bayesian models, an important part in the set-up is prior elicitation. Throughout the article, we decided to use only diffuse priors, as these have not presented any problems and the posterior distribution seems to be insensitive to minor changes. In this case, we consider, for example, $b_{0} = g_{0} = 0$ , $B_{0} = G_{0} = 100 I$ , $\tilde{n} = 3$ and $\tilde{s} = 0.1$ . But in the presence of prior information, we believe that a subjective prior could be used instead, and that it would not undermine the algorithm presented here.

Moreover, all other full conditional distributions are relatively easy to generate posterior samples. The instructions to use our R code are available in the online supplementary material of this article.

3 Bayesian quantile regression for continuous data with a point mass distribution at zero

Going back to the problem of studying expenditures in durable goods in a determined period, we have that a proportion of households might not have made any purchase for these kinds of items, having zero as the total of expenditures. This result was analysed using two different approaches in Tobin (1958) and Cragg (1971). Here we aim to combine these ideas, while making inference about the conditional quantiles of the continuous element of this type of problem. In the biometrics literature, Moulton and Halsey (1995) proposed an extension of Cragg's model, considering that an observed zero, or a lower detection limit, can be either from the point mass distribution or from the continuous distribution, being in the latter case a censored observation. Chai and Bailey (2008) considered a similar formulation, but changed the distribution of the continuous part and also the link function to model the probability $p_{i}$ . In the survival analysis literature, an analogous situation occurs when there is a proportion of the sample which is believed to be cured from a possible illness. As a consequence, at the end of study from all the observations that are not affected by the event of interest, some are assumed to be cured and others are considered censored observations.

In this case, we should rewrite the density in (2.1) as

g (y) = [p + (1 - p) F (0)] I (y = 0) + (1 - p) f (y) I (y > 0),

(3.1)

where

F (.)

is the cdf of the continuous part.

The latent variable $C$ , which indicates the censoring mechanism, is non-observable for all zero observations. Therefore, we need to use a data augmentation algorithm in this case. Our complete cases are ${Y_{i}, C_{i}, v_{i}}$ , from which we only observe $Y_{i}$ . We discussed in the previous section the role of latent variable $v_{i}$ and how to update it in the MCMC algorithm.

For the complete cases, define the sets $C = {y_{i} : y_{i} = 0 and c_{i} = 1}$ , $D = {y_{i} : y_{i} = 0 and c_{i} = 0}$ and $K = {y_{i} : y_{i} > 0}$ of censored observations, non-censored but with response equal to zero and observations greater than zero, respectively. Then the likelihood function for $ξ = (β (τ), γ, σ)$ , without writing the conditional parameters for $F (0)$ and $f (y_{i} | v_{i})$ for notational simplicity, can be written as

L (ξ) = \prod_{y_{i} \in D} η^{- 1} (z_{i}^{'} γ) \prod_{y_{i} \in C} (1 - η^{- 1} (z_{i}^{'} γ)) F (0) \prod_{y_{i} \in K} (1 - η^{- 1} (z_{i}^{'} γ)) f (y_{i} | v_{i}) f (v_{i}) .

It is important to note that $F (0)$ varies for every observation in $C$ , depends on the parameters of the asymmetric Laplace distribution Yu and Zhang (2005) and can be defined as

F (0; β (τ), σ) = \{\begin{matrix} τ exp \{- \frac{(1 - τ) x^{'} β (τ)}{σ}\}, & if x^{'} β (τ) ⩾ ​ 0, \\ 1 - (1 - τ) exp \{\frac{τ x^{'} β (τ)}{σ}\}, & if x^{'} β (τ) < 0 . \end{matrix}

In our quantile regression models, for the censored observations, we consider the idea proposed by Chib (1992) and adapted for quantile regression by Kozumi and Kobayashi (2011), which samples $Y^{★}$ , the assumed latent variable using

Y^{★} \sim {TN}_{(- \infty, 0]} (x^{'} β (τ) + θ v, ψ^{2} σ v),

(3.2)

where

{TN}_{[a, b]} (μ, σ^{2})

denotes a truncated normal distribution, with mean

μ

and variance

σ^{2}

, in the interval

[a, b]

If we define the variable $C$ as a censoring indicator, where $C$ is equal to one when $Y$ is censored and zero otherwise, we have that the probability of being censored given that $Y$ is zero is

P (C = 1 | Y = 0) = \frac{(1 - p) F (0)}{p + (1 - p) F (0)} .

(3.3)

To complete the specification for this censoring variable, we have that for non zero observations $P (C = 1 | Y > 0) = 0$ and $P (C = 0 | Y > 0) = 1$ . If we calculate the probability in (3.3) when $σ = 1$ and for two possible values for $x^{'} (β)$ , $- 1$ and 1, considering different values for $p$ , we obtain the curves in Figure 1. Comparing the two possible values, we have that the probability of being censored is greater when the linear predictor is equal to $- 1$ , for the same $τ$ . It is easy to show that this probability is limited between zero and $1 - p$ , and this characteristic is depicted in the plot. Also, for the same $p$ , this probability is increasing in $τ$ , which is controlled by $F (0)$ .

Figure 1:

Plots for the probability of being censored for different $τ$ ’s as a function of the probability $p = P (Y = 0)$ . (a) $x^{'} β (τ) = - 1$ , (b) $x^{'} β (τ) = 1$

However, instead of evaluating $F (0)$ for every censored observation, we can replace the censored observation for its estimate based on (3.2). By doing this, our likelihood function is similar to the likelihood for the two-part model of the previous section, and the MCMC described there can be used to update the parameters here as well. In fact, one only needs to include the posterior distribution of $C$ in (2.3), after inputing the censored observations, to complete the MCMC algorithm. And its full conditional posterior distribution is

C_{i} | y, β (τ), σ, v \sim Ber (\frac{(1 - p_{i}) F (0)}{p_{i} + (1 - p_{i}) F (0)}),

which we can easily update in the MCMC algorithm. This update depends on the current values of the parameter related to the probability

p

, but also to the parameters related to the continuous part. This interesting property grants that a certain observation, with its response value equal to zero, can have different probabilities of being censored given its conditional quantile estimate and its conditional probability of being equal to zero.

A posterior estimate of the probability of being censored for each observation can be calculated as

P (C_{i} = 1 | Y, v, β (τ), σ, γ) = \sum_{k = b + 1}^{M} \frac{C_{i}^{(k)}}{M - b}, \forall i : y_{i} \in C \cup D,

(3.4)

where

C_{i}^{(k)}

is the

k

th term of the Markov chain for the censoring indicator of the

i

th observation,

M

is the length of the chain and

b

is the length of the burn-in period.

4 Simulation study for the censoring probability

In this section, we are concerned in checking the performance of our model regarding its capability of making statements about the probability of being censored given that a certain observation has its response value equal to zero. In order to accomplish that, we replicate a study where we know which observations are censored and also which ones are not censored, between those with zero as their response value, and we compute this probability of interest for each group.

We consider a model with just two covariates and the following structure as

\begin{matrix} log (\frac{p_{i}}{1 - p_{i}}) & = γ_{0} + γ_{1} z_{i 1} + γ_{2} z_{i 2}, \\ Y_{i} & = β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + ε_{i}, \end{matrix}

where

ε_{i} \sim N (0, 0 . 5^{2})

β_{0} = - 0.5

β_{1} = 0

and

β_{2} = 1.5

. We draw

x_{ij}

j = 1, 2,

from a uniform distribution, and we set

x_{ij} = z_{ij}

, that is, we use the same covariates for both parts of the model. We set our study with a scenario where there is a distinct difference between the observations that belong to the point mass distribution and the observations from the continuous part. Therefore, to achieve this goal, we use relatively large absolute values for

γ_{1}

and

γ_{2}

, and 10 and

- 10

, respectively. By defining these values, our intention is to give greater probability

p_{i}

for observations which are true zeros. On average, 50% of the sample is classified as a true zero in the start and then another 10% is classified as zero after being censored. Our sample size in this study is 500, and we report our results for 1 000 replications of this model. For the prior hyperparameters, we assumed the

b_{0} = g_{0} = 0

and

B_{0} = G_{0} = 100 I

, where

I

is the identity matrix, and for

σ

, we considered

IG (3 / 2, 0.1 / 2)

, as in Kozumi and Kobayashi (2011). All the results were based after discarding the first 500 posterior samples and considering the next 1 500 draws from the posterior, calculating the posterior mean for each parameter in these draws.

For each simulation, we calculate the probability of being censored for all observations with $y_{i} = 0$ . Then we summarize all probabilities with the mean for the group of censored observations and also the group of non-censored observations, respectively as

ζ_{C} = \sum_{y_{i} \in C} \frac{P (C_{i} = 1 | Y_{i} = 0)}{n_{C}} and ζ_{D} = \sum_{y_{i} \in D} \frac{P (C_{i} = 1 | Y_{i} = 0)}{n_{D}},

where

n_{C}

is the number of censored observations with

y_{i} = 0

and

n_{D}

is the number of non-censored observations with

y_{i} = 0

, and

P (C_{i} = 1 | Y_{i} = 0)

is calculated according to (3.4).

Figure 2:

Density estimate for the mean posterior probabilities of being censored for censored observations and also non-censored observations, based on 1 000 replicates for each scenario using a Gaussian smoothing kernel

In Figure 2, we plot the estimated density for the 1 000 $ζ_{C}$ and $ζ_{D}$ obtained in the simulation study, considering three different $τ$ ’s: 0.25, 0.50 and 0.75. We can see that for non-censored observations, the density is mostly concentrated before 0.25, and that this concentration is even more acute for the lower quantiles. For the censored group, we see that the probabilities of being censored are definitely larger than those for the non-censored, even though for the conditional quantile 0.25 is mostly concentrated around 0.30, but varying widely. For the conditional quantile 0.75, this probability is mainly concentrated around 0.70, therefore giving, on average, greater probability to censored observations in this case. Overall, we note that probabilities of being censored increase with the conditional quantile of interest, which was already discussed in Section 3.

Besides the interest in estimation of the probability of being censored, it is important to check whether the uncertainty about whether zero observations are censored or not undermine the estimation process of other parameters. There is just one note about the length of the MCMC chains in this simulation example, as we acknowledge that these sample sizes are small for the algorithm with a Metropolis–Hastings step, but we believe that this compromise was necessary due to time constraints in order to get the results in this simulation study.

Figure 3:

Density estimate for parameters $β_{2}$ , $γ_{1}$ and $γ_{2}$ , based on 1 000 values for each parameter using a Gaussian smoothing kernel

The density of the 1 000 posterior estimates of $γ_{1}$ , $γ_{2}$ and $β_{2}$ is depicted in Figure 3. We are able to check that for all three quantiles, the estimation of these parameters is not affected by the uncertainty of the censoring process. In general, these parameters are reasonably estimated. We note, however, some interesting results about the variance of these estimates. Due to the smaller error of the censoring process for $τ = 0.75$ , since the mean probability of being censored is greater for this quantile, $β_{2}$ is better estimated, with the mean value closer to 1.5 and also with less variation around this value. On the other hand, the estimates for $γ_{2}$ show a similar performance, but in the 0.25 quantile. This is probably due to the fact that, in this quantile, the observations which are true zeros get greater estimates for $p_{i}$ . Therefore, in the iteration process the latent variable $C_{i}$ , which is responsible to discriminate between the censored and non-censored observations, will be updated as 1 in fewer opportunities. Nevertheless, this interesting feature is less highlighted in the estimates of $γ_{1}$ .

Ultimately, we do not add another distribution for the error in this study, instead of the normal distribution, as we believe that such a change would not provide more information about the effectiveness of our model. Moreover, considering a larger sample would be helpful, as we expect that estimate errors would decrease with larger samples, but since the results were already satisfactory for this sample size, we decided not to continue any further. And adding more variables to each part of the model would definitely make the estimation process more difficult, but then we understand that this is a complication with which the algorithms could deal separately.

5 Applications

We exemplify our model with two applications. First, we consider the data from Mroz (1987), which was used for illustration purposes of the Tobit quantile regression model in Kozumi and Kobayashi (2011). And second, we present data about expenditures with durable goods in Brazil between 2008 and 2009, motivated by earlier considerations about this type of data by Tobin (1958) and Cragg (1971).

5.1 Labour supply data

Analysing empirical models about female labour supply, Mroz (1987) collected data about 753 married women, aged between 30 and 60 years old. The response variable of interest here is the number of hours worked for pay during the year of 1975, measured in 100h. In the sample, which was collected from the ‘Panel Study of Income Dynamics’, there are 325 women who did not work in that year, so their response variable is equal to zero. In Kozumi and Kobayashi (2011), they used this dataset to exhibit the Tobit quantile regression model, where these zero observations are assumed to be left censored. In our model, we combine the probability of being censored with the probability of being equal to zero to provide a more comprehensive study of these women who did not work that year. For covariates, we select non-wife income ( $x_{1}$ ), years of education ( $x_{2}$ ), actual years of labour market experience ( $x_{3}$ ), wife's age ( $x_{4}$ ), the number of children 5 years old or younger ( $x_{5}$ ) and the number of children between 6 and 18 years old ( $x_{6}$ ). After standardizing all covariates, we consider the following models for the probability and the conditional quantiles, without $τ$ in the formulas to abbreviate its presentation,

\begin{matrix} log (\frac{p_{i}}{1 - p_{i}}) & = γ_{0} + γ_{1} x_{i 1} + γ_{2} x_{i 2} + γ_{3} x_{i 3} + γ_{4} x_{i 4} + γ_{5} x_{i 5} + γ_{6} x_{i 6}, \\ Q_{Y_{i}} (τ | x_{i}) & = β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + β_{3} x_{i 3} + β_{4} x_{i 4} + β_{5} x_{i 5} + β_{6} x_{i 6} . \end{matrix}

In Figure 4, we present the densities of the probabilities of being censored given that a certain observation is equal to zero for distinct $τ$ ’s. As mentioned earlier, these probabilities are $τ$ -dependent, therefore the very contrasting shapes of these densities are expected. We estimate that for $τ = 0.10$ , these probabilities are mostly concentrated under 0.25, while for $τ = 0.90$ , they are well spread between 0 and 1, with mean value equal to 0.53. These differences between quantiles in the probabilities of being censored generate some variation in the posterior credible intervals for a few $γ$ ’s, as we can see in Figure 5. For instance, for $γ_{3}$ , there is a notable increase in absolute value for greater $τ$ . Also, the effect of the number of children aged older than 6 years is estimated to be different than zero and negative just for greater quantiles, $τ = 0.8$ and $τ = 0.9$ , while the coefficient for non-wife income could be considered greater than zero just for lower quantiles. Overall, the coefficients for $x_{1}$ , $x_{4}$ and $x_{5}$ are estimated to be positive. Nevertheless, the other variables have a negative effect in the probability of the hours of work being equal to zero.

Figure 4:

Density estimates of the probability of being censored for 325 women who did not work in Mroz (1987), for $τ = 0.1, 0.3, 0.5, 0.7 and 0.9$

Figure 5:

Posterior mean and 95% credible interval for $γ_{i} (τ)$ , $i = 1, \dots, 6$

Besides that, we are able to compare the variation for each variable in the model according to the probability of being censored, which is a new factor that we add with our model for this type of analysis. For example, taking into account the variables income which is not due to the wife and years of experience, there are interesting results when we compare groups with different values for this probability for two values of $τ$ .

Figure 6:

Density of the variable nwifeinc separated in two groups according to their relative censoring probability in comparison with the mean probability in a given quantile, for $τ = 0.1 and 0.5$ . Solid lines are for the group above the mean probability of being censored, while dashed lines are for the group below

In Figure 6, we oppose the distribution of non-wife's income for women who have probability of being censored below and above the average for a specific quantile. For $τ = 0.10$ , the mean probability is equal to 0.07 and for $τ = 0.50$ is 0.21. While there is not a conceivable difference in distribution in the quantile 0.10 between those two groups, for $τ = 0.50$ there is a noticeable change in the distribution for women with lower than average probability of being censored in comparison with women with higher probability. Those women more inclined to be classified as true zeros, with lower than average probability of being censored, have a distribution of income more fat in the upper tail. We should note that this income, which is not due to the wife, could be seen as a greater incentive to work zero hours.

Figure 7:

Density of the variable educ separated in two groups according to their relative censoring probability in comparison with the mean probability in a given quantile, for $τ = 0.5 and 0.9$ . Solid lines are for the group below the mean probability of being censored, while dashed lines are for the group above

Moreover, an analogous result is obtained with the variable years of experience and it is depicted in Figure 7, but now considering quantiles 0.5 and 0.9. Considering the probabilities of being censored for $τ = 0.50$ , there is not a distinguishable difference between the distribution of years of experience from the group with probabilities below the average probability of being censored in this quantile in comparison with women above the mean probability. Additionally, for $τ = 0.9$ , where the mean probability of being censored is 0.53, the group of women with probability below this is less experienced in terms of years of experience in the labour market, since their distribution is concentrated under 20 years of experience. The distribution for the other group of women peaks around 10 years, but it has a sizable part over 20 years of experience.

5.2 Durable goods expenditures in Brazil

Figure 8:

Distribution of expenditures with durable goods, with a point mass at zero, in reais

We present here a second illustration of our model using data about household expenditures in Brazil, from the ‘Consumer Expenditure Survey’ that ocurred between 2008 and 2009, which is a national survey that interviewed more than 50 000 households around the country and it is available, in Portuguese, at http://www.ibge.gov.br/home/estatistica/populacao/condicaodevida/pof/2008_2009_encaa/microdados.shtm. Due to computational limitations, we select the sample from a specific state, Maranhão, to preserve to some extent the complex sampling scheme and also because this state had some similarities to the whole country data. After making this selection, we have 2 240 observations, from which 1 062 had zero expenditure with durable goods in the period. In Figure 8, it is possible to observe the distribution of these expenditures with a mass point at zero.

We include as covariates, gender ( $x_{1}$ : 0 = male, 1 = female), race ( $x_{2}$ 0: white, 1: non-white), age in years ( $x_{3}$ ), years of education ( $x_{4}$ ) and indicator variable if the individual has credit card ( $x_{5}$ : 0 = yes, 1 = no). Similar to the previous application, we again use a logistic model to analyse the probability of being equal to zero, where we add all variables and an intercept term. Also, all variables are included to explain the continuous part of the model. We decided to model a transformed response variable, namely $\sqrt{Y}$ , instead of $Y$ due to information gain we achieve with this transformation, as before using this transformation, in the analysis for greater quantiles, all zero observations were being considered censored. Although this characteristic is not completely lost, as we see in the results presented here, it is vastly improved after using the square root of the expenditures with durable goods. Using a trimmed notation, removing the $τ$ ’s in the formulas, the following models are considered

\begin{matrix} log (\frac{p_{i}}{1 - p_{i}}) & = γ_{0} + γ_{1} x_{i 1} + γ_{2} x_{i 2} + γ_{3} x_{i 3} + γ_{4} x_{i 4} + γ_{5} x_{i 5}, \end{matrix}

(5.1)

\begin{matrix} Q_{\sqrt{Y_{i}}} (τ | x_{i}) & = β_{0} + β_{1} x_{i 1} + β_{2} x_{i 2} + β_{3} x_{i 3} + β_{4} x_{i 4} + β_{5} x_{i 5} . \end{matrix}

(5.2)

We base our conclusions after running the MCMC obtaining 50 000 samples from the posterior distribution of the parameters of interest, from which we discarded the first 10 000 observations for burn-in purposes and later considered every 40th draw. As our estimator, we calculated the posterior mean for each parameter. Using these estimators, we see in Table 1, for the model of the probability $p_{i}$ , that the only variables with credible intervals that do not contain zero are the indicators for credit card and years of education, where the former has a positive effect and the latter a negative effect in the probability of the expenditures being equal to zero, when $τ = 0.5$ . In comparison, we estimate that the odds of having zero expenditures in a given period for a person who does not have a credit card is 2.2 times the odds of a person with a credit card.

Table 1:

Mean posterior estimates and 90% credible intervals for $γ$ in model 5.1, for τ = 0.50

Variable	Estimate	Credible interval
Intercept	$-$ 2.14	[ $-$ 2.88; $-$ 1.49]
Gender	0.02	[ $-$ 0.16; 0.18]
Race	0.16	[ $-$ 0.05; 0.36]
Age	0.08	[ $-$ 0.02; 0.18]
Education	$-$ 0.12	[ $-$ 0.23; $-$ 0.01]
Credit card	0.79	[0.50; 1.09]

Figure 9:

Mean posterior estimates and 90% credible intervals for $β (τ)$

The posterior mean estimates for the continuous part are shown in Figure 9. Using quantile regression models, we are able to check if some variables have positive or negative effects in just some parts of the conditional distribution of the response variable. In this example, we find that difference between expenditures of men and women is different from zero only in the upper tail, or for $τ = 0.8$ and $τ = 0.9$ , with a negative coefficient in this case. Similarly, the estimates for years of education and race are only different from zero for $τ > 0.3$ , where education has a positive effect in the expenditure, and race is estimated to have a negative effect.

Analogous to the previous application, we can also compare the probability of being censored for different predictor variables. For $τ = 0.5$ , in Figure 10(a) we see that the probability of being censored for people with credit cards is definitely greater than for people without one. This is, in fact, in concordance with the estimates in Table 1. Furthermore, we are able to check this information for other quantiles as well. As we notice in Figure 10(b), where we plot the probabilities of being censored for all observations with response value equal to zero, varying $τ$ from 0.1 to 0.7, we recognize that, in general, the probabilities for the group with a credit card has a much faster increase when we move to greater quantiles.

We do not show the probabilities for $τ = 0.8$ and $τ = 0.9$ because when adding these values to the plot, it was more difficult to make the comparison between those groups, since for these quantiles, most observations are considered censored when their response is zero. This happens due to weight $F (0)$ has in the calculus in the probability of being censored. This factor in the probability depends on the scale of the data and that is the reason we used the square root transformation before starting the modelling process. We remember that in the previous application, there was no need for a transformation. We should note that for quantile regression models, one can use the equivariance to monotone transformations property for the quantile function, in order to get the estimated quantiles in the original scale of the data. But for the two-part model, one problem with this particular setting, where most observations are considered censored for some quantiles, is that the posterior learning for the parameters $γ$ is compromised when this happens. In this example, we tested other transformations, but we left the square root just to emphasize that this feature should be carefully analysed in each application.

Figure 10:

Comparisons for the probabilities of being censored given the indicator variable credit card: (a) density estimates for $τ = 0.5$ , yes = solid line, no = dashed line and (b) profiles of the probabilities for $τ = 0.1, 0.2, \dots, 0.7$

6 Final remarks

In this article, we extended the Bayesian Tobit quantile regression model to include the probability of a response variable equal to zero being censored, instead of considering them censored observations from the beginning, using the asymmetric Laplace distribution in the likelihood to analyse the conditional quantiles of the continuous part of the model. We used a two-part model to study the probability of $Y = 0$ , denoting a point mass distribution at zero, while also building a linear predictor for the conditional quantiles of the continuous distributions for the response variable, meaning all observations where $Y > 0$ and when $Y$ is considered to be censored.

Although this extension is based on the results obtained initially by Santos and Bolfarine (2015), we added an important feature in the modelling framework, which is the analysis of zero observations where the probability of being censored is taken into account for all zero observations. In the cited paper, the authors are only interested in measuring the difference between zero and non zero observations, while also studying the conditional quantiles for all observations in the continuous part of the model. Here, we showed how one can use the evidence in these two parts of the model to weight if a zero observation should be deemed censored given a fixed quantile or part of a point mass distribution at zero.

We showed how the probability of being censored, in our model, is dependent on the quantile of interest, providing more information about whether one observation should be considered indeed censored or not, for instance, comparing the profiles of probabilities from different observations and checking their variation for smaller or greater $τ$ ’s. In comparison with the method proposed by Moulton and Halsey (1995), where this probability is not $τ$ dependent, we are able to vary this value according to the information in the data. Here, $p$ , which is the probability of $Y = 0$ , should be interpreted as the proportion of the data assumed to be generated from a point mass distribution when we analyse a certain quantile of interest.

We demonstrated our method in a simulated dataset, where we can control which observations are censored and which are non-censored, and we obtained contrasting estimates of the probabilities of being censored for these two types of observations when $Y = 0$ . We illustrated our findings with a well-known dataset in the econometrics literature about female labour supply, where we showed how this probability of being censored given some covariates affects the model for different $τ$ ’s. We also exemplified our model considering the problem of analysing expenditures with durable goods in Brazil. Again, we found interesting results for the probabilities of being censored given the indicator variable of credit card.

In quantile regression models, it is always important to discuss the possibility of crossing quantiles, that is, if the expected ordering of $Q_{y} (τ_{1} | x) < Q_{y} (τ_{2} | x)$ , when $τ_{1} < τ_{2}$ , is observed. In fact, a few proposals in the literature, such as Tokdar and Kadane (2011), Yang and Tokdar (2016) and Rodrigues and Fan (2016)) design their conditional quantile regression model respecting this constraint. In our model, because we fit several quantile regression models separately and we are adding a probability of being censored in the process of estimating these conditional quantiles, it is necessary to examine whether this two-part model setting amplifies the probability of this issue. Using the inputing mechanism described in Section 3 for the latent variable $C$ , the two-part proposal does not increase the likelihood of crossing quantiles more than the censoring process. This is due to the fact that probability of being censored as presented in Section 3 observes the monotonicity constraint of the quantile function.

It is important to mention that our model could also be used in the survival analysis framework, when there is an assumption of a cure fraction in the study. Minor modifications would be necessary just to change the left censoring at zero for right censoring in this case. For future research, we are currently developing variable selection methods that try to share the information across both parts of the model and also when the variable selection method is quantile dependent.

Footnotes

Acknowledgments

This research was supported by the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under Grants 2012/20267-9 and 2013/04419-6.

References

Alhamzawi

(2016) Bayesian elastic net Tobit quantile regression. Communications in Statistics: Simulation and Computation , 45, 2409–27.

Alhamzawi

(2012) Variable selection in quantile regression via Gibbs sampling. Journal of Applied Statistics , 39, 799–813.

Alhamzawi

(2013) Conjugate priors and variable selection for Bayesian quantile regression. Computational Statistics and Data Analysis , 64, 209–19.

Alhamzawi

(2015) Bayesian Tobit quantile regression using g-prior distribution with ridge parameter. Journal of Statistical Computation and Simulation , 85, 2903–18.

Chai

Bailey

(2008) Use of log-skew-normal distribution in analysis of continuous data with a discrete component at zero. Statistics in Medicine , 27, 3643–55.

Chib

(1992) Bayes inference in the Tobit censored regression model. Journal of Econometrics , 51, 79–99.

Cragg

(1971) Some statistical models for limited dependent variables with application to the demand for durable goods. Econometrica , 39, 829–44.

Dagpunar

(1989) An easily implemented generalised inverse Gaussian generator. Communications in Statistics: Simulation and Computation , 18, 703–10.

Elsner

Kossin

Jagger

(2008) The increasing intensity of the strongest tropical cyclones. Nature , 455, 92–5.

10.

Harris

(1950) Regression using minimum absolute deviations. The American Statis- tician , 4, 14–5.

11.

Kneib

(2013) Beyond mean regression. Statistical Modelling , 13, 275–303.

12.

Kobayashi

(2016). Bayesian endogenous Tobit quantile regression. Bayesian Analysis . \doi 10.1214/16-BA996

13.

Koenker

(2005) Quantile regression . Cambridge University Press.

14.

Koenker

Bassett

(1978) Regression quantiles. Econometrica , 46, 33–50.

15.

Kozumi

Kobayashi

(2011) Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation , 81, 1565–78.

16.

Kozumi

and Kobayashi

(2012) Bayesian analysis of quantile regression for censored dynamic panel data. Computational Statistics , 27, 359–80.

17.

Lum

and Gelfand

(2012) Spatial quantile multiple regression using the asymmetric Laplace process. Bayesian Analysis , 7, 1–24.

18.

Luo

Lian

Tian

(2012) Bayesian quantile regression for longitudinal data models. Journal of Statistical Computation and Simulation , 82, 1635–49.

19.

Moulton

Halsey

(1995) A mixture model with detection limits for regression analyses of antibody response to vaccine. Biometrics , 51, 1570–78.

20.

Mroz

(1987) The sensitivity of an empirical model of married women's hours of work to economic and statistical assumptions. Econometrica , 55, 765–99.

21.

Rodrigues

Fan

(2016) Regression adjustment for noncrossing Bayesian quantile regression. Journal of Computational and Graphical Statistics . \doi 10.1080/ 10618600.2016.1172016

22.

Santos

Bolfarine

(2015) Bayesian analysis for zero-or-one inflated proportion data using quantile regression. Journal of Statistical Computation and Simulation , 85, 3579–93.

23.

Sriram

Ramamoorthi

Ghosh

(2013) Posterior consistency of Bayesian quantile regression based on the misspecified asy- mmetric Laplace density. Bayesian Analysis , 8, 479–504.

24.

Tobin

(1958) Estimation for relationship for limited dependent variables. Econometrica , 26, 24–36.

25.

Tokdar

Kadane

(2011) Simultaneous linear quantile regression: A semiparametric Bayesian approach. Bayesian Analysis , 6, 1–22.

26.

Wagner

(1959) Linear programming techniques for regression analysis. Journal of the American Statistical Association , 54, 206–12.

27.

Yang

Tokdar

(2016) Joint estimation of quantile planes over arbitrary predictor spaces. Journal of the American Statistical Association . \doi 10.1080/01621459.2016. 1192545

28.

Stander

(2003) Quantile regression: Applications and current research areas. Journal of the Royal Statistical Society, Series D (The Statistician) , 53, 331–50.

29.

Moyeed

(2001) Bayesian quantile regression. Statistics and Probability Letters , 54, 437–47.

30.

Stander

(2007) Bayesian analysis of a Tobit quantile regression model. Journal of Econometrics , 137, 260–76.

31.

Zhang

(2005) A three-parameter asymmetric Laplace distribution and its extension. Communications in Statistics: Theory and Methods , 34, 1867–79.

32.

Yue

Hong

(2012) Bayesian Tobit quantile regression model for medical expenditure panel survey data. Statistical Modelling , 12, 323–46.

33.

Yue

Rue

(2011) Bayesian inference for additive mixed quantile regression models. Computational Statistics and Data Analysis , 55, 84–96.

34.

Zhao

Lian

(2015) Bayesian Tobit quantile regression with single-index models. Journal of Statistical Computation and Simulation , 85, 1247–63.

Bayesian quantile regression analysis for continuous data with a discrete component at zero

Abstract

Keywords

1 Introduction

2 Two-part model review

Plots for the probability of being censored for different τ ’s as a function of the probability p = P ( Y = 0 ) . (a) x ′ β ( τ ) = − 1 , (b) x ′ β ( τ ) = 1

Figure 2:

Density estimate for the mean posterior probabilities of being censored for censored observations and also non-censored observations, based on 1 000 replicates for each scenario using a Gaussian smoothing kernel

Density estimate for parameters β 2 , γ 1 and γ 2 , based on 1 000 values for each parameter using a Gaussian smoothing kernel

5.1 Labour supply data

Figure 4:

Density estimates of the probability of being censored for 325 women who did not work in Mroz (1987), for τ = 0.1 , 0.3 , 0.5 , 0.7 and 0.9

Posterior mean and 95% credible interval for γ i ( τ ) , i = 1 , … , 6

Figure 8:

Distribution of expenditures with durable goods, with a point mass at zero, in reais

Mean posterior estimates and 90% credible intervals for γ in model 5.1, for τ = 0.50

Mean posterior estimates and 90% credible intervals for β ( τ )

Comparisons for the probabilities of being censored given the indicator variable credit card: (a) density estimates for τ = 0.5 , yes = solid line, no = dashed line and (b) profiles of the probabilities for τ = 0.1 , 0.2 , … , 0.7

Footnotes

Acknowledgments

References

Plots for the probability of being censored for different $τ$ ’s as a function of the probability $p = P (Y = 0)$ . (a) $x^{'} β (τ) = - 1$ , (b) $x^{'} β (τ) = 1$

Density estimate for parameters $β_{2}$ , $γ_{1}$ and $γ_{2}$ , based on 1 000 values for each parameter using a Gaussian smoothing kernel

Density estimates of the probability of being censored for 325 women who did not work in Mroz (1987), for $τ = 0.1, 0.3, 0.5, 0.7 and 0.9$

Posterior mean and 95% credible interval for $γ_{i} (τ)$ , $i = 1, \dots, 6$

Mean posterior estimates and 90% credible intervals for $γ$ in model 5.1, for τ = 0.50

Mean posterior estimates and 90% credible intervals for $β (τ)$

Comparisons for the probabilities of being censored given the indicator variable credit card: (a) density estimates for $τ = 0.5$ , yes = solid line, no = dashed line and (b) profiles of the probabilities for $τ = 0.1, 0.2, \dots, 0.7$