\SweaveOpts{results=hide, prefix.string=vigpics/oprobit}
\include{zinput}
%\VignetteIndexEntry{Ordinal Probit Regression for Ordered Categorical Dependent Variables}
%\VignetteDepends{Zelig, MASS}
%\VignetteKeyWords{model, probit,regression,ordinal, ordered,categorical}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@

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

\section{{\tt oprobit}: Ordinal Probit Regression for Ordered
Categorical Dependent Variables}\label{oprobit}

Use the ordinal probit regression model if your dependent variables
are ordered and categorical.  They may take on either integer values
or character strings.  For a Bayesian implementation of this model,
see \Sref{oprobit.bayes}.  

\subsubsection{Syntax}

\begin{verbatim}
> z.out <- zelig(as.factor(Y) ~ X1 + X2, model = "oprobit", data = mydata)
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}
If {\tt Y} takes discrete integer values, the {\tt as.factor()}
command will order it automatically.  If {\tt Y} takes on values
composed of character strings, such as ``strongly agree'', ``agree'',
and ``disagree'', {\tt as.factor()} will order the values in the order
in which they appear in {\tt Y}.  You will need to replace your
dependent variable with a factored variable prior to estimating the
model through {\tt zelig()}.  See \Sref{factors} for more information
on creating ordered factors and Example \ref{ord.fact.p} below.

\subsubsection{Example}
\begin{enumerate}
\item {Creating An Ordered Dependent Variable} \label{ord.fact.p}

Load the sample data:  
<<Example.data>>=
 data(sanction)
@ 
Create an ordered dependent variable: 
<<Example.factor>>=
 sanction$ncost <- factor(sanction$ncost, ordered = TRUE,
                           levels = c("net gain", "little effect", 
                           "modest loss", "major loss"))
@
Estimate the model:
<<Example.zelig>>=
 z.out <- zelig(ncost ~ mil + coop, model = "oprobit", data = sanction)
 summary(z.out)
@ 
Set the explanatory variables to their observed values:  
<<Example.setx>>=
 x.out <- setx(z.out, fn = NULL)
@ 
Simulate fitted values given {\tt x.out} and view the results:
<<Example.sim>>=
 s.out <- sim(z.out, x = x.out)
 summary(s.out)
@
%plot does not work but is nnot included in the demo 
\begin{center}
<<label=ExamplePlot,fig=true,echo=false>>= 
# plot(s.out)
@ 
\end{center}

\item {First Differences}

Using the sample data \texttt{sanction}, let us estimate the empirical model and return the coefficients:
<<FirstDifferences.zelig>>=
 z.out <- zelig(as.factor(cost) ~ mil + coop, model = "oprobit", 
                 data = sanction)
@ 
<<FirstDifferences.summary>>=
summary(z.out)
@ 
Set the explanatory variables to their means, with {\tt mil} set
to 0 (no military action in addition to sanctions) in the baseline
case and set to 1 (military action in addition to sanctions) in the
alternative case:
<<FirstDifferences.setx>>=
 x.low <- setx(z.out, mil = 0)
 x.high <- setx(z.out, mil = 1)
@ 
Generate simulated fitted values and first differences, and view the results:
<<FirstDifferences.sim>>=
 s.out <- sim(z.out, x = x.low, x1 = x.high)
@ 
<<FirstDifferences.summary.sim>>= 
summary(s.out)
@
\begin{center}
<<label=FirstDifferencesPlot,fig=true,echo=true>>= 
 plot(s.out)
@ 
\end{center}
\end{enumerate}

\subsubsection{Model}
  Let $Y_i$ be the ordered categorical dependent variable for
  observation $i$ that takes one of the integer values from $1$ to $J$
  where $J$ is the total number of categories.
\begin{itemize}
\item The \emph{stochastic component} is described by an unobserved continuous
  variable, $Y^*_i$, which follows the normal distribution with mean
  $\mu_i$ and unit variance
  \begin{equation*}
    Y_i^* \; \sim \; N(\mu_i, 1). 
  \end{equation*}
  The observation mechanism is 
  \begin{equation*}
    Y_i \; = \; j \quad {\rm if} \quad \tau_{j-1} \le Y_i^* \le \tau_j
    \quad {\rm for} \quad j=1,\dots,J.
  \end{equation*}
  where $\tau_k$ for $k=0,\dots,J$ is the threshold parameter with the
  following constraints; $\tau_l < \tau_m$ for all $l<m$ and
  $\tau_0=-\infty$ and $\tau_J=\infty$.
  
  Given this observation mechanism, the probability for each category,
  is given by
  \begin{equation*}
    \Pr(Y_i = j) \; = \; \Phi(\tau_{j} \mid \mu_i) - \Phi(\tau_{j-1} \mid
    \mu_i) \quad {\rm for} \quad j=1,\dots,J
  \end{equation*}
  where $\Phi(\mu_i)$ is the cumulative distribution function for the
  Normal distribution with mean $\mu_i$ and unit variance.
  
