#  bqror: An R package for Bayesian Quantile Regression in Ordinal Models

by Prajual Maheshwari and Mohammad Arshad Rahman

###### Abstract

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 $OR_{I}$ model) is estimated by a combination of Gibbs sampling and
Metropolis-Hastings algorithm. Whereas an ordinal model with exactly 3
outcomes (labeled $OR_{II}$ 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.

Keywords: Asymmetric Laplace, Bayesian quantile regression, Gibbs sampling,
marginal likelihood, Markov chain Monte Carlo, ordered choice model.

## 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}$ model111In contrast to Rahman
(2016), our definition of $OR_{I}$ model includes an ordinal model with
exactly 3 outcomes. 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 estimated through a combination of
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_{II}$ 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:

$z_{i}=x^{\prime}_{i}\beta_{p}+\epsilon_{i},\hskip
54.2025pt\forall\;i=1,\cdots,n,$ (1)

where $x^{\prime}_{i}$ is a $1\times k$ vector of covariates, $\beta_{p}$ is a
$k\times 1$ vector of unknown parameters at the $p$-th quantile,
$\epsilon_{i}$ follows an AL distribution i.e., $\epsilon_{i}\sim 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,

$\gamma_{p,j-1}<z_{i}\leq\gamma_{p,j}\;\Rightarrow\;\emph{$y_{i}$ = j},\hskip
54.2025pt\forall\;i=1,\cdots,n;\;j=1,\cdots,J,$ (2)

where
$\gamma_{p}=(\gamma_{p,0}=-\infty,\gamma_{p,1},\ldots,\gamma_{p,J-1},\gamma_{p,J}=\infty)$
is the cut-point vector and J denotes the number of outcomes or categories.
Typically, the cut-point $\gamma_{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},\cdots,y_{n})^{\prime}$,
the likelihood function for the ordinal quantile model can be written as,

$f(y|\Theta_{p})=\prod_{i=1}^{n}\prod_{j=1}^{J}P(y_{i}=j|\Theta_{p})^{I(y_{i}=j)}=\prod_{i=1}^{n}\prod_{j=1}^{J}\bigg{[}F_{AL}(\gamma_{p,j}-x^{\prime}_{i}\beta_{p})-F_{AL}(\gamma_{p,j-1}-x^{\prime}_{i}\beta_{p})\bigg{]}^{I(y_{i}=j)},$
(3)

where $\Theta_{p}=(\beta_{p},\gamma_{p})$, $F_{AL}(\cdot)$ 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,

$z_{i}=x^{\prime}_{i}\beta_{p}+\theta w_{i}+\tau\sqrt{w_{i}}\,u_{i},\hskip
36.135pt\forall\;i=1,\cdots,n,$ (4)

where $\epsilon_{i}=\theta w_{i}+\tau\sqrt{w_{i}}\,u_{i}\sim AL(0,1,p)$,
$w_{i}\sim\mathcal{E}(1)$ is mutually independent of $u_{i}\sim N(0,1)$, $N$
and $\mathcal{E}$ denotes normal and exponential distributions, respectively;
$\theta=(1-2p)/[p(1-p)]$ and $\tau=\sqrt{2/[p(1-p)]}$. Based on this
formulation, we can write the conditional distribution of the latent variable
as $z_{i}|\beta_{p},w_{i}\sim N(x^{\prime}_{i}\beta_{p}+\theta
w_{i},\tau^{2}w_{i})$ for $i=1,\ldots,n$. This allows access to the properties
of normal distribution which helps in constructing efficient MCMC algorithms.

###  $\mathrm{OR_{I}}$ model

