\SweaveOpts{results=hide, prefix.string=vigpics/oprobitBayes}
\include{zinput}
%\VignetteIndexEntry{Bayesian Ordered Probit Regression}
%\VignetteDepends{Zelig}
%\VignetteKeyWords{model, bayes,regression,ordinal, ordered,categorical}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@

<<loadLibrary, echo=F,results=hide>>=
library(Zelig)
@

\section{\texttt{oprobit.bayes}: Bayesian Ordered Probit Regression}

\label{oprobit.bayes}

Use the ordinal probit regression model if your dependent variables are ordered and 
categorical.  They may take either integer values or character strings.  The model 
is estimated using a Gibbs sampler with data augmentation.  For a 
maximum-likelihood implementation of this models, see \Sref{oprobit}.

\subsubsection{Syntax}
\begin{verbatim}
> z.out <- zelig(Y ~ X1 + X2, model = "oprobit.bayes", data = mydata)
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}


\subsubsection{Additional Inputs}

{\tt zelig()} accepts the following arguments to monitor the Markov
chain:  
\begin{itemize}
\item \texttt{burnin}: number of the initial MCMC iterations to be 
 discarded (defaults to 1,000).

\item \texttt{mcmc}: number of the MCMC iterations after burnin
(defaults 10,000).

\item \texttt{thin}: thinning interval for the Markov chain. Only every 
 \texttt{thin}-th draw from the Markov chain is kept. The value of
\texttt{mcmc} must be divisible by this value. The default value is 1.

\item{\texttt{tune}}: tuning parameter for the Metropolis-Hasting step.
The default value is \texttt{NA} which corresponds to 0.05 divided by
the number of categories in the response variable.

\item \texttt{verbose}: defaults to {\tt FALSE}  If \texttt{TRUE},
the progress of the sampler (every $10\%$) is printed to the screen.

\item \texttt{seed}: seed for the random number generator. The default 
is \texttt{NA} which corresponds to a random seed 12345.

\item \texttt{beta.start}: starting values for the Markov 
chain, either a scalar or vector with length equal to the number 
of estimated coefficients. The default is \texttt{NA}, which uses the
maximum likelihood estimates as the starting values.  

\end{itemize}

Use the following parameters to specify the model's priors:  
\begin{itemize}
\item \texttt{b0}: prior mean for the coefficients, either a numeric 
vector or a scalar. If a scalar value, that value will be the prior
mean for all the coefficients. The default is 0.

\item \texttt{B0}: prior precision parameter for the coefficients,
either a square matrix (with dimensions equal to the number of
coefficients) or a scalar. If a scalar value, that value times an
identity matrix will be the prior precision parameter. The default is
0 which leads to an improper prior.
\end{itemize}

\noindent Zelig users may wish to refer to \texttt{help(MCMCoprobit)} 
for more information.

\input{coda_diag}

\subsubsection{Examples}

\begin{enumerate}
\item {Basic Example} \\
Attaching the sample  dataset:
<<BasicExample.data>>=
 data(sanction)
@ 
Estimating ordered probit regression using \texttt{oprobit.bayes}:
<<BasicExample.zelig>>=
 z.out <- zelig(ncost ~ mil + coop, model = "oprobit.bayes",
                  data = sanction, verbose=TRUE)
@ 

Creating an ordered dependent variable:
<<BasicExample.factor>>=
sanction$ncost <- factor(sanction$ncost, ordered = TRUE,
                         levels = c("net gain", "little effect", 
                         "modest loss", "major loss"))
@ 

Checking for convergence before summarizing the estimates:
<<BasicExample.heidel>>=
heidel.diag(z.out$coefficients)
@ 
<<BasicExample.raftery>>= 
raftery.diag(z.out$coefficients)
@  
<<BasicExample.summary>>=
summary(z.out) 
@  
Setting values for the explanatory variables to their sample averages:
<<BasicExample.setx>>=
 x.out <- setx(z.out)
