{'fraction_non_alphanumeric': 0.04962391325583667, 'fraction_numerical': 0.01709485200742405, 'mean_word_length': 4.405808580858086, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
bqror: An R package for Bayesian Quantile Regression in Ordinal Models

This article describes an R package bqror that estimates Bayesian quantile regression for ordinal models introduced in Rahman (2016). The paper classifies ordinal models into two types and offers computationally efficient, yet simple, Markov chain Monte Carlo (MCMC) algorithms for estimating ordinal quantile regression. The generic ordinal model with 3 or more outcomes (labeled ORI model) is estimated by a combination of Gibbs sampling and Metropolis-Hastings algorithm. Whereas an ordinal model with exactly 3 outcomes (labeled ORII model) is estimated using Gibbs sampling only. In line with the Bayesian literature, we suggest using marginal likelihood for comparing alternative quantile regression models and explain how to compute the same. The models and their estimation procedures are illustrated via multiple simulation studies and implemented in two applications. The article also describes several other functions contained within the bqror package, which are necessary for estimation, inference, and assessing model fit.


Introduction
Quantile regression defines the conditional quantiles of a continuous dependent variable as a function of the covariates without assuming any error distribution (Koenker and Bassett, 1978). The method is robust and offers several advantages over least squares regression such as desirable equivariance properties, invariance to monotone transformation of the dependent variable, and robustness against outliers (Koenker, 2005;Davino et al., 2014;Furno and Vistocco, 2018). The Bayesian approach to quantile regression was proposed by Yu and Moyeed (2001), where they constructed a working likelihood based on an asymmetric Laplace (AL) distribution (Yu and Zhang, 2005). The approach was further refined and simplified in Kozumi and Kobayashi (2011). Most statistical softwares now have built-in procedures for estimating quantile regression with continuous responses. In R software and the within the Classical approach, the package quantreg provides functions for estimating the conditional quantiles of linear and non-linear parametric and non-parametric models (Koenker, 2021). The R package lqmm considers estimation and inference of quantile mixed-effect models (Geraci and Bottai, 2014;Geraci, 2014). Besides, the R package SIQR by Zu and Yu (2021) implements single-index quantile regression models via an iterative local linear approach. In comparison, the R package bayesQR adopts the Bayesian approach for estimating quantile regression with cross section data.
Quantile regression with discrete outcomes is more complex because quantiles of discrete data cannot be obtained through a simple inverse operation of the cumulative distribution function (cdf ). Besides, discrete outcome (binary and ordinal) modeling requires location and scale restrictions to uniquely identify the parameters (see Section 2 for details). Kordas (2006) estimated quantile regression with binary outcomes using simulated annealing, while Benoit and Van den Poel (2010) proposed Bayesian binary quantile regression; the estimation procedure of the latter is available in the bayesQR package of R software (Benoit and Van den Poel, 2017). Some recent works on Bayesian quantile regression with binary longitudinal (panel) outcomes are Rahman and Vossmeyer (2019) and Bresson et al. (2021). Extending the quantile framework to ordinal outcomes is more intricate than binary quantile regression due to the difficulty in satisfying the ordering of cut-points while sampling. Rahman (2016) introduced Bayesian quantile analysis of ordinal data and proposed two efficient MCMC algorithms. Since Rahman (2016), ordinal quantile regression has attracted some attention which includes Alhamzawi (2016), Alhamzawi and Ali (2018), Ghasemzadeh et al. (2018), Rahman and Karnawat (2019), Ghasemzadeh et al. (2020), and Tian et al. (2021).
Ordinal outcomes frequently occur in a wide class of applications in economics, finance, marketing, and the social sciences. Here, ordinal regression (e.g. ordinal probit, ordinal logit) is an important tool for modeling, analysis, and inference. Given the prevalence of ordinal models in applications and the recent theoretical developments surrounding ordinal quantile regression, an estimation package is essential so that applied econometricians and statisticians can benefit from a more comprehensive data analysis. The current paper fills this gap and describes the implementation of the bqror package (version 1.4.0) for estimating Bayesian ordinal quantile regression. The package offers two MCMC algorithms to exploit the gains in computation. Ordinal model with 3 or more outcomes and a certain identification scheme is referred to as OR I model 1 and estimated through a combination of 1 In contrast to Rahman (2016), our definition of OR I model includes an ordinal model with exactly 3 outcomes.
Gibbs sampling (Geman and Geman, 1984;Casella and George, 1992) and Metropolis-Hastings (MH) algorithm (Chib and Greenberg, 1995). The method is implemented in the function quantregOR1. For ordinal models with exactly 3 outcomes, an alternative parameter identification scheme can be employed for a simpler and faster algorithm based on Gibbs sampling only. As in Rahman (2016), this is labeled as OR I I model and implemented in the function quantregOR2. Both functions, quantregOR1 and quantregOR2, report the posterior mean, posterior standard deviation, and 95% posterior credible interval of the model parameters. To compare alternative quantile regression models, we recommend the use of marginal likelihood over the deviance information criterion (DIC) as reported in Rahman (2016). This is because the "Bayesian approach" to compare models is via the marginal likelihood (Chib, 1995;Chib and Jeliazkov, 2001). So, the bqror package also provides the marginal likelihood with technical details for computation described in this paper. Additionally, the package includes functions to calculate the covariate effects and example codes to produce trace plots of MCMC draws . Lastly, this paper utilizes the bqror package to demonstrate the estimation of quantile ordinal models on simulated data and real-life applications from Rahman (2016).

Quantile regression in ordinal models
Ordinal outcomes are common in a wide class of applications in economics, finance, marketing, social sciences, statistics in medicine, and transportation. In a typical study, the observed outcomes are ordered and categorical; so for the purpose of analysis scores/numbers are assigned to each outcome. For example, in a study on public opinion about offshore drilling (Mukherjee and Rahman, 2016), responses may be recorded as follows: 1 for 'strongly oppose', 2 for 'somewhat oppose', 3 for 'somewhat support', and 4 for 'strongly support'. The numbers have an ordinal meaning but have no cardinal interpretation. We cannot interpret a score of 2 as twice the support compared to a score of 1, or the difference in support between 2 and 3 is the same as that between 3 and 4. With ordinal outcomes, the primary modeling objective is to express the probability of outcomes as a function of covariates. Ordinal models that have been extensively studied and employed in applications include the ordinal probit and ordinal logit models (Johnson and Albert, 2000;Greene and Hensher, 2010), but they only give information about the average probability of outcomes conditional on the covariates.
Quantile regression with ordinal outcomes can be estimated using the monotone equivariance property and provides information on the probability of outcomes at different quantiles. In the spirit of Albert and Chib (1993), the ordinal quantile regression model can be presented in terms of an underlying latent variable z i as follows: where x i is a 1 × k vector of covariates, β p is a k × 1 vector of unknown parameters at the p-th quantile, i follows an AL distribution i.e., i ∼ AL(0, 1, p), and n denotes the number of observations. Although z i is latent (or unobserved), it is related to the observed discrete response y i through the following relationship, where γ p = (γ p,0 = −∞, γ p,1 , . . . , γ p,J−1 , γ p,J = ∞) is the cut-point vector and J denotes the number of outcomes or categories. Typically, the cut-point γ p,1 is fixed at 0 to anchor the location of the distribution required for parameter identification (Jeliazkov and Rahman, 2012). Given the observed data y = (y 1 , · · · , y n ) , the likelihood function for the ordinal quantile model can be written as, where Θ p = (β p , γ p ), F AL (·) denotes the cdf of an AL distribution and I(y i = j) is an indicator function, which equals 1 if y i = j and 0 otherwise. While Bayesian quantile regression formulates a working likelihood by assuming an AL error distribution, working directly with the AL distribution is not convenient for MCMC sampling. Therefore, the latent formulation of the ordinal quantile model (equation 1) is expressed in the normal-exponential mixture form (Kozumi and Kobayashi, 2011) as follows, This was done to address a comment from an anonymous referee that a 3 outcome ordinal model can also be estimated using a combination of Gibbs and MH algorithms. However, we recommend using the Gibbs sampling algorithm for estimating an ordinal model with exactly 3 outcomes since it is faster and more efficient.
and E denotes normal and exponential distributions, respectively; θ = (1 − 2p)/[p(1 − p)] and τ = 2/[p(1 − p)]. Based on this formulation, we can write the conditional distribution of the latent variable as z i |β p , w i ∼ N(x i β p + θw i , τ 2 w i ) for i = 1, . . . , n. This allows access to the properties of normal distribution which helps in constructing efficient MCMC algorithms.

OR I model
The term "OR I model" is ascribed to an ordinal model in which the number of outcomes (J) is equal to or greater than 3 (see Footnote 1), location restriction is imposed by setting γ p,1 = 0, and scale restriction is achieved via constant variance (See Rahman (2016); for a given value of p, variance of a standard AL distribution is constant). The location and scale restrictions are necessary to uniquely identify the parameters (see Jeliazkov et al., 2008;Jeliazkov and Rahman, 2012, for further details and a pictorial representation).
In estimation of OR I model, we need to preserve the ordering of cut-points (γ p,0 = −∞ < γ p,1 < γ p,2 < . . . < γ p,J−1 < γ p,J = ∞) during the sampling process. This is achieved by using a monotone transformation from a compact set to the real line. Many such transformations are available (e.g., log-ratios of category bin widths, arctan, arcsin), but the bqror package utilizes the logarithmic transformation (Albert and Chib, 2001;Jeliazkov et al., 2008;Rahman, 2016), This cut-points (γ p,1 , γ p,2 , · · · , γ p,J−1 ) can be obtained from a one-to-one mapping to (δ p,2 , · · · , δ p,J−1 ). With all the modeling ingredients in place, we employ the Bayes' theorem and express the joint posterior distribution as proportional to the product of the likelihood and prior distributions. As in Rahman (2016), we employ independent normal priors: β p ∼ N(β p0 , B p0 ), δ p ∼ N(δ p0 , D p0 ) in the bqror package. The augmented joint posterior distribution for the OR I model can thus be written as, where in the likelihood function of the second line, we use the fact that the observed y i is independent of (β p , w) given (z i , δ p ). This follows from equation (2) which shows that y i given (z i , δ p ) is determined with probability 1. In the third line, we specify the conditional distribution of the latent variable and the prior distribution on the parameters.
The conditional posterior distributions are derived from the augmented joint posterior distribution (i.e. equation 6), and the parameters are sampled as per Algorithm 1. This algorithm is implemented in the bqror package. The parameter β p is sampled from an updated multivariate normal distribution and the latent weight w is sampled element-wise from a generalized inverse Gaussian (GIG) distribution. The cut-point vector δ p is sampled marginally of (z, w) using a random-walk MH algorithm. Lastly, the latent variable z is sampled element-wise from a truncated normal (TN) distribution.
Algorithm 1: Sampling in OR I model.
• Sample w i |β p , z i ∼ GIG (0.5,λ i ,η), for i = 1, · · · , n, where, • Sample δ p |y, β marginally of w (latent weight) and z (latent data), by generating δ p using a random-walk chain δ p = δ p + u, where u ∼ N(0 J−2 , ι 2D ), ι is a tuning parameter andD denotes negative inverse Hessian, obtained by maximizing the log-likelihood with respect to δ p . Given the current value of δ p and the proposed draw δ p , return δ p with probability, otherwise repeat the old value δ p . The variance of u may be tuned as needed for appropriate step size and acceptance rate.

OR II model
The term "OR II model" is used for an ordinal model with exactly 3 outcomes (i.e., J = 3) where both location and scale restrictions are imposed by fixing cut-points. Since there are only 2 cut-points and both are fixed at some values, the scale of the distribution needs to be free. Therefore, a scale parameter σ p is introduced and the quantile ordinal model rewritten as follows: where σ p i ∼ AL(0, σ p , p), (γ 1 , γ 2 ) are fixed, and recall γ 0 = −∞ and γ 3 = ∞ . In this formulation, the conditional mean of z i is dependent on σ p which is problematic for Gibbs sampling. So, we define a new variable ν i = σ p w i and rewrite the model in terms of ν i . In this representation, z i |β p , σ p , ν i ∼ N(x i β p + θν i , τ 2 σ p ν i ), the conditional mean is free of σ p and the model is conducive to Gibbs sampling.
The next step is to specify the prior distributions required for Bayesian inference. We follow Rahman (2016) and assume β p ∼ N(β p0 , B p0 ), σ p ∼ IG(n 0 /2, d 0 /2), and ν i ∼ E (σ p ); where IG stands for inverse-gamma distribution. These are the default prior distributions in the bqror package. Employing the Bayes' theorem, the augmented joint posterior distribution can be expressed as, where the derivations in each step are analogous to those for the OR I model.
The augmented joint posterior distribution (i.e., equation 8) can be utilized to derive the conditional posterior distributions and the parameters are sampled as presented in Algorithm 2. This involves sampling β p from an updated multivariate normal distribution and sampling σ p from an updated IG distribution. The latent weight ν is sampled element-wise from a GIG distribution and similarly, the latent variable z is sampled element-wise from a truncated normal distribution.
Algorithm 2: Sampling in OR II model.

Marginal likelihood
The article by Rahman (2016) employed DIC (Spiegelhalter et al., 2002;Gelman et al., 2013) for model comparison. However, in the Bayesian framework, alternative models are typically compared using marginal likelihood or Bayes factor (Poirier, 1995;Greenberg, 2012). Therefore, we prefer using marginal likelihood (or Bayes factor) for comparing two or more regression models at a given quantile.
Consider a model M s with parameter vector Θ s . Let f (y|M s , Θ s ) be its sampling density, and π(Θ s |M s ) be the prior distribution; where s = 1, . . . , S. Then, the marginal likelihood for the model M s is given by the expression, m(y) = f (y|Θ s )π(Θ s ) dΘ s , where for notational simplicity we have suppressed the conditioning on M s . The Bayes factor is the ratio of marginal likelihoods. So, for any two models M q versus M r , the Bayes factor, B qr = , can be easily computed once we have the marginal likelihoods. Chib (1995) and later Chib and Jeliazkov (2001) showed that a simple and stable estimate of marginal likelihood can be obtained from the MCMC outputs. The approach is based on the recognition that the marginal likelihood can be written as the product of likelihood function and prior density over the posterior density. So, the marginal likelihood m(y|M s ) for model M s is expressed as, Chib (1995) refers to equation (9) as the basic marginal likelihood identity since it holds for all values in the parameter space, but typically computed at a high density point (e.g., mean, mode) denoted Θ * to minimize estimation variability. The likelihood ordinate f (y|M s , Θ * ) is directly available from the model and the prior density π(Θ * |M s ) is assumed by the researcher. The novel part is the computation of posterior ordinate π(Θ * |y), which is estimated using the MCMC outputs obtained from conditional posterior densities. Since marginal likelihood is often a large number, it is convenient to express it on the logarithmic scale. An estimate of the logarithm of marginal likelihood is given by, where, once again, we have dropped the conditioning on M s . The following two subsections explain the computation of the marginal likelihood for the OR I and OR I I quantile regression models.

Marginal likelihood for OR I model
We know from Section 2.1 that the MCMC algorithm for estimating the OR I model is defined by the following conditional posterior densities: π(β p |z, w), π(δ p |β p , y), π(w|β p , z), and π(z|β p , δ p , w, y). The conditional posteriors for β p , w, and z have a known form, but that of δ p is not tractable and is sampled using an MH algorithm. Consequently, we adopt the approach of Chib and Jeliazkov (2001) to calculate the marginal likelihood for the OR I model.
To simplify the computational process (specifically, to keep the computation over a reasonable dimension), we estimate the marginal likelihood marginally of the latent variables (w, z). Moreover, we decompose the posterior ordinate as, where Θ * = (β * p , δ * p ) denotes a high density point. By placing the intractable posterior ordinate first, we avoid the MH step in the reduced run -the process of running an MCMC sampler with one or more parameters fixed at some value -of the MCMC sampler. We first estimate π(δ * p |y) and then the reduced conditional posterior ordinate π(β * p |δ * p , y). To obtain an estimate of π(δ * p |y), we need to express it in a computationally convenient form. The parameter δ p is sampled using an MH step, which requires specifying a proposal density. Let q(δ p , δ p |β p , w, z, y) denote the proposal density for the transition from δ p to δ p , and let, denote the probability of making the move. In the context of the model, f (y|β p , δ p ) is the likelihood given by equation (3), π(β p ) and π(δ p ) are normal prior distributions (i.e., β p ∼ N(β p0 , B p0 ) and δ p ∼ N(δ p0 , D p0 ) as specified in Section 2.1), and the proposal density q(δ p , δ p |β p , w, z, y) is normal given by f N (δ p |δ p , ι 2D ) (see Algorithm 1 in Section 2.1). There are two points to be noted about the proposal density. First, the conditioning on (β p , w, z, y) is only for generality and not necessary as illustrated by the use of a random-walk proposal density. Second, since our MCMC sampler utilizes a random-walk proposal density, the second ratio on the right hand side of equation (11) reduces to 1.
We closely follow the derivation in Chib and Jeliazkov (2001) and arrive at the following expression of the posterior ordinate, where E 1 represents expectation with respect to the distribution π(β p , δ p , w, z|y) and E 2 represents expectation with respect to the distribution π(β p , w, z|δ * p , y) × q(δ * p , δ p |β p , w, z, y). The quantities in equation (12) can be estimated using MCMC technique. To estimate the numerator, we take the from the complete MCMC run and take an average of the quantity α MH (δ p , δ * p |β p , w, z, y) q(δ p , δ * p |β p , w, z, y), where α MH (·) is given by equation (11) with δ p replaced by δ * p , and q(δ p , δ * p |β p , w, z, y) is given by the normal density f N (δ * p |δ p , ι 2D ). The estimation of the quantity in the denominator is tricky! This requires generating an additional sample (say of H iterations) from the reduced conditional densities: π(β p |w, z), π(w|β p , z), and π(z|β p , δ * p , w, y), where note that δ p is fixed at δ * p in the MCMC sampling, and thus corresponds to a reduced run. Moreover, at each iteration, we generate The resulting quadruplet of draws {β With the numerator and denominator now available, the posterior ordinate π(δ * p |y) is estimated as, The computation of the posterior ordinate π(β * p |δ * p , y) is trivial. We have the sample of H draws {w (h) , z (h) } from the reduced run, which are marginally of β p from the distribution π(w, z|δ * p , y). These draws are utilized to estimate the posterior ordinate as, Substituting the two density estimates given by equations (13) and (14) into equation (10), an estimate of the logarithm of marginal likelihood for OR I model is obtained as, where the likelihood f (y|β * p , δ * p ) (see equation 3) and prior densities are evaluated at Θ * = (β * p , δ * p ).

Marginal likelihood for OR I I model
We know from Section 2.2 that the OR I I model is estimated by Gibbs sampling and hence we follow Chib (1995) to compute the marginal likelihood. The Gibbs sampler consists of four conditional posterior densities given by π(β p |σ p , ν, z), π(σ p |β p , ν, z), π(ν|β p , σ p , z), and π(z|β p , σ p , ν, y). However, the variables (ν, z) are latent. So, we integrate them out and write the posterior ordinate as where the terms on the right hand side can be written as, and Θ * = (β * p , σ * p ) denotes a high density point, such as the mean or median. The posterior ordinate π(β * p |y) can be estimated as the ergodic average of the conditional posterior density with the posterior draws of (σ p , ν, z). Therefore, π(β * p |y) is estimated as, The term π(σ * p |β * p , y) is a reduced conditional density ordinate and can be estimated with the help of a reduced run. So, we generate an additional sample (say another G iterations) of {ν (g) , z (g) } from π(ν, z|β * p , y) by successively sampling from π(σ p |β * p , ν, z), π(ν|β * p , σ p , z), and π(z|β * p , σ p , ν, y), where note that β p is fixed at β * p in each conditional density. Next, we use the draws {ν (g) , z (g) } to compute, which is a simulation consistent estimate of π(σ * p |β * p , y). Substituting the two density estimates (equations 16 and 17) into equation (10), we obtain an estimate of the logarithm of marginal likelihood, where the likelihood function and prior densities are evaluated at Θ * = (β * p , σ * p ). Here, the likelihood function has the expression, where the cut-points γ are known and fixed for identification reasons as explained in Section 2.2.

Simulation studies
This section explains the data generating process for simulation studies, the functions offered in the R package bqror, and how to use these functions for estimation of ordinal quantile models on the simulated data.
Here, U and AL denote a uniform and an asymmetric Laplace distributions, respectively. The z values are continuous and are classified into 4 categories based on the cut-points (0, 2, 4) to generate ordinal discrete values of y, the outcome variable. We follow the above procedure to generate 3 data sets each with 500 observations (i.e., n = 500). The 3 data sets correspond to the quantile p equaling 0.25, 0.50, and 0.75, and are stored as data25j4, data50j4, and data75j4, respectively. Note that the last two letters in the name of the data (i.e., j4) denote the number of unique outcomes in the y variable.
Functions and Estimation: We now describe the major functions for Bayesian quantile estimation of OR I model, demonstrate their usage, and note the inputs and outputs of each function.

quantregOR1
The function quantregOR1 is the primary function for estimating Bayesian quantile regression in ordinal models with 3 or more outcomes (i.e., OR I model) and implements Algorithm 1. In the code snippet below, we first read the data and then do the following: define the ordinal response variable (y) and covariates (x), specify the number of covariates (k) and number of outcomes (J), and set the prior means and covariances for β p and δ p . We then call the quantregOR1 function and specify the following inputs: ordinal outcome variable (y), covariate matrix including a column of ones (xMat), prior mean (b0) and prior covariance matrix (B0) for the regression coefficients, prior mean (d0) and prior covariance matrix (D0) for the transformed cut-points, burn-in size (burn), post burn-in size (mcmc), quantile (p), and the tuning factor (tune) to adjust MH acceptance rate. The last input verbose, when set to (FALSE) TRUE will (not) print the outputs.
In the code below, we use a sufficiently diffuse normal prior β p ∼ N(0 k , 10 * I k ) where 0 k is a zero vector and I k is an identity matrix of dimension k × k. The prior distribution can be further diffused (i.e., less informative) by increasing 10 to say 100. Besides, the prior variance for δ p needs to be small, such as 0.25 * I J−2 since the distribution is on the transformed cut-points. If there is a need for prior elicitation, one approach is to estimate the model on a training sample and then use the results to design prior distributions (See Greenberg, 2012, for examples).

covEffectOR1
The function covEffectOR1 computes the average covariate effect, marginally of the parameters and the remaining covariates, for different outcomes of OR I model at specified quantiles. While a demonstration of this function is best understood in a real-life study and is presented in the application section, here we present the mechanics behind the computation of average covariate effect.
Suppose, we want to compute the average covariate effect when the l-th covariate {x i,l } is set to the values a and b, denoted as {x a i,l } and {x b i,l }, respectively. We split the covariate and parameter vectors as follows: x a i = (x a i,l , x i,−l ), x b i = (x b i,l , x i,−l ), and β = (β p,l , β p,−l ), where −l in the subscript denotes all covariates (parameters) except the l-th covariate (parameter). We are interested in the distribution of the difference {Pr(y i = j|x b i,l ) − Pr(y i = j|x a i,l )} for 1 ≤ j ≤ J, marginalized over {x i,−l } and (β p , δ p ), given the data y = (y 1 , · · · , y n ) . We marginalize the covariates using their empirical distribution and integrate the parameters based on the posterior distributions.
To obtain draws from the distribution {Pr(y i = j|x b i,l ) − Pr(y i = j|x a i,l )}, we use the method of composition (See Chib and Jeliazkov, 2006;Jeliazkov et al., 2008;Rahman and Vossmeyer, 2019;Bresson et al., 2021, for additional details). In this process, we randomly select an individual, extract the corresponding sequence of covariate values, draw a value (β p , δ p ) from their posterior distributions, and finally evaluate {Pr(y for q = b, a and 1 ≤ j ≤ J. This process is repeated for all remaining individuals and other MCMC draws from the posterior distribution. Finally, the average covariate effect (ACE) for outcome j (1 ≤ j ≤ J) is calculated as the mean of the difference in pointwise probabilities, where (β (m) p , δ (m) p ) is an MCMC draw of (β p , δ p ), M is the number of post burn-in MCMC draws.

logMargLikeOR1
The logarithm of marginal likelihood for OR I model is computed using MCMC outputs from the complete and reduced runs as explained in Section 3. It is reported as a part of the model output, but can also be obtained by calling the model output as shown below. In addition, the package provides logLik.bqrorOR1, an extractor function for displaying the log of marginal likelihood for OR I models. This function is built over the generic function logLik of stats.

R> modelORI$logMargLike
[1] -545.56 devianceOR1 The function devianceOR1 computes the DIC , the effective number of parameters denoted p D , and the deviance (at the posterior mean) for Bayesian quantile regression in OR I model. Additional details on the three quantities can be found in Gelman et al. (2013).
The code below displays the usage of the function. The inputs are covariate matrix including a column of ones (xMat), MCMC draws of β p (betadraws) and δ p (deltadraws), cut-off value (typical values used are 0.05 or 0.10) to determine the lag length and form batches, and the verbose option, which prints the output when set to TRUE (else set to FALSE). Note that if xMat contains names of covariates, then they appear in the summary output; else the summary displays the parameters with numbers in the subscript (such as below).

OR I I model: data, functions, and estimation
Data Generation: The data generation process for the OR I I model closely resembles that of OR I model. In particular, 500 observations are generated for each value of p from the regression model: −4, 6, 5), (x 2 , x 3 ) ∼ U(0, 1) and i ∼ AL(0, σ = 1, p) for i = 1, . . . , n. The continuous values of z are classified based on the cut-points (0, 3) to generate 3 discrete values for y, the outcome variable. Once again, we choose p equal to 0.25, 0.50, and 0.75 to generate three samples from the model, which are referred to as data25j3, data50j3, and data75j3, respectively. The last two letters in the names of the data (i.e., j3) denote the number of unique outcomes in the y variable.

quantregOR2
The function quantregOR2 implements Algorithm 2 and is the main function and for estimating Bayesian quantile regression in OR I I model i.e., an ordinal model with exactly 3 outcomes.
In the code snippet below, we first read the data, define the required quantities, and then call the quantregOR2 for estimating the quantile model. The inputs to the function are as follows: the ordinal outcome variable (y), covariate matrix including a column of ones (xMat), prior mean (b0) and prior covariance matrix (B0) for β p , prior shape (n0) and scale (d0) parameters for σ p , second cut-point (gamma2), burn-in size (burn), post burn-in size (mcmc), quantile (p), and the verbose option which when set to (FALSE) TRUE will (not) print the outputs. Once again, in the code below, we use relatively diffuse prior distributions on (β p , σ p ) to allow the data to speak for itself.

covEffectOR2
The function covEffectOR2 computes the average covariate effect, marginally of the parameters and remaining covariates, for the 3 outcomes of OR I I model at the specified quantiles. Its usage is demonstrated in the application section; however, below we explain the process for computing the average covariate effect. The underlying principle is identical to that of OR I model.
Suppose, we are interested in computing the average covariate effect for the l-th covariate {x i,l } for two different values a and b, and split the covariate and parameter vectors as: given the data y = (y 1 , · · · , y n ) . We again employ the method of composition i.e., randomly select an individual, extract the corresponding sequence of covariate values, draw a value (β p , σ p ) from their posterior distributions, and lastly evaluate {Pr(y for q = b, a, and 1 ≤ j ≤ J = 3. This process is repeated for other individuals and remaining Gibbs draws to compute the ACE for outcome j = 1, 2, 3 as below, where (β (g) is a Gibbs draw of (β p , σ p ) and G is the number of post burn-in Gibbs draws.

logMargLikeOR2
The logarithm of marginal likelihood for OR I I model is computed using Gibbs outputs from the complete and reduced runs as explained in Section 3. It is reported as a part of the model output, but can also be obtained by calling the model output as shown below. For two or more models at the same quantile, the model with a higher log marginal likelihood indicates a better model fit.

R> modelORII$logMargLike
The function devianceOR2 computes the DIC, the effective number of parameters denoted p D , and the deviance (at the posterior mean) for Bayesian quantile regression in OR I I model. Additional details on the three quantities can be found in Gelman et al. (2013).
The code below displays the usage of the function, where the inputs are covariate matrix including a column of ones (xMat), MCMC draws of β p (betadraws) and σ p (sigmadraws), cut-off value (typical values used are 0.05 or 0.10) to determine the lag length and form batches, and the verbose option, which prints the output when set to TRUE (else set to FALSE). Note that if xMat contains names of covariates, then they appear in the summary output; else the summary displays the parameters with numbers in the subscript (such as below).

trace plots
Trace plots are useful for assessing the convergence of MCMC draws and the codes below show how to get the trace plots for (β p , σ p ).

Applications
To illustrate the usage of bqror package on real data applications, we estimate and recreate the results for the educational attainment and tax policy applications presented in Rahman (2016). While the educational attainment study displays the implementation of ordinal quantile regression in OR I model, the tax policy study highlights the use of ordinal quantile regression in OR I I model. Data for both the applications are included as a part of the bqror package.

Educational attainment
In this application, we study the effect of family background, individual level variables, and age cohort on educational attainment of 3923 individuals using data from the National Longitudinal Study of Youth (NLSY, 1979) (Jeliazkov et al., 2008;Rahman, 2016). The dependent variable in the model, education degrees, has four categories: (i) Less than high school, (ii) High school degree, (iii) Some college or associate's degree and (iv) College or graduate degree. A frequency distribution of educational attainment is presented in Figure 1. The independent variables in the model include intercept, square root of family income, mother's education, father's education, mother's working status, gender, race, To estimate a Bayesian ordinal quantile model on educational attainment data, we load the bqror package, prepare the required inputs and feed them into the quantregOR1 function. Specifically, we define the outcome variable, covariate matrix (with covariates in order as in Rahman, 2016), specify the prior distributions on (β p , δ p ), burn-in size, number of post burn-in MCMC iterations, and choose a quantile (p = 0.5 for this illustration).
Next, we utilize the covEffectOR1 function to compute the average covariate effect for a $10,000 increase in family income on the four categories of educational attainment. In general, the calculation of average covariate effect requires creation of either one or two new covariate matrices depending 2 The results reported here are slightly different from those presented in Rahman (2016). This difference in results, aside from lesser number of MCMC draws, is due to a different approach in sampling from the GIG distribution. Rahman (2016) employed the ratio of uniforms method to sample from the GIG distribution (Dagpunar, 2007), while the current paper utilizes the rgig function in the GIGrvg package that overcomes the disadvantages associated with sampling using the ratio of uniforms method (see GIGrvg documentation for further details). Also, see Devroye (2014) for an efficient sampling technique from a GIG distribution. whether the covariate is continuous or indicator (binary), respectively. If the covariate is continuous (such as family income and lets assume this is the lth variable), then a modified covariate matrix is created by adding $10,000 to each observation of family income. This is our xMod2 in the code snippet below and corresponds to the notation x b i,l for i = 1, · · · , n in Subsection 4.1. While the matrix xMod1 (in the code below) is simply the covariate or design matrix xMat (no change here) and so with reference to Subsection 4.1, x a i,l = x i,l for i = 1, · · · , n. When the covariate of interest is an indicator variable, then the xMod1 requires modification as explained in the tax policy application.
We then call the covEffectOR1 function and supply the inputs to get the results.

Tax policy
In this application, the objective is to analyze the factors that affects public opinion on the proposal to raise federal taxes for couples (individuals) earning more than $250,000 ($200,000) per year in the United States (US). The background of this proposal was to extend the "Bush Tax" cuts for the lower and middle income classes, but restore higher rates for the richer class. Such a policy is considered "progrowth," since the motivation is to promote growth in the US economy by augmenting consumption among the low-middle income families. After extensive debate, the proposed policy received a two year extension and formed a part of the "Tax Relief, Unemployment Insurance Reauthorization, and Job Creation Act of 2010".
The data for the study was taken from the 2010-2012 American National Election Studies (ANES) on the Evaluations of Government and Society Study 1 (EGSS 1) and contains 1,164 observations. The dependent variable in the model, individual's opinion on tax increase, has 3 categories: Oppose, Neither favor nor oppose, and Favor (See Figure 3). The covariates included in the model are the intercept, indicator variables for income above $75,000, bachelors' degree, post-bachelors' degree, computer ownership, cell phone ownership, and white race. To fit a Bayesian ordinal quantile model on public opinion about federal tax increase, we load the bqror package, prepare the data, and provide the inputs into the quantregOR2 function. Specifically, we define the outcome variable, covariate matrix (with covariates in order as in Rahman, 2016), specify the prior distributions for (β p , σ p ), choose the second cut-off value, burn-in size, number of post burn-in MCMC iterations, and quantile value (p = 0.5 for this illustration).
R> b0 <-array(rep(0, k), dim = c(k, 1)) R> B0 = 1*diag(k) R> n0 <-5 R> d0 <-8 R> FedTax <-quantregOR2(y = yOrd, x = xMat, b0, B0, + n0, d0, gamma2 = 3, burn = 1125, mcmc = 4500, + p = 0.5, TRUE) [ Finally, we utilize the covEffectOR2 function to demonstrate the calculation of average covariate effect within the OR I I framework. Below, we compute the average covariate effect for computer ownership (assume this is the lth variable) on the 3 categories of public opinion about the tax policy. Here, the covariate is an indicator variable since you may either own a computer (coded as 1) or not (coded as 0). So, we need to create two modified covariate matrices. In the code snippet below, the first matrix xMat1 and the second matrix xMat2 are created by replacing the column on computer ownership with a column of zeros and ones, respectively. With reference to notations in Subsection 4.1 and Subsection 4.2, x a i,l = 0 and x b i,l = 1 for i = 1, · · · , n. We then call the covEffectOR2 function and supply the inputs to get the results.

Conclusion
A wide class of applications in economics, finance, marketing, and the social sciences have dependent variables which are ordinal in nature (i.e., they are discrete and ordered) and are characterized by an underlying continuous variable. Modeling and analysis of such variables have been typically confined to ordinal probit or ordinal logit models, which offers information on the conditional mean of the outcome variable given the covariates. However, a recently proposed method by Rahman (2016) allows Bayesian quantile modeling of ordinal data and thus presents the tool for a more comprehensive analysis and inference. The prevalence of ordinal responses in applications is well known and hence a software package that allows Bayesian quantile analysis with ordinal data will be of immense interest to applied researchers from different fields, including economics and statistics.
The current paper presents an implementation of the bqror package -the only package available for estimation and inference of Bayesian quantile regression in ordinal models. The package offers two MCMC algorithms for estimating ordinal quantile models. An ordinal quantile model with 3 or more outcomes is estimated by a combination of Gibbs sampling and MH algorithm, while estimation of an ordinal quantile model with exactly 3 outcomes utilizes a simpler and computationally faster algorithm that relies solely on Gibbs sampling. For both forms of ordinal quantile models, the bqror package also provides functions to calculate the covariate effects (for continuous as well as binary regressors) and measures for model comparison -marginal likelihood and the DIC. The paper explains how to compute the marginal likelihood from the MCMC outputs and recommends its use over the DIC for model comparison. In addition, this paper demonstrates usage of all the functions for estimation and analysis of Bayesian quantile regression with ordinal data on simulation studies, and the educational attainment and tax policy applications from Rahman (2016).