\item The \emph{systematic component} is given by
  \begin{equation*}
    \mu_i \; = \; x_i \beta
  \end{equation*}
  where $x_i$ is the vector of explanatory variables and $\beta$ is
  the vector of coefficients.
\end{itemize}

\subsubsection{Quantities of Interest} 

\begin{itemize}
\item The expected values ({\tt qi\$ev}) for the ordinal probit model
  are simulations of the predicted probabilities for each category:
\begin{equation*}
    E(Y_i = j) \; = \; \Pr(Y_i = j) \; = \; \Phi(\tau_{j} \mid \mu_i)
    - \Phi(\tau_{j-1} \mid  \mu_i) \quad {\rm for} \quad j=1,\dots,J, 
\end{equation*}
given draws of $\beta$ from its posterior.
  
\item The predicted value ({\tt qi\$pr}) is the observed value of
  $Y_i$ given the underlying standard normal distribution described by
  $\mu_i$.

\item The difference in each of the predicted probabilities ({\tt
    qi\$fd}) is given by
  \begin{equation*}
    \Pr(Y=j \mid x_1) \;-\; \Pr(Y=j \mid x) \quad {\rm for} \quad
    j=1,\dots,J.
  \end{equation*}

\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 treated observations in 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 treated observations in 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 \texttt{z.out <- zelig(y \~\,
  x, model = "oprobit", data)}, then you may examine the available
information in \texttt{z.out} by using \texttt{names(z.out)},
see the {\tt coefficients} by using {\tt z.out\$coefficients}, and
a default summary of information through \texttt{summary(z.out)}.
Other elements available through the {\tt \$} operator are listed
below.

\begin{itemize}
\item From the {\tt zelig()} output object {\tt z.out}, you may
  extract:
   \begin{itemize}
   \item {\tt coefficients}: the named vector of coefficients.   
   \item {\tt fitted.values}: an $n \times J$ matrix of the in-sample
     fitted values.
   \item {\tt predictors}: an $n \times (J-1)$ matrix of the linear
     predictors $x_i \beta_j$.
   \item {\tt residuals}: an $n \times (J-1)$ matrix of the residuals. 
   \item {\tt zeta}: a vector containing the estimated class boundaries. 
   \item {\tt df.residual}: the residual degrees of freedom.  
   \item {\tt df.total}: the total degrees of freedom.
   \item {\tt rss}: the residual sum of squares.  
   \item {\tt y}: an $n \times J$ matrix of the dependent variables.
   \item {\tt zelig.data}: the input data frame if {\tt save.data = TRUE}.  
   \end{itemize}

\item From {\tt summary(z.out)}, you may extract:
\begin{itemize}
  \item {\tt coef3}: a table of the coefficients with their associated
    standard errors and $t$-statistics.
  \item {\tt cov.unscaled}: the variance-covariance matrix. 
  \item {\tt pearson.resid}: an $n \times (m-1)$ matrix of the Pearson residuals.  
\end{itemize}

 \item From the {\tt sim()} output object {\tt s.out}, you may extract
   quantities of interest arranged as arrays.  Available quantities
   are:

   \begin{itemize}
   \item {\tt qi\$ev}: the simulated expected probabilities for the
     specified values of {\tt x}, indexed by simulation $\times$
     quantity $\times$ {\tt x}-observation (for more than one {\tt
       x}-observation).
   \item {\tt qi\$pr}: the simulated predicted values drawn from the
     distribution defined by the expected probabilities, indexed by
     simulation $\times$ {\tt x}-observation.
   \item {\tt qi\$fd}: the simulated first difference in the predicted
     probabilities for the values specified in {\tt x} and {\tt x1},
     indexed by simulation $\times$ quantity $\times$ {\tt
       x}-observation (for more than one {\tt x}-observation).
   \item {\tt qi\$att.ev}: the simulated average expected treatment
     effect for the treated from conditional prediction models.  
   \item {\tt 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}
\input{citeZelig}
\subsection*{See also}
The ordinal probit function is part of the VGAM package by Thomas Yee \citep{YeeHas03}. In addition, advanced users may wish to refer to \texttt{help(vglm)} 
in the VGAM library.  Additional documentation is available at
\hlink{http://www.stat.auckland.ac.nz/\~\,yee}{http://www.stat.auckland.ac.nz/~yee}.Sample data are from \cite{Martin92}
\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
%%% TeX-master: t
%%% End: 