@ 
Simulating quantities of interest from the posterior distribution given:
\texttt{x.out}.
<<BasicExample.sim>>=
 s.out1 <- sim(z.out, x = x.out)
 summary(s.out1)
@
\item {Simulating First Differences} \\
Estimating the first difference (and risk ratio) in the probabilities of
incurring different level of cost when there is no military action versus 
military action while all the other variables held at their 
default values.

<<FirstDifferences.setx>>=
 x.high <- setx(z.out, mil=0)
 x.low <- setx(z.out, mil=1)
@ 
<<FirstDifferences.sim>>= 
s.out2 <- sim(z.out, x = x.high, x1 = x.low)
 summary(s.out2)
@ 
\end{enumerate}

\subsubsection{Model}
Let $Y_{i}$ be the ordered categorical dependent variable for
observation $i$ which takes an integer value $j=1, \ldots, J$.

\begin{itemize}
\item The \emph{stochastic component} is described by an unobserved 
continuous variable, $Y_i^*$, 
\begin{eqnarray*}
Y_{i}^*  \sim \textrm{Normal}(\mu_i, 1).
\end{eqnarray*}
Instead of $Y_i^*$, we observe categorical variable $Y_i$,
\begin{eqnarray*}
Y_i = j \quad \textrm{ if } \tau_{j-1} \le Y_i^* \le \tau_j \textrm{
for } j=1,\ldots, J.
\end{eqnarray*}
where $\tau_j$ for $j=0,\ldots, J$ are the threshold parameters with
the following constraints, $\tau_l < \tau_m$ for $l < m$, and
$\tau_0=-\infty, \tau_J=\infty$.

The probability of observing $Y_i$ equal to category $j$ is,
\begin{eqnarray*}
\Pr(Y_i=j) &=& \Phi(\tau_j \mid \mu_i)-\Phi(\tau_{j-1} \mid \mu_i) 
\textrm{ for } j=1,\ldots, J
\end{eqnarray*}
where $\Phi(\cdot \mid \mu_i)$ is the cumulative distribution function
of the Normal distribution with mean $\mu_i$ and variance 1.

\item The \emph{systematic component} is given by

\begin{eqnarray*}
\mu_{i}= x_i \beta,
\end{eqnarray*}
where $x_{i}$ is the vector of $k$ explanatory variables for 
observation $i$ and $\beta$ is the vector of coefficients.

\item The \emph{prior} for $\beta$ is given by
\begin{eqnarray*}
\beta \sim \textrm{Normal}_k\left(  b_{0},B_{0}^{-1}\right)
\end{eqnarray*}
where $b_{0}$ is the vector of means for the $k$ explanatory variables
and $B_{0}$ is the $k \times k$ precision matrix (the inverse of a
variance-covariance matrix).
\end{itemize}

\subsubsection{Quantities of Interest}

\begin{itemize}
\item The expected values (\texttt{qi\$ev}) for the ordered probit model are
the predicted probability of belonging to each category:
\begin{eqnarray*}
\Pr(Y_i=j)= \Phi(\tau_j \mid x_i \beta)-\Phi(\tau_{j-1} \mid x_i \beta),
\end{eqnarray*}
given the posterior draws of $\beta$ and threshold parameters $\tau$
from the MCMC iterations.

\item The predicted values (\texttt{qi\$pr}) are the observed values of 
$Y_i$ given the observation scheme and the posterior draws of $\beta$
and cut points $\tau$ from the MCMC iterations.

\item The first difference (\texttt{qi\$fd}) in category $j$ for the 
ordered probit model is defined as
\begin{eqnarray*}
\text{FD}_j=\Pr(Y_i=j\mid X_{1})-\Pr(Y_i=j\mid X).
\end{eqnarray*}

\item The risk ratio (\texttt{qi\$rr}) in category $j$ is defined as
\begin{eqnarray*}
\text{RR}_j=\Pr(Y_i=j\mid X_{1})\ /\ \Pr(Y_i=j\mid X).
\end{eqnarray*}