The term “$\mathrm{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 $\gamma_{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 $\mathrm{OR_{I}}$ model, we need to preserve the ordering of
cut-points
($\gamma_{p,0}=-\infty<\gamma_{p,1}<\gamma_{p,2}<\ldots<\gamma_{p,J-1}<\gamma_{p,J}=\infty$)
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),

$\delta_{p,j}=\ln(\gamma_{p,j}-\gamma_{p,j-1}),\qquad 2\leq j\leq J-1.$ (5)

This cut-points $(\gamma_{p,1},\gamma_{p,2},\cdots,\gamma_{p,J-1})$ can be
obtained from a one-to-one mapping to $(\delta_{p,2},\cdots,\delta_{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: $\beta_{p}\sim N(\beta_{p0},B_{p0})$, $\delta_{p}\sim
N(\delta_{p0},D_{p0})$ in the bqror package. The augmented joint posterior
distribution for the $\mathrm{OR_{I}}$ model can thus be written as,

$\begin{split}\pi(z,\beta_{p},\delta_{p},w|y)&\propto
f(y|z,\beta_{p},\delta_{p},w)\;\pi(z|\beta_{p},w)\;\pi(w)\;\pi(\beta_{p})\;\pi(\delta_{p}),\\\
&\propto\Big{\\{}\prod_{i=1}^{n}f(y_{i}|z_{i},\delta_{p})\Big{\\}}\;\pi(z|\beta_{p},w)\;\pi(w)\;\pi(\beta_{p})\;\pi(\delta_{p}),\\\
&\propto\prod_{i=1}^{n}\bigg{\\{}\prod_{j=1}^{J}1\\{\gamma_{p,j-1}<z_{i}<\gamma_{p,j}\\}\;N(z_{i}|x^{\prime}_{i}\beta_{p}+\theta
w_{i},\tau^{2}w_{i})\;\mathcal{E}(w_{i}|1)\bigg{\\}}\\\
&\qquad\times\;N(\beta_{p}|\beta_{p0},B_{p0})\;N(\delta_{p}|\delta_{p0},D_{p0}).\end{split}$
(6)

where in the likelihood function of the second line, we use the fact that the
observed $y_{i}$ is independent of $(\beta_{p},w)$ given $(z_{i},\delta_{p})$.
This follows from equation (2) which shows that $y_{i}$ given
$(z_{i},\delta_{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 $\beta_{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
$\delta_{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 $\mathrm{OR_{I}}$ model.

  * •

Sample $\beta_{p}|z,w$ $\sim$ $N(\tilde{\beta}_{p},\tilde{B}_{p})$, where,

  * $\tilde{B}^{-1}_{p}=\bigg{(}\sum_{i=1}^{n}\frac{x_{i}x^{\prime}_{i}}{\tau^{2}w_{i}}+B_{p0}^{-1}\bigg{)}$ and $\tilde{\beta}_{p}=\tilde{B}_{p}\bigg{(}\sum_{i=1}^{n}\frac{x_{i}(z_{i}-\theta w_{i})}{\tau^{2}w_{i}}+B_{p0}^{-1}\beta_{p0}\bigg{)}$.

  * •

Sample $w_{i}|\beta_{p},z_{i}$ $\sim$
$GIG\,(0.5,\tilde{\lambda}_{i},\tilde{\eta})$, for $i=1,\cdots,n$, where,

  * $\tilde{\lambda}_{i}=\Big{(}\frac{z_{i}-x^{\prime}_{i}\beta_{p}}{\tau}\Big{)}^{2}$ and $\tilde{\eta}=\Big{(}\frac{\theta^{2}}{\tau^{2}}+2\Big{)}$.

  * •

Sample $\delta_{p}|y,\beta$ marginally of $w$ (latent weight) and $z$ (latent
data), by generating $\delta_{p}^{\prime}$ using a random-walk chain
$\delta^{\prime}_{p}=\delta_{p}+u$, where $u\sim N(0_{J-2},\iota^{2}\hat{D})$,
$\iota$ is a tuning parameter and $\hat{D}$ denotes negative inverse Hessian,
obtained by maximizing the log-likelihood with respect to $\delta_{p}$. Given
the current value of $\delta_{p}$ and the proposed draw $\delta^{\prime}_{p}$,
return $\delta^{\prime}_{p}$ with probability,

$\alpha_{MH}(\delta_{p},\delta^{\prime}_{p})=\min\bigg{\\{}1,\frac{\;f(y|\beta_{p},\delta^{\prime}_{p})\;\pi(\beta_{p},\delta^{\prime}_{p})}{f(y|\beta_{p},\delta_{p})\;\pi(\beta_{p},\delta_{p})}\bigg{\\}};$

otherwise repeat the old value $\delta_{p}$. The variance of $u$ may be tuned
as needed for appropriate step size and acceptance rate.

  * •

Sample $z_{i}|y,\beta_{p},\gamma_{p},w$ $\sim$
$TN_{(\gamma_{p,j-1},\gamma_{p,j})}(x^{\prime}_{i}\beta_{p}+\theta
w_{i},\tau^{2}w_{i})$ for $i=1,\cdots,n$, where $TN$ denotes a truncated
normal distribution and $\gamma_{p}$ is obtained via $\delta_{p}$ using
equation (5).

###  $\mathrm{OR_{II}}$ model

The term “$\mathrm{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 $\sigma_{p}$ is introduced and the quantile
ordinal model rewritten as follows:

$\begin{split}&z_{i}=x^{\prime}_{i}\beta_{p}+\sigma_{p}\epsilon_{i}=x^{\prime}_{i}\beta_{p}+\sigma_{p}\theta
w_{i}+\sigma_{p}\tau\sqrt{w_{i}}\,u_{i},\hskip
79.49744pt\forall\;i=1,\cdots,n,\\\
&\gamma_{j-1}<z_{i}\leq\gamma_{j}\;\Rightarrow\;\emph{$y_{i}$ = j},\hskip
122.85876pt\forall\;i=1,\cdots,n;\;j=1,2,3,\end{split}$ (7)

where $\sigma_{p}\,\epsilon_{i}\sim AL(0,\sigma_{p},p)$,
($\gamma_{1},\gamma_{2}$) are fixed, and recall $\gamma_{0}=-\infty$ and
$\gamma_{3}=\infty$ . In this formulation, the conditional mean of $z_{i}$ is
dependent on $\sigma_{p}$ which is problematic for Gibbs sampling. So, we
define a new variable $\nu_{i}=\sigma_{p}w_{i}$ and rewrite the model in terms
of $\nu_{i}$. In this representation, $z_{i}|\beta_{p},\sigma_{p},\nu_{i}\sim
N(x^{\prime}_{i}\beta_{p}+\theta\nu_{i},\tau^{2}\sigma_{p}\nu_{i})$, the
conditional mean is free of $\sigma_{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 $\beta_{p}\sim
N(\beta_{p0},B_{p0})$, $\sigma_{p}\sim IG(n_{0}/2,d_{0}/2)$, and
$\nu_{i}\sim\mathcal{E}(\sigma_{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,

$\begin{split}\pi(z,\beta_{p},\nu,\sigma_{p}|y)&\propto
f(y|z,\beta_{p},\nu,\sigma_{p})\;\pi(z|\beta_{p},\nu,\sigma_{p})\;\pi(\nu|\sigma_{p})\;\pi(\beta_{p})\;\pi(\sigma_{p}),\\\
&\propto\Big{\\{}\prod_{i=1}^{n}f(y_{i}|z_{i},\sigma_{p})\Big{\\}}\;\pi(z|\beta_{p},\nu,\sigma_{p})\;\pi(\nu|\sigma_{p})\;\pi(\beta_{p})\;\pi(\sigma_{p}),\\\
&\propto\bigg{\\{}\prod_{i=1}^{n}\prod_{j=1}^{3}1(\gamma_{j-1}<z_{i}<\gamma_{j})\;N(z_{i}|x^{\prime}_{i}\beta_{p}+\theta\nu_{i},\tau^{2}\sigma_{p}\nu_{i})\;\mathcal{E}(\nu_{i}|\sigma_{p})\bigg{\\}}\\\
&\qquad\times\;N(\beta_{p}|\beta_{p0},B_{p0})\;IG(\sigma_{p}|n_{0}/2,d_{0}/2),\end{split}$
(8)

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 $\beta_{p}$ from
an updated multivariate normal distribution and sampling $\sigma_{p}$ from an
updated IG distribution. The latent weight $\nu$ 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 $\mathrm{OR_{II}}$ model.

  * •

Sample $\beta_{p}|z,\sigma_{p},\nu$ $\sim$
$N(\tilde{\beta}_{p},\tilde{B}_{p})$, where,

  * $\tilde{B}^{-1}_{p}=\bigg{(}\sum_{i=1}^{n}\frac{x_{i}x^{\prime}_{i}}{\tau^{2}\sigma_{p}\nu_{i}}+B_{p0}^{-1}\bigg{)}$ and $\tilde{\beta}_{p}=\tilde{B}_{p}\bigg{(}\sum_{i=1}^{n}\frac{x_{i}(z_{i}-\theta\nu_{i})}{\tau^{2}\sigma_{p}\nu_{i}}+B_{p0}^{-1}\beta_{p0}\bigg{)}$.

  * •

Sample $\sigma_{p}|z,\beta_{p},\nu$ $\sim$ $IG(\tilde{n}/2,\tilde{d}/2)$,
where,

  * $\tilde{n}=(n_{0}+3n)$ and $\tilde{d}=\sum_{i=1}^{n}(z_{i}-x^{\prime}_{i}\beta_{p}-\theta\nu_{i})^{2}/\tau^{2}\nu_{i}+d_{0}+2\sum_{i=1}^{n}\nu_{i}$.

  * •

Sample $\nu_{i}|z_{i},\beta_{p},\sigma_{p}$ $\sim$
$GIG(0.5,\tilde{\lambda_{i}},\tilde{\eta})$, for $i=1,\cdots,n$, where,

  * $\tilde{\lambda_{i}}=\frac{(z_{i}-x^{\prime}_{i}\beta_{p})^{2}}{\tau^{2}\sigma_{p}}$ and $\tilde{\eta}=\Big{(}\frac{\theta^{2}}{\tau^{2}\sigma_{p}}+\frac{2}{\sigma_{p}}\Big{)}$.

  * •

Sample $z_{i}|y,\beta_{p},\sigma_{p},\nu_{i}$ $\sim$
$TN_{(\gamma_{j-1},\gamma_{j})}(x^{\prime}_{i}\beta_{p}+\theta\nu_{i},\tau^{2}\sigma_{p}\nu_{i})$
for $i=1,\cdots,n$, and $j=1,2,3$.

## 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 $\mathcal{M}_{s}$ with parameter vector $\Theta_{s}$. Let
$f(y|\mathcal{M}_{s},\Theta_{s})$ be its sampling density, and
$\pi(\Theta_{s}|\mathcal{M}_{s})$ be the prior distribution; where
$s=1,\ldots,S$. Then, the marginal likelihood for the model $\mathcal{M}_{s}$
is given by the expression, $m(y)=\int
f(y|\Theta_{s})\pi(\Theta_{s})\,d\Theta_{s}$, where for notational simplicity
we have suppressed the conditioning on $\mathcal{M}_{s}$. The Bayes factor is
the ratio of marginal likelihoods. So, for any two models $\mathcal{M}_{q}$
versus $\mathcal{M}_{r}$, the Bayes factor,
$B_{qr}=\frac{m(y|\mathcal{M}_{q})}{m(y|\mathcal{M}_{r})}=\frac{\int
f(y|\mathcal{M}_{q},\Theta_{q})\;\pi(\Theta_{q}|\mathcal{M}_{q})\;d\Theta_{q}}{\int
f(y|\mathcal{M}_{r},\Theta_{r})\;\pi(\Theta_{r}|\mathcal{M}_{r})\;d\Theta_{r}}$,
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|\mathcal{M}_{s})$ for
model $\mathcal{M}_{s}$ is expressed as,

$m(y|\mathcal{M}_{s})=\frac{f(y|\mathcal{M}_{s},\Theta_{s})\,\pi(\Theta_{s}|\mathcal{M}_{s})}{\pi(\Theta_{s}|\mathcal{M}_{s},y)}.$
(9)

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 $\Theta^{\ast}$ to minimize
estimation variability. The likelihood ordinate
$f(y|\mathcal{M}_{s},\Theta^{\ast})$ is directly available from the model and
the prior density $\pi(\Theta^{\ast}|\mathcal{M}_{s})$ is assumed by the
researcher. The novel part is the computation of posterior ordinate
$\pi(\Theta^{\ast}|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,

$\ln\hat{m}(y)=\ln
f(y|\Theta^{\ast})+\ln\pi(\Theta^{\ast})-\ln\hat{\pi}(\Theta^{\ast}|y),$ (10)

where, once again, we have dropped the conditioning on $\mathcal{M}_{s}$. The
following two subsections explain the computation of the marginal likelihood
for the $OR_{I}$ and $OR_{II}$ 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:
$\pi(\beta_{p}|z,w)$, $\pi(\delta_{p}|\beta_{p},y)$, $\pi(w|\beta_{p},z)$, and
$\pi(z|\beta_{p},\delta_{p},w,y)$. The conditional posteriors for $\beta_{p}$,
$w$, and $z$ have a known form, but that of $\delta_{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,

$\pi(\beta_{p}^{\ast},\delta_{p}^{\ast}|y)=\pi(\delta_{p}^{\ast}|y)\pi(\beta_{p}^{\ast}|\delta_{p}^{\ast},y),$

where $\Theta^{\ast}=(\beta_{p}^{\ast},\delta_{p}^{\ast})$ 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 $\pi(\delta_{p}^{\ast}|y)$ and then the reduced conditional posterior
ordinate $\pi(\beta_{p}^{\ast}|\delta_{p}^{\ast},y)$.

To obtain an estimate of $\pi(\delta_{p}^{\ast}|y)$, we need to express it in
a computationally convenient form. The parameter $\delta_{p}$ is sampled using
an MH step, which requires specifying a proposal density. Let
$q(\delta_{p},\delta_{p}^{\prime}|\beta_{p},w,z,y)$ denote the proposal
density for the transition from $\delta_{p}$ to $\delta_{p}^{\prime}$, and
let,

$\alpha_{MH}(\delta_{p},\delta^{\prime}_{p})=\min\bigg{\\{}1,\frac{\;f(y|\beta_{p},\delta^{\prime}_{p})\;\pi(\beta_{p})\pi(\delta^{\prime}_{p})}{f(y|\beta_{p},\delta_{p})\;\pi(\beta_{p})\pi(\delta_{p})}\times\frac{q(\delta^{\prime}_{p},\delta_{p}|\beta_{p},w,z,y)}{q(\delta_{p},\delta^{\prime}_{p}|\beta_{p},w,z,y)}\bigg{\\}},$
(11)

denote the probability of making the move. In the context of the model,
$f(y|\beta_{p},\delta_{p})$ is the likelihood given by equation (3),
$\pi(\beta_{p})$ and $\pi(\delta_{p})$ are normal prior distributions (i.e.,
$\beta_{p}\sim N(\beta_{p0},B_{p0})$ and $\delta_{p}\sim
N(\delta_{p0},D_{p0})$ as specified in Section 2.1), and the proposal density
$q(\delta_{p},\delta^{\prime}_{p}|\beta_{p},w,z,y)$ is normal given by
$f_{N}(\delta^{\prime}_{p}|\delta_{p},\iota^{2}\hat{D})$ (see Algorithm 1 in
Section 2.1). There are two points to be noted about the proposal density.
First, the conditioning on $(\beta_{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,

$\pi(\delta_{p}^{\ast}|y)=\frac{E_{1}\\{\alpha_{MH}(\delta_{p},\delta_{p}^{\ast}|\beta_{p},w,z,y)\,q(\delta_{p},\delta_{p}^{\ast}|\beta_{p},w,z,y)\\}}{E_{2}\\{\alpha_{MH}(\delta_{p}^{\ast},\delta_{p}|\beta_{p},w,z,y)\\}},$
(12)

where $E_{1}$ represents expectation with respect to the distribution
$\pi(\beta_{p},\delta_{p},w,z|y)$ and $E_{2}$ represents expectation with
respect to the distribution $\pi(\beta_{p},w,z|\delta_{p}^{\ast},y)\times
q(\delta_{p}^{\ast},\delta_{p}|\beta_{p},w,z,y)$. The quantities in equation
(12) can be estimated using MCMC technique. To estimate the numerator, we take
the draws $\\{\beta_{p}^{(m)},\delta_{p}^{(m)},w^{(m)},z^{(m)}\\}_{m=1}^{M}$
from the _complete MCMC run_ and take an average of the quantity
$\alpha_{MH}(\delta_{p},\delta_{p}^{\ast}|\beta_{p},w,z,y)\,q(\delta_{p},\delta_{p}^{\ast}|\beta_{p},w,z,y)$,
where $\alpha_{MH}(\cdot)$ is given by equation (11) with
$\delta^{\prime}_{p}$ replaced by $\delta_{p}^{\ast}$, and
$q(\delta_{p},\delta_{p}^{\ast}|\beta_{p},w,z,y)$ is given by the normal
density $f_{N}(\delta_{p}^{\ast}|\delta_{p},\iota^{2}\hat{D})$.

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: $\pi(\beta_{p}|w,z)$, $\pi(w|\beta_{p},z)$, and
$\pi(z|\beta_{p},\delta_{p}^{\ast},w,y)$, where note that $\delta_{p}$ is
fixed at $\delta_{p}^{\ast}$ in the MCMC sampling, and thus corresponds to a
_reduced run_. Moreover, at each iteration, we generate

$\delta_{p}^{(h)}\sim
q(\delta_{p}^{\ast},\delta_{p}|\beta_{p}^{(h)},w^{(h)},z^{(h)},y)\equiv
f_{N}(\delta_{p}|\delta_{p}^{\ast},\iota^{2}\hat{D}).$

The resulting quadruplet of draws
$\\{\beta_{p}^{(h)},\delta_{p}^{(h)},w^{(h)},z^{(h)}\\}$, as required, is a
sample from the distribution $\pi(\beta_{p},w,z|\delta_{p}^{\ast},y)\times
q(\delta_{p}^{\ast},\delta_{p}|\beta_{p},w,z,y)$. With the numerator and
denominator now available, the posterior ordinate $\pi(\delta_{p}^{\ast}|y)$
is estimated as,

$\hat{\pi}(\delta_{p}^{\ast}|y)=\frac{M^{-1}\sum_{m=1}^{M}\\{\alpha_{MH}(\delta_{p}^{(m)},\delta_{p}^{\ast}|\beta_{p}^{(m)},w^{(m)},z^{(m)},y)\,q(\delta_{p}^{(m)},\delta_{p}^{\ast}|\beta_{p}^{(m)},w^{(m)},z^{(m)},y)\\}}{H^{-1}\sum_{h=1}^{H}\\{\alpha_{MH}(\delta_{p}^{\ast},\delta_{p}^{(h)}|\beta_{p}^{(h)},w^{(h)},z^{(h)},y)\\}}.$
(13)

The computation of the posterior ordinate
$\pi(\beta_{p}^{\ast}|\delta_{p}^{\ast},y)$ is trivial. We have the sample of
$H$ draws $\\{w^{(h)},z^{(h)}\\}$ from the reduced run, which are marginally
of $\beta_{p}$ from the distribution $\pi(w,z|\delta_{p}^{\ast},y)$. These
draws are utilized to estimate the posterior ordinate as,

$\hat{\pi}(\beta_{p}^{\ast}|\delta_{p}^{\ast},y)=H^{-1}\sum_{h=1}^{H}\pi(\beta_{p}^{\ast}|\delta_{p}^{\ast},w^{(h)},z^{(h)},y).$
(14)

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,

$\ln\hat{m}(y)=\ln
f(y|\beta_{p}^{\ast},\delta_{p}^{\ast})+\ln\Big{[}\pi(\beta_{p}^{\ast})\pi(\delta_{p}^{\ast})\Big{]}-\ln\Big{[}\hat{\pi}(\beta_{p}^{\ast}|y)\,\hat{\pi}(\delta_{p}^{\ast}|\beta_{p}^{\ast},y)\Big{]},$
(15)

where the likelihood $f(y|\beta_{p}^{\ast},\delta_{p}^{\ast})$ (see equation
3) and prior densities are evaluated at
$\Theta^{\ast}=(\beta_{p}^{\ast},\delta_{p}^{\ast})$.

### Marginal likelihood for $OR_{II}$ model

We know from Section 2.2 that the $OR_{II}$ 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
$\pi(\beta_{p}|\sigma_{p},\nu,z)$, $\pi(\sigma_{p}|\beta_{p},\nu,z)$,
$\pi(\nu|\beta_{p},\sigma_{p},z)$, and $\pi(z|\beta_{p},\sigma_{p},\nu,y)$.
However, the variables $(\nu,z)$ are latent. So, we integrate them out and
write the posterior ordinate as
$\pi(\beta_{p}^{\ast},\sigma_{p}^{\ast}|y)=\pi(\beta_{p}^{\ast}|y)\pi(\sigma_{p}^{\ast}|\beta_{p}^{\ast},y)$,
where the terms on the right hand side can be written as,

$\begin{split}\pi(\beta_{p}^{\ast}|y)&=\int\pi(\beta_{p}^{\ast}|\sigma_{p},\nu,z,y)\,\pi(\sigma_{p},\nu,z|y)\,d\sigma_{p}\,d\nu\,dz,\\\
\pi(\sigma_{p}^{\ast}|\beta_{p}^{\ast},y)&=\int\pi(\sigma_{p}^{\ast}|\beta_{p}^{\ast},\nu,z,y)\,\pi(\nu,z|\beta_{p}^{\ast},y)\,dv\,dz,\end{split}$

and $\Theta^{\ast}=(\beta_{p}^{\ast},\sigma_{p}^{\ast})$ denotes a high
density point, such as the mean or median.

The posterior ordinate $\pi(\beta_{p}^{\ast}|y)$ can be estimated as the
ergodic average of the conditional posterior density with the posterior draws
of $(\sigma_{p},\nu,z)$. Therefore, $\pi(\beta_{p}^{\ast}|y)$ is estimated as,

$\hat{\pi}(\beta_{p}^{\ast}|y)=G^{-1}\sum_{g=1}^{G}\pi(\beta_{p}^{\ast}|\sigma_{p}^{(g)},\nu^{(g)},z^{(g)},y).$
(16)

The term $\pi(\sigma_{p}^{\ast}|\beta_{p}^{\ast},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
$\\{\nu^{(g)},z^{(g)}\\}$ from $\pi(\nu,z|\beta_{p}^{\ast},y)$ by successively
sampling from $\pi(\sigma_{p}|\beta_{p}^{\ast},\nu,z)$,
$\pi(\nu|\beta_{p}^{\ast},\sigma_{p},z)$, and
$\pi(z|\beta_{p}^{\ast},\sigma_{p},\nu,y)$, where note that $\beta_{p}$ is
fixed at $\beta_{p}^{\ast}$ in each conditional density. Next, we use the
draws $\\{\nu^{(g)},z^{(g)}\\}$ to compute,

$\hat{\pi}(\sigma_{p}^{\ast}|\beta_{p}^{\ast},y)=G^{-1}\sum_{g=1}^{G}\pi(\sigma_{p}^{\ast}|\beta_{p}^{\ast},\nu^{(g)},z^{(g)},y).$
(17)

which is a simulation consistent estimate of
$\pi(\sigma_{p}^{\ast}|\beta_{p}^{\ast},y)$.

Substituting the two density estimates (equations 16 and 17) into equation
(10), we obtain an estimate of the logarithm of marginal likelihood,

$\ln\hat{m}(y)=\ln
f(y|\beta_{p}^{\ast},\sigma_{p}^{\ast})+\ln\Big{[}\pi(\beta_{p}^{\ast})\pi(\sigma_{p}^{\ast})\Big{]}-\ln\Big{[}\hat{\pi}(\beta_{p}^{\ast}|y)\,\hat{\pi}(\sigma_{p}^{\ast}|\beta_{p}^{\ast},y)\Big{]},$
(18)

where the likelihood function and prior densities are evaluated at
$\Theta^{\ast}=(\beta_{p}^{\ast},\sigma_{p}^{\ast})$. Here, the likelihood
function has the expression,

$f(y|\beta_{p}^{\ast},\sigma_{p}^{\ast})=\prod_{i=1}^{n}\prod_{j=1}^{3}\bigg{[}F_{AL}\left(\frac{\gamma_{j}-x^{\prime}_{i}\beta_{p}^{\ast}}{\sigma_{p}^{\ast}}\right)-F_{AL}\left(\frac{\gamma_{j-1}-x^{\prime}_{i}\beta_{p}^{\ast}}{\sigma_{p}^{\ast}}\right)\bigg{]}^{I(y_{i}=j)},$

where the cut-points $\gamma$ 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.

###  $OR_{I}$ model: data, functions, and estimation

_Data Generation_ : The simulated data for estimation of $OR_{I}$ model is
generated from the regression model: $z_{i}=x^{\prime}_{i}\beta+\epsilon_{i}$,
where $\beta=(-4,5,6)$, $(x_{2},x_{3})\sim U(0,1)$, and $\epsilon_{i}\sim
AL(0,\sigma=1,p)$ for $i=1,\ldots,n$. 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 $\beta_{p}$ and
$\delta_{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 $\beta_{p}\sim
N(0_{k},10\ast I_{k})$ where $0_{k}$ is a zero vector and $I_{k}$ is an
identity matrix of dimension $k\times k$. The prior distribution can be
further diffused (i.e., less informative) by increasing 10 to say 100.
Besides, the prior variance for $\delta_{p}$ needs to be small, such as
$0.25\ast 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).

R> library(’bqror’)R> data("data25j4")R> y <\- data25j4$yR> xMat <\-
data25j4$xR> k <\- dim(xMat)[2]R> J <\- dim(as.array(unique(y)))[1]R> b0 <\-
array(rep(0, k), dim = c(k, 1))R> B0 <\- 10*diag(k)R> d0 <\- array(0, dim =
c(J-2, 1))R> D0 <\- 0.25*diag(J - 2)R> modelORI <\- quantregOR1(y = y, x =
xMat, b0, B0, d0, D0,\+ burn = 1125, mcmc = 4500, p = 0.25, tune = 1,\+
verbose = TRUE)[1] Number of burn-in draws: 1125[1] Number of retained draws:
4500[1] Summary of MCMC draws: Post Mean Post Std Upper Credible Lower
Crediblebeta_1 -3.7244 0.4249 -2.9322 -4.6169beta_2 4.9349 0.5368 6.0444
3.9346beta_3 6.1428 0.5858 7.2996 5.0465delta_1 0.7326 0.1090 0.9380
0.5086delta_2 0.7632 0.0888 0.9363 0.5795[1] MH acceptance rate: 32.16%[1] Log
of Marginal Likelihood: -545.5[1] DIC: 1080.16

The output summarizes the MCMC draws through the posterior mean, posterior
standard deviation, and 95% posterior credible (or probability) interval for
the quantile regression coefficients $\beta$ and transformed cut-points
$\delta$. The output also displays the MH acceptance rate for $\delta$
(preferred acceptance rate is around 30 percent), the logarithm of marginal
likelihood, and the DIC. For two competing models at the same quantile, the
model with a higher (lower) log marginal likelihood (DIC) provides a better
model fit. Besides, after estimation, one may also use the command
modelORI$summary to extract and print the summary output.

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_{i,l}^{a}\\}$ and $\\{x_{i,l}^{b}\\}$, respectively. We split the
covariate and parameter vectors as follows:
$x_{i}^{a}=(x_{i,l}^{a},\,x_{i,-l})$, $x_{i}^{b}=(x_{i,l}^{b},\,x_{i,-l})$,
and $\beta=(\beta_{p,l},\,\beta_{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_{i,l}^{b})-\Pr(y_{i}=j|x_{i,l}^{a})\\}$ for $1\leq j\leq J$,
marginalized over $\\{x_{i,-l}\\}$ and $(\beta_{p},\delta_{p})$, given the
data $y=(y_{1},\cdots,y_{n})^{\prime}$. 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_{i,l}^{b})-\Pr(y_{i}=j|x_{i,l}^{a})\\}$, 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 $(\beta_{p},\delta_{p})$ from their
posterior distributions, and finally evaluate
$\\{\Pr(y_{i}=j|x_{i,l}^{b},x_{i,-l},\beta_{p},\delta_{p})-\Pr(y_{i}=j|x_{i,l}^{a},x_{i,-l},\beta_{p},\delta_{p})\\}$,
where

$\begin{split}\Pr(y_{i}=j|x_{i}^{l},\beta_{p},\delta_{p})&=F_{AL}(\gamma_{p,j}-x_{i,l}^{q}\,\beta_{p,l}-x^{\prime}_{i,-l}\,\beta_{p,-l})-F_{AL}(\gamma_{p,j-1}-x_{i,l}^{q}\,\beta_{p,l}-x^{\prime}_{i,-l}\,\beta_{p,-l}),\end{split}$

for $q=b,a$ and $1\leq j\leq 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\leq j\leq J$) is
calculated as the mean of the difference in pointwise probabilities,

$ACE_{y=j}\approx\frac{1}{M}\frac{1}{n}\sum_{m=1}^{M}\sum_{i=1}^{n}\Big{[}\Pr(y_{i}=j|x_{i,j}^{b},x_{i,-j},\beta_{p}^{(m)},\delta_{p}^{(m)})-\Pr(y_{i}=j|x_{i,j}^{a},x_{i,-j},\beta_{p}^{(m)},\delta_{p}^{(m)})\Big{]},$
(19)

where $(\beta_{p}^{(m)},\delta_{p}^{(m)})$ is an MCMC draw of
$(\beta_{p},\delta_{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).

As shown in the code below, the function requires the following inputs:
observed ordinal outcomes (y), covariate matrix including a column of ones
(xMat), MCMC draws of $\beta_{p}$ (betadraws) and $\delta_{p}$ (deltadraws),
posterior means of $\beta_{p}$ (postMeanbeta) and $\delta_{p}$
(postMeandelta), burn-in size (burn), post burn-in size (mcmc), and the
quantile (p).

R> data("data25j4")R> y <\- data25j4$yR> xMat <\- data25j4$xR> p <\- 0.25R>
mcmc <\- 4500R> burn <\- 0.25*mcmcR> nsim <\- burn + mcmcR> deltadraws <\-
modelORI$deltadrawsR> betadraws <\- modelORI$betadrawsR> postMeanbeta <\-
modelORI$postMeanbetaR> postMeandelta <\- modelORI$postMeandeltaR> allQuantDIC
<\- devianceOR1(y, xMat, betadraws, deltadraws,\+ postMeanbeta, postMeandelta,
burn, mcmc, p)R> allQuantDIC$DIC[1] 1080.159R> allQuantDIC$pd[1] 9.5463R>
allQuantDIC$devpostmean[1] 1061.066

qrnegLogLikensumOR1

The function qrnegLogLikensumOR1, as shown in the code below, requires the
following inputs: ordinal outcome variable (y), covariate matrix including a
column of ones (xMat) , a single draw of $\beta_{p}$ (betaOne) and
$\delta_{p}$ (deltaOne), and the quantile (p). The outputs from this function
are negative of the log-likelihood for each individual and negative of the sum
log-likelihood. The latter when evaluated at postMeanbeta and postMeandelta is
used to calculate DIC and may also be utilized to calculate the Akaike
information criterion (AIC) and Bayesian information criterion (BIC).

R> data("data25j4")R> y <\- data25j4$yR> xMat <\- data25j4$xR> p <\- 0.25R>
betaOne <\- c(0.3990094, 0.8168991, 2.8034963)R> deltaOne <\- c(-0.002570995,
1.044481071)R> output <\- qrnegLogLikensumOR1(y, xMat, betaOne, deltaOne, p)R>
output$negsumlogl[1] 663.5475

infactorOR1

The function infactorOR1 utilizes the batch-means method to compute the
inefficiency factor of $(\beta_{p},\delta_{p})$ based on the MCMC samples
(Greenberg, 2012). Inefficiency factor measures the efficiency of MCMC
algorithm and indicates the cost of working with correlated samples. For
example, an inefficiency factor of 3 implies that it takes 3 correlated draws
to get one independent draw. As such, low inefficiency factor indicates better
mixing and a more efficient MCMC algorithm.

The code below displays the usage of the function. The inputs are covariate
matrix including a column of ones (xMat), MCMC draws of $\beta_{p}$
(betadraws) and $\delta_{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).

R> betadraws <\- modelORI$betadrawsR> deltadraws <\- modelORI$deltadrawsR>
inefficiency <\- infactorOR1(xMat, betadraws, deltadraws, 0.3,verbose=TRUE)[1]
Summary of Inefficiency Factor: Inefficiencybeta_1 2.8711beta_2 2.9982beta_3
3.7465delta_1 5.0465delta_2 2.5527

trace plots

Trace plots are useful for assessing the convergence of MCMC draws. The codes
below illustrate how to obtain the trace plots of MCMC draws for
$(\beta_{p},\delta_{p})$.

R> mcmc <\- 4500R> burn <\- round(0.25*mcmc)R> nsim <\- mcmc + burnR> par(mar
= c(2,5,2,5), mfrow = c(3,2))R> layout(matrix(c(1,2,3,4,5,6), nrow = 3,ncol =
2, byrow = TRUE),\+ widths=c(10,10,10), heights=c(10,10,10,10,10))R>
plot(modelORI$betadraws[1,(burn+1):nsim], type = ’l’, ylab = ’’,\+ xlab = ’’,
main = expression(paste(beta[1])))R> plot(modelORI$betadraws[2,(burn+1):nsim],
type = ’l’, ylab = ’’,\+ xlab = ’’, main = expression(paste(beta[2])))R>
plot(modelORI$betadraws[3,(burn+1):nsim], type = ’l’, ylab = ’’,\+ xlab = ’’,
main = expression(paste(beta[3])))R>
plot(modelORI$deltadraws[1,(burn+1):nsim], type = ’l’, ylab = ’’,\+ xlab = ’’,
main = expression(paste(delta[1])))R>
plot(modelORI$deltadraws[2,(burn+1):nsim], type = ’l’, ylab = ’’,\+ xlab = ’’,
main = expression(paste(delta[2])))

###  $OR_{II}$ model: data, functions, and estimation

_Data Generation_ : The data generation process for the $OR_{II}$ model
closely resembles that of $OR_{I}$ model. In particular, 500 observations are
generated for each value of $p$ from the regression model:
$z_{i}=x^{\prime}_{i}\beta+\epsilon_{i}$, where $\beta=(-4,6,5)$,
$(x_{2},x_{3})\sim U(0,1)$ and $\epsilon_{i}\sim AL(0,\sigma=1,p)$ for
$i=1,\ldots,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_{II}$ 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 $\beta_{p}$, prior shape (n0) and scale (d0)
parameters for $\sigma_{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 $(\beta_{p},\sigma_{p})$ to
allow the data to speak for itself.

R> library(’bqror’)R> data("data25j3")R> y <\- data25j3$yR> xMat <\-
data25j3$xR> k <\- dim(xMat)[2]R> b0 <\- array(rep(0, k), dim = c(k, 1))R> B0
<\- 10*diag(k)R> n0 <\- 5R> d0 <\- 8R> modelORII <\- quantregOR2(y = y, x =
xMat, b0, B0, n0,\+ d0, gamma2 = 3, burn = 1125, mcmc = 4500,\+ p = 0.25,
verbose = TRUE)[1] Number of burn-in draws: 1125[1] Number of retained draws:
4500[1] Summary of MCMC draws: Post Mean Post Std Upper Credible Lower
Crediblebeta_1 -3.9418 0.4690 -3.1029 -4.9207beta_2 5.8494 0.5150 6.8679
4.8611beta_3 4.7777 0.5359 5.9270 3.7950sigma 0.8977 0.0761 1.0551 0.7576[1]
Log of Marginal Likelihood: -404.57[1] DIC: 790.41

The summary output, as displayed above, are the posterior mean, posterior
standard deviation, and 95% posterior credible (or probability) interval for
$(\beta_{p},\sigma_{p})$. The output also exhibits the logarithm of marginal
likelihood and the DIC. After estimation, one may also use the command
modelORII$summary to extract and print the summary output.

covEffectOR2

The function covEffectOR2 computes the average covariate effect, marginally of
the parameters and remaining covariates, for the 3 outcomes of $OR_{II}$ 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:
$x_{i}^{a}=(x_{i,l}^{a},\,x_{i,-l})$, $x_{i}^{b}=(x_{i,l}^{b},\,x_{i,-l})$,
$\beta=(\beta_{p,l},\,\beta_{p,-l})$. We are interested in the distribution of
the difference $\\{\Pr(y_{i}=j|x_{i,l}^{b})-\Pr(y_{i}=j|x_{i,l}^{a})\\}$ for
$1\leq j\leq J=3$, marginalized over $\\{x_{i,-l}\\}$ and
$(\beta_{p},\sigma_{p})$, given the data $y=(y_{1},\cdots,y_{n})^{\prime}$. We
again employ the method of composition i.e., randomly select an individual,
extract the corresponding sequence of covariate values, draw a value
$(\beta_{p},\sigma_{p})$ from their posterior distributions, and lastly
evaluate
$\\{\Pr(y_{i}=j|x_{i,l}^{b},x_{i,-l},\beta_{p},\sigma_{p})-\Pr(y_{i}=j|x_{i,l}^{a},x_{i,-l},\beta_{p},\sigma_{p})\\}$,
where

$\begin{split}\Pr(y_{i}=j|x_{i}^{l},\beta_{p},\sigma_{p})&=F_{AL}\bigg{(}\frac{\gamma_{p,j}-x_{i,l}^{q}\,\beta_{p,l}-x^{\prime}_{i,-l}\,\beta_{p,-l}}{\sigma_{p}}\bigg{)}-F_{AL}\bigg{(}\frac{\gamma_{p,j-1}-x_{i,l}^{q}\,\beta_{p,l}-x^{\prime}_{i,-l}\,\beta_{p,-l}}{\sigma_{p}}\bigg{)}\end{split}$

for $q=b,a$, and $1\leq j\leq 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,

$ACE_{y=j}\approx\frac{1}{G}\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}\Big{[}\Pr(y_{i}=j|x_{i,l}^{b},x_{i,-l},\beta_{p}^{(g)},\sigma_{p}^{(g)})-\Pr(y_{i}=j|x_{i,l}^{a},x_{i,-l},\beta_{p}^{(g)},\sigma_{p}^{(g)})\Big{]},$
(20)

where $(\beta_{p}^{(g)},\sigma_{p}^{(g)})$ is a Gibbs draw of
$(\beta_{p},\sigma_{p})$ and $G$ is the number of post burn-in Gibbs draws.

logMargLikeOR2

The logarithm of marginal likelihood for $OR_{II}$ 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[1] -404.57

devianceOR2

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_{II}$ model. Additional details on the three
quantities can be found in Gelman et al. (2013).

The code below shows the implementation and inputs to the function which are
as follows: observed ordinal outcomes (y), covariate matrix including a column
ones (xMat), Gibbs draws of $\beta_{p}$ (betadraws), Gibbs draws of
$\sigma_{p}$ (sigmadraws), cut-point vector (gammaCp) including
($-\infty,\infty$), posterior mean of $\beta_{p}$ (postMeanbeta) and
$\sigma_{p}$ (postMeansigma), burn-in size (burn), post burn-in size (mcmc),
and quantile (p).

R> library(’bqror’)R> data("data25j3")R> y <\- data25j3$yR> xMat <\-
data25j3$xR> p <\- 0.25R> betadraws <\- modelORII$betadrawsR> sigmadraws <\-
modelORII$sigmadrawsR> gammaCp <\- c(-Inf, 0, 3, Inf)R> postMeanbeta <\-
modelORII$postMeanbetaR> postMeansigma <\- modelORII$postMeansigmaR>
allQuantDIC <\- devianceOR2(y, xMat, betadraws, sigmadraws, gammaCp,\+
postMeanbeta, postMeansigma,\+ burn = 500, mcmc = 4500, p)R>
allQuantDIC$DIC[1] 790.5451R> allQuantDIC$pd[1] 3.926267R>
allQuantDIC$devpostmean[1] 782.6926

qrnegLogLikeOR2

The function qrnegLogLikeOR2 computes the negative sum of log-likelihood for
the $OR_{II}$ model assuming the errors follow $AL(0,\sigma_{p},p)$. As shown
in the code below, the inputs to the function are: observed ordinal outcomes
(y), covariate matrix including a column ones (xMat), cut-point vector
(gammaCp) including ($-\infty,\infty$), a single draw of $\beta_{p}$ (betaOne)
and $\sigma_{p}$ (sigmaOne), and quantile (p). When evaluated at posterior
means (i.e., postMeanbeta and postMeansigma), the output (i.e., negative sum
of the log-likelihood) is utilized in calculation of DIC and may also be used
to compute the AIC and BIC.

R> library(’bqror’)R> data("data25j3")R> y <\- data25j3$yR> xMat <\-
data25j3$xR> p <\- 0.25R> gammaCp <\- c(-Inf, 0, 3, Inf)R> betaOne <\-
c(1.7201671, 1.9562172, 0.8334668)R> sigmaOne <\- 0.9684741R> output <\-
qrnegLogLikeOR2(y, xMat, gammaCp, betaOne, sigmaOne, p)R> output[1] 567.2358

infactorOR2

The function infactorOR1 utilizes the batch-means method to compute the
inefficiency factor of $(\beta_{p},\sigma_{p})$ based on the Gibbs draws
(Greenberg, 2012).

The code below displays the usage of the function, where the inputs are
covariate matrix including a column of ones (xMat), MCMC draws of $\beta_{p}$
(betadraws) and $\sigma_{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).

R> betadraws <\- modelORII$betadrawsR> sigmadraws <\- modelORII$sigmadrawsR>
inefficiency <\- infactorOR2(xMat, betadraws, sigmadraws, 0.1,\+
verbose=TRUE)[1] Summary of Inefficiency Factor: Inefficiencybeta_1
3.3937beta_2 3.0315beta_3 3.1486sigma 3.4146

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 $(\beta_{p},\sigma_{p})$.

R> mcmc <\- 4500R> burn <\- round(0.25*mcmc)R> nsim <\- mcmc + burnR> par(mar
= c(2,5,2,5), mfrow = c(2,2))R> layout(matrix(c(1,2,3,4), nrow = 2,ncol = 2,
byrow = TRUE),\+ widths=c(10,10,10), heights=c(10,10,10,10,10))R>
plot(modelORII$betadraws[1,(burn+1):nsim], type = ’l’, ylab = ’’,\+ xlab = ’’,
main = expression(paste(beta[1])))R>
plot(modelORII$betadraws[2,(burn+1):nsim], type = ’l’, ylab = ’’,\+ xlab = ’’,
main = expression(paste(beta[2])))R>
plot(modelORII$betadraws[3,(burn+1):nsim], type = ’l’, ylab = ’’,\+ xlab = ’’,
main = expression(paste(beta[3])))R>
plot(modelORII$sigmadraws[1,(burn+1):nsim], type = ’l’, ylab = ’’,\+ xlab =
’’, main = expression(paste(sigma[1])))

## 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_{II}$ 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, indicator variables to point
whether the youth lived in an urban area or South at the age of 14, and three
indicator variables to indicate the individual’s age in 1979 (serves as a
control for age cohort effects).

Figure 1:  Frequency distribution for educational attainment. The four
categories, denoted by Cat 1, Cat 2, Cat 3, and Cat 4 correspond to less than
high school, high school degree, some college or associate’s degree, and
college or graduate degree, respectively. For each vertical bar, the number
inside the box (at the top) denote the number of observations (percentage) for
that category.

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 $(\beta_{p},\delta_{p})$, burn-in size, number of post burn-
in MCMC iterations, and choose a quantile ($p=0.5$ for this illustration).

R> library(bqror)R> data("Educational_Attainment")R> data <\-
na.omit(Educational_Attainment)

R> data$fam_income_sqrt <\- sqrt(data$fam_income)R> cols <\-
c("mother_work","urban","south",\+
"father_educ","mother_educ","fam_income_sqrt","female",\+
"black","age_cohort_2","age_cohort_3","age_cohort_4")R> x <\- data[cols]R>
x$intercept <\- 1R> xMat <\- x[,c(12,6,5,4,1,7,8,2,3,9,10,11)]R> k <\-
dim(xMat)[2]R> yOrd <\- data$dep_edu_levelR> J <\-
dim(as.array(unique(yOrd)))[1]

R> b0 <\- array(rep(0, k), dim = c(k, 1))R> B0 <\- 1*diag(k)R> d0 <\- array(0,
dim = c(J-2, 1))R> D0 <\- 0.25*diag(J - 2)R> p <\- 0.5R> EducAtt <\-
quantregOR1(y = yOrd, x = xMat, b0, B0,\+ d0, D0, burn = 1125, mcmc = 4500, p,
tune=1, TRUE)[1] Number of burn-in draws: 1125[1] Number of retained draws:
4500[1] Summary of MCMC draws: Post Mean Post Std Upper Credible Lower
Credibleintercept -3.2231 0.2276 -2.7795 -3.6797fam_income_sqrt 0.2796 0.0258
0.3321 0.2300mother_educ 0.1219 0.0189 0.1601 0.0855father_educ 0.1858 0.0154
0.2166 0.1565mother_work 0.0697 0.0798 0.2231 -0.0913female 0.3493 0.0788
0.5066 0.1992black 0.4372 0.1021 0.6367 0.2406urban -0.0826 0.0951 0.1046
-0.2728south 0.0819 0.0874 0.2542 -0.0865age_cohort_2 -0.0365 0.1229 0.1998
-0.2777age_cohort_3 -0.0500 0.1260 0.1988 -0.2874age_cohort_4 0.5005 0.1328
0.7493 0.2343delta_1 0.8979 0.0286 0.9596 0.8434delta_2 0.5455 0.0346 0.6105
0.4776[1] MH acceptance rate: 27%[1] Log of Marginal Likelihood: -4925.14[1]
DIC: 9788.14

The posterior results222The 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. from the MCMC draws are
summarized above. In the summary, we report the posterior mean and posterior
standard deviation of the parameters ($\beta_{p},\delta_{p}$), and their 95%
posterior credible interval. Additionally, the summary displays the MH
acceptance rate of $\delta_{p}$, the logarithm of marginal likelihood, and the
DIC. Figure 2 presents the trace plots of the MCMC, which can be obtained
using codes similar to those demonstrated in the simulation studies.

Figure 2:  Trace plots of the MCMC draws in the educational attainment study.

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
whether the covariate is continuous or indicator (binary), respectively. If
the covariate is continuous (such as family income and lets assume this is the
$l$th 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_{i,l}^{b}$ for $i=1,\cdots,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_{i,l}^{a}=x_{i,l}$ for $i=1,\cdots,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.

R> xMat1 <\- xMatR> xMat2 <\- xMatR> xMat2$fam_income_sqrt <\-
sqrt((xMat1$fam_income_sqrt)^2 + 10)R> EducAttCE <\- covEffectOR1(EducAtt,
yOrd, xMat1, xMat2, p = 0.5, verbose = TRUE)[1] Summary of Covariate Effect:
Covariate EffectCategory_1 -0.0315Category_2 -0.0132Category_3
0.0179Category_4 0.0268

The results shows that for a $10,000 increase in family income, at the 50th
quantile the probability of obtaining _less than high school_ ( _high school
degree_) decreases by 3.15 percent (1.32 percent), while the probability of
achieving _some college or associate’s degree_ (_College or graduate degree_)
increases by 1.79 percent (2.68 percent).

### Tax policy

Figure 3: Frequency distribution for individual’s opinion on tax increase. The
3 categories, denoted by Cat 1, Cat 2, and Cat 3, correspond to oppose,
neither favor nor oppose, and favor the tax increase. For each vertical bar,
the number inside the box (at the top) denote the number of observations
(percentage) for that category.

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 “pro-growth,” 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 $(\beta_{p},\sigma_{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> library(bqror)R> data("Policy_Opinion")R> data <\- na.omit(Policy_Opinion)

R> cols <\- c("Intercept","AgeCat","IncomeCat",\+
"Bachelors","Post.Bachelors","Computers","CellPhone",\+ "White")R> x <\-
data[cols]R> k <\- dim(x)[2]R> xMat <\- x[,c(1,2,3,4,5,6,7,8)]R> yOrd <\-
data$y

The results (See Footnote 2) from the MCMC draws are summarized below. In the
summary, we report the posterior mean, posterior standard deviation, and the
95% posterior credible interval of the parameters ($\beta_{p},\sigma_{p}$).
Additionally, the summary displays the logarithm of marginal likelihood and
the DIC. Figure 4 presents the trace plots of the MCMC draws, which can be
obtained using codes similar to those presented in the simulation studies.

R> b0 <\- array(rep(0, k), dim = c(k, 1))R> B0 = 1*diag(k)R> n0 <\- 5R> d0 <\-
8R> FedTax <\- quantregOR2(y = yOrd, x = xMat, b0, B0,\+ n0, d0, gamma2 = 3,
burn = 1125, mcmc = 4500,\+ p = 0.5, TRUE)[1] Number of burn-in draws :
1125[1] Number of retained draws : 4500[1] Summary of MCMC draws : Post Mean
Post Std Upper Credible Lower CredibleIntercept 1.9869 0.4454 2.8546
1.0936AgeCat 0.2477 0.2980 0.8481 -0.3401IncomeCat -0.4968 0.3364 0.1492
-1.1767Bachelors 0.0684 0.3678 0.7790 -0.6423Post.Bachelors 0.4726 0.4646
1.4041 -0.4232Computers 0.7204 0.3653 1.4561 0.0005CellPhone 0.8595 0.4010
1.6283 0.0606White 0.0814 0.3762 0.8093 -0.6704sigma 2.2268 0.1434 2.5164
1.9497[1] Log of Marginal Likelihood: -1173.97[1] DIC: 2334.76

Figure 4:  Trace plots for the MCMC draws in the tax policy study.

Finally, we utilize the covEffectOR2 function to demonstrate the calculation
of average covariate effect within the $OR_{II}$ framework. Below, we compute
the average covariate effect for computer ownership (assume this is the $l$th
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_{i,l}^{a}=0$ and $x_{i,l}^{b}=1$ for
$i=1,\cdots,n$. We then call the covEffectOR2 function and supply the inputs
to get the results.

R> xMat1 <\- xMatR> xMat1$Computers <\- 0R> xMat2 <\- xMatR> xMat2$Computers
<\- 1R> FedTaxCE <\- covEffectOR2(FedTax, yOrd, xMat1, xMat2, gamma2 = 3, p =
0.5,\+ verbose = TRUE)[1] Summary of Covariate Effect: Covariate
EffectCategory_1 -0.0397Category_2 -0.0330Category_3 0.0727

The results on covariate effect shows that at the 50th quantile, ownership of
computer decreases probability for the first category _Oppose_ (_Neither favor
nor oppose_) by 3.97 percent (3.30 percent), and increases the probability for
the third category _Favor_ by 7.27 percent.

## 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).

## References

  * Albert and Chib (1993) J. Albert and S. Chib.  Bayesian analysis of binary and polychotomous response data.  _Journal of the American Statistical Association_ , 88(422):669–679, 1993.  doi: https://doi.org/10.1080/01621459.1993.10476321.
  * Albert and Chib (2001) J. Albert and S. Chib.  Sequential ordinal modeling with applications to survival data.  _Biometrics_ , 57:829–836, 2001.  doi: https://www.jstor.org/stable/3068422.
  * Alhamzawi (2016) R. Alhamzawi.  Bayesian model selection in ordinal quantile regression.  _Computational Statistics and Data Analysis_ , 103:68–78, 2016.  doi: https://doi.org/10.1016/j.csda.2016.04.014.
  * Alhamzawi and Ali (2018) R. Alhamzawi and H. T. M. Ali.  Bayesian quantile regression for ordinal longitudinal data.  _Journal of Applied Statistics_ , 45(5):815–828, 2018.  doi: https://doi.org/10.1080/02664763.2017.1315059.
  * Benoit and Van den Poel (2010) D. F. Benoit and D. Van den Poel.  Binary quantile regression: A Bayesian approach based on the asymmetric Laplace distribution.  _Journal of Applied Econometrics_ , 27(7):1174–1188, 2010.  doi: https://doi.org/10.1002/jae.1216.
  * Benoit and Van den Poel (2017) D. F. Benoit and D. Van den Poel.  bayesqr: A Bayesian approach to quantile regression.  _Journal of Statistical Software_ , 76(7), 2017.  URL https://cran.r-project.org/web/packages/bayesQR/index.html.  R package version 2.3.
  * Bresson et al. (2021) G. Bresson, G. Lacroix, and M. A. Rahman.  Bayesian panel quantile regression for binary outcomes with correlated random effects: An application on crime recidivism in Canada.  _Empirical Economics_ , 60(1):227–259, 2021.  doi: https://doi.org/10.1007/s00181-020-01893-5.
  * Casella and George (1992) G. Casella and E. I. George.  Explaining the Gibbs sampler.  _The American Statistician_ , 46(3):167–174, 1992\.  doi: https://doi.org/10.1080/00031305.1992.10475878.
  * Chib (1995) S. Chib.  Marginal likelihood from the Gibbs output.  _Journal of the American Statistical Association_ , 90(432):1313–1321, 1995.  doi: https://doi.org/10.1080/01621459.1995.10476635.
  * Chib and Greenberg (1995) S. Chib and E. Greenberg.  Understanding the Metropolis-Hastings algorithm.  _The American Statistician_ , 49:327–335, 1995.  doi: https://doi.org/10.1080/00031305.1995.10476177.
  * Chib and Jeliazkov (2001) S. Chib and I. Jeliazkov.  Marginal likelihood from the Metropolis-Hastings output.  _Journal of the American Statistical Association_ , 96(453):270–281, 2001.  doi: https://doi.org/10.1198/016214501750332848.
  * Chib and Jeliazkov (2006) S. Chib and I. Jeliazkov.  Inference in semiparametric dynamic models for binary longitudinal data.  _Journal of the American Statistical Association_ , 101(474):685–700, 2006.  doi: https://doi.org/10.1198/016214505000000871.
  * Dagpunar (2007) J. Dagpunar.  _Simulations and Monte Carlo: With Applications in Finance and MCMC_.  John Wiley & Sons Ltd., UK, 2007.  doi: https://onlinelibrary.wiley.com/doi/book/10.1002/9780470061336.
  * Davino et al. (2014) C. Davino, M. Furno, and D. Vistocco.  _Quantile Regression: Theory and Applications_.  John Wiley & Sons, Chichester, 2014.  doi: https://doi.org/10.1002/9781118752685.
  * Devroye (2014) L. Devroye.  Random variate generation for the generalized inverse Gaussian distribution.  _Statistics and Computing_ , 24(2):239–246, 2014\.  doi: https://doi.org/10.1007/s11222-012-9367-z.
  * Furno and Vistocco (2018) M. Furno and D. Vistocco.  _Quantile Regression: Estimation and Simulation_.  John Wiley & Sons, New Jersey, 2018.  doi: https://doi.org/10.1002/9781118863718.
  * Gelman et al. (2013) A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin.  _Bayesian Data Analysis_.  Chapman & Hall, New York, 2013.  doi: https://doi.org/10.1201/b16018.
  * Geman and Geman (1984) S. Geman and D. Geman.  Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 6(6):721–741, 1984.  doi: https://doi.org/10.1109/TPAMI.1984.4767596.
  * Geraci (2014) M. Geraci.  Linear quantile mixed models: The lqmm package for Laplace quantile regression.  _Journal of Statistical Software_ , 57(13), 2014.  URL https://cran.r-project.org/web/packages/lqmm/index.html.  R package version 1.5.5.
  * Geraci and Bottai (2014) M. Geraci and M. Bottai.  Linear quantile mixed models.  _Statistics and Computing_ , 24(461-479), 2014.  doi: https://doi.org/10.1007/s11222-013-9381-9.
  * Ghasemzadeh et al. (2018) S. Ghasemzadeh, M. Ganjali, and T. Baghfalaki.  Bayesian quantile regression for analyzing ordinal longitudinal responses in the presence of non-ignorable missingness.  _METRON_ , 76(3):321–348, 2018.  doi: https://doi.org/10.1007/s40300-018-0136-4.
  * Ghasemzadeh et al. (2020) S. Ghasemzadeh, M. Ganjali, and T. Baghfalaki.  Bayesian quantile regression for joint modeling of longitudinal mixed ordinal continuous data.  _Communications in Statistics $-$ Simulation and Computation_, 49(2):1375–1395, 2020.  doi: https://doi.org/10.1080/03610918.2018.1484482.
  * Greenberg (2012) E. Greenberg.  _Introduction to Bayesian Econometrics_.  2nd Edition, Cambridge University Press, New York, 2012.  doi: https://doi.org/10.1017/CBO9780511808920.
  * Greene and Hensher (2010) W. H. Greene and D. A. Hensher.  _Modeling Ordered Choices: A Primer_.  Cambridge University Press, Cambridge, 2010.  doi: https://doi.org/10.1017/CBO9780511845062.
  * Jeliazkov and Rahman (2012) I. Jeliazkov and M. A. Rahman.  Binary and ordinal data analysis in economics: Modeling and estimation.  In X. S. Yang, editor, _Mathematical Modeling with Multidisciplinary Applications_ , pages 123–150. John Wiley & Sons Inc., New Jersey, 2012.  doi: https://doi.org/10.1002/9781118462706.ch6.
  * Jeliazkov et al. (2008) I. Jeliazkov, J. Graves, and M. Kutzbach.  Fitting and comparison of models for multivariate ordinal outcomes.  _Advances in Econometrics: Bayesian Econometrics_ , 23:115–156, 2008.  doi: https://doi.org/10.1016/S0731-9053(08)23004-5.
  * Johnson and Albert (2000) V. E. Johnson and J. H. Albert.  _Ordinal Data Modeling_.  Springer, New York, 2000.  doi: https://doi.org/10.1007/b98832.
  * Koenker (2005) R. Koenker.  _Quantile Regression_.  Cambridge University Press, Cambridge, 2005.  doi: https://doi.org/10.1257/jep.15.4.143.
  * Koenker (2021) R. Koenker.  _quantreg: Quantile Regression_ , 2021.  URL https://cran.r-project.org/web/packages/quantreg/index.html.  R package version 5.85.
  * Koenker and Bassett (1978) R. Koenker and G. Bassett.  Regression quantiles.  _Econometrica_ , 46(1):33–50, 1978.  doi: https://doi.org/10.2307/1913643.
  * Kordas (2006) G. Kordas.  Smoothed binary regression quantiles.  _Journal of Applied Econometrics_ , 21(3):387–407, 2006.  doi: https://doi.org/10.1002/jae.843.
  * Kozumi and Kobayashi (2011) H. Kozumi and G. Kobayashi.  Gibbs sampling methods for Bayesian quantile regression.  _Journal of Statistical Computation and Simulation_ , 81(11):1565–1578, 2011.  doi: https://doi.org/10.1080/00949655.2010.496117.
  * Mukherjee and Rahman (2016) D. Mukherjee and M. A. Rahman.  To drill or not to drill? an econometric analysis of US public opinion.  _Energy Policy_ , 91:341–351, 2016.  doi: 10.1016/j.enpol.2015.11.023.
  * Poirier (1995) D. J. Poirier.  _Intermediate Statistics and Econometrics_.  The MIT Press, Cambridge, 1995.
  * Rahman (2016) M. A. Rahman.  Bayesian quantile regression for ordinal models.  _Bayesian Analysis_ , 11(1):1–24, 2016.  doi: https://doi.org/10.1214/15-BA939.
  * Rahman and Karnawat (2019) M. A. Rahman and S. Karnawat.  Flexible Bayesian quantile regression in ordinal models.  _Advances in Econometrics_ , 40B:211–251, 2019.  doi: https://doi.org/10.1108/S0731-90532019000040B011.
  * Rahman and Vossmeyer (2019) M. A. Rahman and A. Vossmeyer.  Estimation and applications of quantile regression for binary longitudinal data.  _Advances in Econometrics_ , 40B:157–191, 2019.  doi: https://doi.org/10.1108/S0731-90532019000040B009.
  * Spiegelhalter et al. (2002) D. J. Spiegelhalter, N. G. Best, B. P. Carlin, and A. Van Der Linde.  Bayesian measures of model complexity and fit.  _Journal of the Royal Statistical Society – Series B_ , 64(4):583–639, 2002.  doi: https://doi.org/10.1111/1467-9868.00353.
  * Tian et al. (2021) Y.-Z. Tian, M.-L. Tang, W.-S. Chan, and M.-Z. Tian.  Bayesian bridge-randomized penalized quantile regression for for ordinal longitudinal data, with application to firm’s bond ratings.  _Computational Statistics_ , 36:1289–1319, 2021.  doi: https://doi.org/10.1007/s00180-020-01037-4.
  * Yu and Moyeed (2001) K. Yu and R. A. Moyeed.  Bayesian quantile regression.  _Statistics and Probability Letters_ , 54(4):437–447, 2001.  doi: https://doi.org/10.1016/S0167-7152(01)00124-9.
  * Yu and Zhang (2005) K. Yu and J. Zhang.  A three parameter asymmetric Laplace distribution and its extensions.  _Communications in Statistics – Theory and Methods_ , 34(9-10):1867–1879, 2005.  doi: https://doi.org/10.1080/03610920500199018.
  * Zu and Yu (2021) T. Zu and Y. Yu.  SIQR: An R package for single-index quantile regression.  _The R Journal_ , 2021.  doi: 10.32614/RJ-2021-092.  URL https://journal.r-project.org/archive/2021/RJ-2021-092/index.html.  R package version 2.3.

_Prajual Maheshwari
Quantative Analyst, Dolat Capital, India
B.S. (Economics, 2016-20), Department of Economic Sciences
Indian Institute of Technology Kanpur, India
<EMAIL_ADDRESS>
_Mohammad Arshad Rahman
Associate Professor, College of Business
Zayed University, Abu Dhabi, UAE
<EMAIL_ADDRESS>
&
Associate Professor, Department of Economic Sciences
Indian Institute of Technology Kanpur, India
ORCID: https://orcid.org/0000-0001-8434-0042
<EMAIL_ADDRESS><EMAIL_ADDRESS>_