\SweaveOpts{results=hide, prefix.string=vigpics/mlogit}
\include{zinput}
%\VignetteIndexEntry{Multinomial Logistic Regression for Dependent Variables with Unordered Categorical Values}
%\VignetteDepends{Zelig, VGAM}
%\VignetteKeyWords{model, bayes,multinomial, logistic regression}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@

<<loadLibrary, echo=F,results=hide>>=
library(Zelig)
@ 
\section{{\tt mlogit}: Multinomial Logistic Regression for
Dependent Variables with Unordered Categorical Values}\label{mlogit}
Use the multinomial logit distribution to model unordered categorical
variables.  The dependent variable may be in the format of either
character strings or integer values.  See 
%\Sref{mlogit.bayes} 
for a Bayesian version of this model.  

\subsubsection{Syntax}

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

\subsubsection{Input Values}

If the user wishes to use the same formula across all levels, then 
\verb|formula <- as.factor(Y) ~ X1 + X2| may be used. 
If the user wants to use different formula for each level then the following 
syntax should be used:
\begin{verbatim}
formulae <- list(list(id(Y, "apples")~ X1, 
                      id(Y, "bananas")~ X1 + X2)
\end{verbatim}
where Y above is supposed to be a factor variable with levels {apples,bananas,oranges}.  
By default, oranges is the last level and omitted.  (You cannot
specify a different base level at this time.)  
For $J$ equations, there must be $J + 1$ levels.
 
\subsubsection{Examples} \label{ternary}

\begin{enumerate}

\item {The same formula for each level}

Load the sample data:  
<<Examples.data>>=
 data(mexico)
@ 
Estimate the empirical model:  
<<Examples.zelig>>=
 z.out1 <- zelig(as.factor(vote88) ~ pristr + othcok + othsocok, 
                  model = "mlogit", data = mexico)
@ 
Set the explanatory variables to their default values, with {\tt pristr}
(for the strength of the PRI) equal to 1 (weak) in the baseline values, and
equal to 3 (strong) in the alternative values:  
<<Examples.setx>>=
 x.weak <- setx(z.out1, pristr = 1)
 x.strong <- setx(z.out1, pristr = 3)
@ 
Generate simulated predicted probabilities {\tt qi\$ev} and differences in
the predicted probabilities {\tt qi\$fd}:
<<Examples.sim>>=
 s.out1 <- sim(z.out1, x = x.strong, x1 = x.weak)

@ 
<<Examples.summary>>= 
summary(s.out1)
@ 
Generate simulated predicted probabilities {\tt qi\$ev} for the
alternative values:  
<<Examples.weak>>=
 ev.weak <- s.out1$qi$ev + s.out1$qi$fd

@ 
Plot the differences in the predicted probabilities.  
<<Examples.library>>=
 library(vcd)
@
%%%%does not work but the demo gives same error message
%%%funlegend=ternarypoints,legendargs=parpoint
\begin{center}
<<label=TernaryPlot,fig=true, echo=true>>=
ternaryplot(x=s.out1$qi$ev, pch = ".", col = "blue",main = "1988 Mexican Presidential Election")

###ternarypoints(ev.weak, pch = ".", col = "red")
@
\end{center}
\item {Different formula for each level}

Estimate the empirical model:  
<<Levels.zelig>>=
 z.out2 <- zelig(list(id(vote88, "1") ~ pristr + othcok,
                       id(vote88, "2") ~ othsocok), 
                  model = "mlogit", data = mexico)
@ 
Set the explanatory variables to their default values, with {\tt pristr}
(for the strength of the PRI) equal to 1 (weak) in the baseline values, and
equal to 3 (strong) in the alternative values: 
<<Levels.setx>>=
 x.weak <- setx(z.out2, pristr = 1)
 x.strong <- setx(z.out2, pristr = 3)
@ 
Generate simulated predicted probabilities {\tt qi\$ev} and differences in
the predicted probabilities {\tt qi\$fd}:
<<Levels.sim>>=
 s.out1 <- sim(z.out2, x = x.strong, x1 = x.weak)
@
<<Levels.summary>>= 
 summary(s.out1)
@ 
Generate simulated predicted probabilities {\tt qi\$ev} for the
alternative values:  
<<Levels.weak>>=
 ev.weak <- s.out1$qi$ev + s.out1$qi$fd
@ 
Using the vcd package, plot the differences in the predicted probabilities.
\begin{center}
<<label=LevelsTernaryPlot,fig=true,echo=true>>=
ternaryplot(s.out1$qi$ev, pch = ".", col = "blue", main = "1988 Mexican Presidential Election")
###ternarypoints(ev.weak, pch = ".", col = "red")
@
\end{center} 

 
\end{enumerate}

\subsubsection{Model}
Let $Y_i$ be the unordered categorical dependent variable that takes
one of the values from 1 to $J$, where $J$ is the total number of
categories.

\begin{itemize}
\item The stochastic component is given by
  \begin{equation*}
    Y_i \; \sim \; \textrm{Multinomial}(y_{i} \mid \pi_{ij}), 
  \end{equation*}
  where $\pi_{ij}=\Pr(Y_i=j)$ for $j=1,\dots,J$. 
  
\item The systemic component is given by:  
  \begin{equation*}
    \pi_{ij}\; = \; \frac{\exp(x_{i}\beta_{j})}{\sum^{J}_{k = 1}
      \exp(x_{i}\beta_{k})},
  \end{equation*}
  where $x_i$ is the vector of explanatory variables for observation
  $i$, and $\beta_j$ is the vector of coefficients for category $j$.
\end{itemize}

\subsubsection{Quantities of Interest}

\begin{itemize}
\item The expected value ({\tt qi\$ev}) is the predicted probability
  for each category:
\begin{equation*}
  E(Y) \; = \; \pi_{ij}\; = \; \frac{\exp(x_{i}\beta_{j})}{\sum^{J}_{k = 1}
      \exp(x_{i}\beta_{k})}.
\end{equation*}

\item The predicted value ({\tt qi\$pr}) is a draw from the
  multinomial distribution defined by the predicted probabilities.  

\item The first difference in predicted
  probabilities ({\tt qi\$fd}), for each category is given by:
\begin{equation*}
\textrm{FD}_j = \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 ({\tt att.ev}) for the treatment group is 
    \begin{equation*} \frac{1}{n_j}\sum_{i:t_i=1}^{n_j} \left\{ Y_i(t_i=1) -
      E[Y_i(t_i=0)] \right\},
    \end{equation*} 
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 ({\tt att.pr}) for the treatment group is 
    \begin{equation*} \frac{1}{n_j}\sum_{i:t_i=1}^{n_j} \left\{ Y_i(t_i=1) -
      \widehat{Y_i(t_i=0)} \right\},
    \end{equation*} 
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 = "mlogit", 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 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/mlogit}
\input{citeZelig}
\subsection*{See also}
The multinomial logit 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{KinTomWit00}.

\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: 