\item In conditional prediction models, the average expected treatment effect
(\texttt{qi\$att.ev}) for the treatment group in category $j$ is
\begin{eqnarray*}
\frac{1}{n_j}\sum_{i:t_{i}=1}^{n_j} \{
Y_{i}(t_{i}=1)-E[Y_{i}(t_{i}=0)] \},
\end{eqnarray*}
where $t_{i}$ is a binary explanatory variable defining the treatment
($t_{i}=1$) and control ($t_{i}=0$) groups, and $n_j$ is the 
number of observations in the treatment group that belong to category $j$.

\item In conditional prediction models, the average predicted treatment effect
(\texttt{qi\$att.pr}) for the treatment group in category $j$ is
\begin{eqnarray*}
\frac{1}{n_j}\sum_{i:t_{i}=1}^{n_j}[Y_{i}(t_{i}=1)-\widehat{Y_{i}(t_{i}=0)}],
\end{eqnarray*}
where $t_{i}$ is a binary explanatory variable defining the treatment
($t_{i}=1$) and control ($t_{i}=0$) groups, and $n_j$ is the 
number of observations in the treatment group that belong to category $j$.
\end{itemize}

\subsubsection{Output Values}

The output of each Zelig command contains useful information which you may
view. For example, if you run:
\begin{verbatim}
z.out <- zelig(y ~ x, model = "oprobit.bayes", data)
\end{verbatim}

\noindent then you may examine the available information in \texttt{z.out} by
using \texttt{names(z.out)}, see the draws from the posterior
distribution of the \texttt{coefficients} by using
\texttt{z.out\$coefficients}, and view a default summary of
information through \texttt{summary(z.out)}. Other elements available
through the \texttt{\$} operator are listed below.

\begin{itemize}
\item From the \texttt{zelig()} output object \texttt{z.out}, you may extract:

\begin{itemize}
\item \texttt{coefficients}: draws from the posterior distributions
of the estimated coefficients $\beta$ and threshold parameters $\tau$.
Note, element $\tau_1$ is normalized to 0 and is not returned in the 
\texttt{coefficients} object.

   \item {\tt zelig.data}: the input data frame if {\tt save.data = TRUE}.  
\item \texttt{seed}: the random seed used in the model.

\end{itemize}

\item From the \texttt{sim()} output object \texttt{s.out}:

\begin{itemize}
\item \texttt{qi\$ev}: the simulated expected values (probabilities) of 
each of the $J$ categories for the specified values of \texttt{x}.

\item \texttt{qi\$pr}: the simulated predicted values (observed values)
 for the specified values of \texttt{x}.

\item \texttt{qi\$fd}: the simulated first difference in the expected
values of each of the $J$ categories for the values specified in 
\texttt{x} and \texttt{x1}.

\item \texttt{qi\$rr}: the simulated risk ratio for the 
expected values of each of the $J$ categories simulated 
from \texttt{x} and \texttt{x1}.

\item \texttt{qi\$att.ev}: the simulated average expected treatment effect
for the treated from conditional prediction models.

\item \texttt{qi\$att.pr}: the simulated average predicted treatment effect
for the treated from conditional prediction models.
\end{itemize}
\end{itemize}


\subsection*{How to Cite}
\input{cites/oprobit.bayes}
\input{citeZelig}
\subsection*{See also}
Bayesian ordinal probit regression is part of the MCMCpack library by Andrew D. Martin and Kevin M. Quinn \citep{MarQui05}.
The convergence diagnostics are part of the CODA library by Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines \citep{PluBesCowVin05}.


\bibliographystyle{asa}
\bibliography{gk,gkpubs}
<<afterpkgs, echo=FALSE>>=
 after<-search()
 torm<-setdiff(after,before)
 for (pkg in torm)
 detach(pos=match(pkg,search()))
@
 \end{document}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:
