\documentclass[article,shortnames,nojss]{jss}
\usepackage{thumbpdf}
%% need no \usepackage{Sweave.sty}

\SweaveOpts{engine=R,eps=FALSE}
%\VignetteIndexEntry{The VGAM Package for Categorical Data Analysis}
%\VignetteDepends{VGAM}
%\VignetteKeywords{categorical data analysis, Fisher scoring, iteratively reweighted least squares, multinomial distribution, nominal and ordinal polytomous responses, smoothing, vector generalized linear and additive models, VGAM R package}
%\VignettePackage{VGAM}

%% new commands
\newcommand{\sVLM}{\mbox{\scriptsize VLM}}
\newcommand{\sformtwo}{\mbox{\scriptsize F2}}
\newcommand{\pr}{\mbox{$P$}}
\newcommand{\logit}{\mbox{\rm logit}}
\newcommand{\bzero}{{\bf 0}}
\newcommand{\bone}{{\bf 1}}
\newcommand{\bid}{\mbox{\boldmath $d$}}
\newcommand{\bie}{\mbox{\boldmath $e$}}
\newcommand{\bif}{\mbox{\boldmath $f$}}
\newcommand{\bix}{\mbox{\boldmath $x$}}
\newcommand{\biy}{\mbox{\boldmath $y$}}
\newcommand{\biz}{\mbox{\boldmath $z$}}
\newcommand{\biY}{\mbox{\boldmath $Y$}}
\newcommand{\bA}{\mbox{\rm \bf A}}
\newcommand{\bB}{\mbox{\rm \bf B}}
\newcommand{\bC}{\mbox{\rm \bf C}}
\newcommand{\bH}{\mbox{\rm \bf H}}
\newcommand{\bI}{\mbox{\rm \bf I}}
\newcommand{\bX}{\mbox{\rm \bf X}}
\newcommand{\bW}{\mbox{\rm \bf W}}
\newcommand{\bY}{\mbox{\rm \bf Y}}
\newcommand{\bbeta}{\mbox{\boldmath $\beta$}}
\newcommand{\boldeta}{\mbox{\boldmath $\eta$}}
\newcommand{\bmu}{\mbox{\boldmath $\mu$}}
\newcommand{\bnu}{\mbox{\boldmath $\nu$}}
\newcommand{\diag}{ \mbox{\rm diag} }
\newcommand{\Var}{ \mbox{\rm Var} }
\newcommand{\R}{{\textsf{R}}}
\newcommand{\VGAM}{\pkg{VGAM}}


\author{Thomas W.~Yee\\University of Auckland}
\Plainauthor{Thomas W. Yee}

\title{The \pkg{VGAM} Package for Categorical Data Analysis}
\Plaintitle{The VGAM Package for Categorical Data Analysis}

\Abstract{
  Classical categorical regression models such as the multinomial logit and
  proportional odds models are shown to be readily handled by the  vector
  generalized linear and additive model (VGLM/VGAM) framework. Additionally,
  there are natural extensions, such as reduced-rank VGLMs for
  dimension reduction, and allowing covariates that have values
  specific to each linear/additive predictor,
  e.g., for consumer choice modeling. This article describes some of the
  framework behind the \pkg{VGAM} \R{}~package, its usage and implementation
  details.
}
\Keywords{categorical data analysis, Fisher scoring,
  iteratively reweighted least squares,
  multinomial distribution, nominal and ordinal polytomous responses,
  smoothing, vector generalized linear and additive models,
  \VGAM{} \R{} package}
\Plainkeywords{categorical data analysis, Fisher scoring,
  iteratively reweighted least squares, multinomial distribution,
  nominal and ordinal polytomous responses, smoothing,
  vector generalized linear and additive models, VGAM R package}

\Address{
  Thomas W. Yee \\
  Department of Statistics \\
  University of Auckland, Private Bag 92019 \\
  Auckland Mail Centre \\
  Auckland 1142, New Zealand \\
  E-mail: \email{t.yee@auckland.ac.nz}\\
  URL: \url{http://www.stat.auckland.ac.nz/~yee/}
}


\begin{document}


<<echo=FALSE, results=hide>>=
library("VGAM")
ps.options(pointsize = 12)
options(width = 72, digits = 4)
options(SweaveHooks = list(fig = function() par(las = 1)))
options(prompt = "R> ", continue = "+")
@


% ----------------------------------------------------------------------
\section{Introduction}
\label{sec:jsscat.intoduction}


This is a \pkg{VGAM} vignette for categorical data analysis (CDA)
based on~\cite{Yee:2010}.
Any subsequent features (especially non-backward compatible ones)
will appear here.

The subject of CDA is concerned with
analyses where the response is categorical regardless of whether
the explanatory variables are continuous or categorical. It is a
very frequent form of data. Over the years several CDA regression
models for polytomous responses have become popular, e.g., those
in Table~\ref{tab:cat.quantities}. Not surprisingly, the models
are interrelated: their foundation is the multinomial distribution
and consequently they share similar and overlapping properties which
modellers should know and exploit. Unfortunately, software has been
slow to reflect their commonality and this makes analyses unnecessarily
difficult for the practitioner on several fronts, e.g., using different
functions/procedures to fit different models which does not aid the
understanding of their connections.


This historical misfortune can be seen by considering \R{}~functions
for~CDA. From the Comprehensive \proglang{R} Archive Network
(CRAN, \url{http://CRAN.R-project.org/}) there is~\texttt{polr()}
\citep[in \pkg{MASS};][]{Venables+Ripley:2002} for a proportional odds
model and~\texttt{multinom()}
\citep[in~\pkg{nnet};][]{Venables+Ripley:2002} for the multinomial
logit model. However, both of these can be considered `one-off'
modeling functions rather than providing a unified offering for CDA.
The function \texttt{lrm()} \citep[in \pkg{rms};][]{Harrell:2009}
has greater functionality: it can fit the proportional odds model
(and the forward continuation ratio model upon preprocessing). Neither
\texttt{polr()} or \texttt{lrm()} appear able to fit the nonproportional
odds model. There are non-CRAN packages too, such as the modeling
function~\texttt{nordr()} \citep[in \pkg{gnlm};][]{gnlm:2007}, which can fit
the proportional odds, continuation ratio and adjacent categories models;
however it calls \texttt{nlm()} and the user must supply starting values.
In general these \R{} \citep{R} modeling functions are not modular
and often require preprocessing and sometimes are not self-starting.
The implementations can be perceived as a smattering and piecemeal
in nature. Consequently if the practitioner wishes to fit the models
of Table~\ref{tab:cat.quantities} then there is a need to master several
modeling functions from several packages each having different syntaxes
etc. This is a hindrance to efficient CDA.


 
\begin{table}[tt]
\centering
\begin{tabular}{|c|c|l|}
\hline
Quantity & Notation &
%Range of~$j$ &
\VGAM{} family function \\
\hline
%
$\pr(Y=j+1) / \pr(Y=j)$ &$\zeta_{j}$ &
%$1,\ldots,M$ &
\texttt{acat()} \\
%
$\pr(Y=j) / \pr(Y=j+1)$ &$\zeta_{j}^{R}$ &
%$2,\ldots,M+1$ &
\texttt{acat(reverse = TRUE)} \\
%
$\pr(Y>j|Y \geq j)$ &$\delta_{j}^*$ &
%$1,\ldots,M$ & 
\texttt{cratio()} \\
%
$\pr(Y<j|Y \leq j)$ &$\delta_{j}^{*R}$ &
%$2,\ldots,M+1$ &
\texttt{cratio(reverse = TRUE)} \\
%
$\pr(Y\leq j)$ &$\gamma_{j}$ &
%$1,\ldots,M$ &
\texttt{cumulative()} \\
%
$\pr(Y\geq j)$ &$\gamma_{j}^R$&
%$2,\ldots,M+1$ &
\texttt{cumulative(reverse = TRUE)} \\
%
$\log\{\pr(Y=j)/\pr(Y=M+1)\}$ & &
%$1,\ldots,M$ &
\texttt{multinomial()} \\
%
$\pr(Y=j|Y \geq j)$ &$\delta_{j}$ &
%$1,\ldots,M$ &
\texttt{sratio()} \\
%
$\pr(Y=j|Y \leq j)$ &$\delta_{j}^R$ &
%$2,\ldots,M+1$ &
\texttt{sratio(reverse = TRUE)} \\
%
\hline
\end{tabular}
\caption{
Quantities defined in \VGAM{} for a
categorical response~$Y$ taking values $1,\ldots,M+1$.
Covariates \bix{} have been omitted for clarity.
The LHS quantities are~$\eta_{j}$
or~$\eta_{j-1}$ for~$j=1,\ldots,M$ (not reversed)
and~$j=2,\ldots,M+1$ (if reversed), respectively.
All models are estimated by minimizing the deviance.
All except for \texttt{multinomial()} are suited to ordinal~$Y$.
\label{tab:cat.quantities}
}
\end{table}
 



\proglang{SAS} \citep{SAS} does not fare much better than~\R. Indeed,
it could be considered as having an \textit{excess} of options which
bewilders the non-expert user; there is little coherent overriding
structure. Its \code{proc logistic} handles the multinomial logit
and proportional odds models, as well as exact logistic regression
\citep[see][which is for Version~8 of \proglang{SAS}]{stok:davi:koch:2000}.
The fact that the proportional odds model may be fitted by \code{proc
logistic}, \code{proc genmod} and \code{proc probit} arguably leads
to possible confusion rather than the making of connections, e.g.,
\code{genmod} is primarily for GLMs and the proportional odds model is not
a GLM in the classical \cite{neld:wedd:1972} sense. Also, \code{proc
phreg} fits the multinomial logit model, and \code{proc catmod} with
its WLS implementation adds to further potential confusion.


This article attempts to show how these deficiencies can be addressed
by considering the vector generalized linear and additive model
(VGLM/VGAM) framework, as implemented by the author's~\pkg{VGAM}
package for~\R{}. The main purpose of this paper is to demonstrate
how the framework is very well suited to many `classical' regression
models for categorical responses, and to describe the implementation and
usage of~\pkg{VGAM} for such. To this end an outline of this article
is as follows. Section~\ref{sec:jsscat.VGLMVGAMoverview} summarizes
the basic VGLM/VGAM framework. Section~\ref{sec:jsscat.vgamff}
centers on functions for CDA in~\VGAM. Given an adequate framework,
some natural extensions of Section~\ref{sec:jsscat.VGLMVGAMoverview} are
described in Section~\ref{sec:jsscat.othermodels}. Users of \pkg{VGAM}
can benefit from Section~\ref{sec:jsscat.userTopics} which shows how
the software reflects their common theory. Some examples are given in
Section~\ref{sec:jsscat.eg}. Section~\ref{sec:jsscat.implementDetails}
contains selected topics in statistial computing that are
more relevant to programmers interested in the underlying code.
Section~\ref{sec:jsscat.extnUtil} discusses several utilities and
extensions needed for advanced CDA modeling, and the article concludes
with a discussion. This document was run using \pkg{VGAM}~0.7-10
\citep{yee:VGAM:2010} under \R~2.10.0.


Some general references for categorical data providing
background to this article include
\cite{agre:2002},
\cite{fahr:tutz:2001},
\cite{leon:2000},
\cite{lloy:1999},
\cite{long:1997},
\cite{mccu:neld:1989} and
\cite{simo:2003}.
An overview of models for ordinal responses is~\cite{liu:agre:2005},
and a manual for fitting common models found in~\cite{agre:2002}
to polytomous responses with various software is~\cite{thom:2009}.
A package for visualizing categorical data in~\R{} is~\pkg{vcd}
\citep{Meyer+Zeileis+Hornik:2006,Meyer+Zeileis+Hornik:2009}.






% ----------------------------------------------------------------------
\section{VGLM/VGAM overview}
\label{sec:jsscat.VGLMVGAMoverview}


This section summarizes the VGLM/VGAM framework with a particular emphasis
toward categorical models since the classes encapsulates many multivariate
response models in, e.g., survival analysis, extreme value analysis,
quantile and expectile regression, time series, bioassay data, nonlinear
least-squares models, and scores of standard and nonstandard univariate
and continuous distributions. The framework is partially summarized by
Table~\ref{tab:rrvglam.jss.subset}. More general details about VGLMs
and VGAMs can be found in \cite{yee:hast:2003} and \cite{yee:wild:1996}
respectively. An informal and practical article connecting the general
framework with the software is~\cite{Rnews:Yee:2008}.



\subsection{VGLMs}
\label{sec:wffc.appendixa.vglms}

Suppose the observed response \biy{} is a $q$-dimensional vector.
VGLMs are defined as a model for which the conditional distribution
of $\biY$ given explanatory $\bix$ is of the form
\begin{eqnarray}
f(\biy | \bix ; \bB, \phi) ~=~ h(\biy, \eta_1,\ldots, \eta_M, \phi)
\label{gammod}
\end{eqnarray}
for some known function $h(\cdot)$, where $\bB = (\bbeta_1 \,
\bbeta_2 \, \cdots \, \bbeta_M)$ is a $p \times M$ matrix of
unknown regression coefficients,
and the~$j$th linear predictor is
\begin{equation}
\eta_j ~=~ \eta_j(\bix) ~=~ \bbeta_j^{\top} \bix ~=~
\sum_{k=1}^p \beta_{(j)k} \, x_k , ~~~~ j=1,\ldots,M.
\label{gammod2}
\end{equation}
Here $\bix=(x_1,\ldots,x_p)^{\top}$ with $x_1 = 1$ if there is an intercept.
Note that~(\ref{gammod2}) means that \textit{all} the parameters may be
potentially modelled as functions of~\bix. It can be seen that VGLMs are
like GLMs but allow for multiple linear predictors, and they encompass
models outside the small confines of the exponential family.
In~(\ref{gammod}) the quantity~$\phi$ is an optional scaling parameter
which is included for backward compatibility with common adjustments
to overdispersion, e.g., with respect to GLMs.


In general there is no relationship between~$q$ and~$M$: it
depends specifically on the model or distribution to be fitted.
However, for the `classical' categorical regression models of
Table~\ref{tab:cat.quantities} we have~$M=q-1$ since~$q$ is the number
of levels the multi-category response~$Y$ has.





The $\eta_j$ of VGLMs may be applied directly to parameters of a
distribution rather than just to a mean for GLMs. A simple example is
a univariate distribution with a location parameter~$\xi$ and a scale
parameter~$\sigma > 0$, where we may take~$\eta_1 = \xi$ and~$\eta_2 =
\log\,\sigma$. In general, $\eta_{j}=g_{j}(\theta_{j})$ for some parameter
link function~$g_{j}$ and parameter~$\theta_{j}$.
For example, the adjacent categories models in
Table~\ref{tab:cat.quantities} are ratios of two probabilities, therefore
a log link of~$\zeta_{j}^{R}$ or~$\zeta_{j}$ is the default.
In \VGAM{}, there are currently over a dozen links to choose from, of
which any can be assigned to any parameter, ensuring maximum flexibility.
Table~\ref{tab:jsscat.links} lists some of them.



\begin{table}[tt]
\centering
%\ ~~~~ \par
\begin{tabular}{|l|l|l|l|}
\hline
\ \ ~~~~~~~~~~~~ $\boldeta$ & 
Model & Modeling & Reference \\
 & & function & \\
%-------------------------------------------------------------
\hline
\hline
%-------------------------------------------------------------
 &&&\\[-1.1ex]
$\bB_1^{\top} \bix_{1} + \bB_2^{\top} \bix_{2}\ ( = \bB^{\top} \bix)$ &
VGLM & \texttt{vglm()}
&
\cite{yee:hast:2003} \\[1.6ex]
%Yee \& Hastie~(2003) \\[1.6ex]
%-------------------------------------------------------------
\hline
 &&&\\[-1.1ex]
$\bB_1^{\top} \bix_{1} +
 \sum\limits_{k=p_1+1}^{p_1+p_2} \bH_k \, \bif_{k}^{*}(x_k)$ &
%\sum\limits_{k=1}^{p_2} \bH_k \, \bif_k(x_k)$ &
VGAM & \texttt{vgam()}
&
\cite{yee:wild:1996} \\[2.2ex]
%Yee \& Wild~(1996) \\[2.2ex]
%-------------------------------------------------------------
\hline
 &&&\\[-1.1ex]
$\bB_1^{\top} \bix_{1} + \bA \, \bnu$ &
RR-VGLM & \texttt{rrvglm()}
&
\cite{yee:hast:2003} \\[1.8ex]
%Yee \& Hastie~(2003) \\[1.8ex]
%-------------------------------------------------------------
\hline
 &&&\\[-1.1ex]
See \cite{yee:hast:2003} &
Goodman's~RC & \texttt{grc()}
&
%\cite{yee:hast:2003} \\[1.8ex]
\cite{good:1981} \\[1.8ex]
%-------------------------------------------------------------
\hline
\end{tabular}
\caption{
Some of 
the package \VGAM{} and
its framework.
The vector of latent variables $\bnu = \bC^{\top} \bix_2$
where
$\bix^{\top} = (\bix_1^{\top}, \bix_2^{\top})$.
\label{tab:rrvglam.jss.subset}
}
%\medskip
\end{table}






VGLMs are estimated using iteratively reweighted least squares~(IRLS) 
which is particularly suitable for categorical models
\citep{gree:1984}.
All models in this article have a log-likelihood
\begin{equation}
\ell ~=~ \sum_{i=1}^n \, w_i \, \ell_i
\label{eq:log-likelihood.VGAM}
\end{equation}
where the~$w_i$ are known positive prior weights.
Let~$\bix_i$ denote the explanatory vector for the~$i$th observation,
for $i=1,\dots,n$.
Then one can write
\begin{eqnarray}
\boldeta_i &=& \boldeta(\bix_i) ~=~
\left(
\begin{array}{c}
\eta_1(\bix_i) \\
\vdots \\
\eta_M(\bix_i)
\end{array} \right) ~=~
\bB^{\top} \bix_i ~=~ 
\left(
\begin{array}{c}
\bbeta_1^{\top} \bix_i \\
\vdots \\
\bbeta_M^{\top} \bix_i
\end{array} \right)
\nonumber
\\
&=& 
\left(
\begin{array}{cccc}
\beta_{(1)1} & \cdots & \beta_{(1)p} \\
\vdots \\
\beta_{(M)1} & \cdots & \beta_{(M)p} \\
\end{array} \right)
\bix_i ~=~
\left(
\bbeta_{(1)} \; \cdots \; \bbeta_{(p)}
\right)
\bix_i .
\label{eq:lin.pred}
\end{eqnarray}
In IRLS,
an adjusted dependent vector $\biz_i = \boldeta_i + \bW_i^{-1} \bid_i$
is regressed upon a large (VLM) model matrix, with
$\bid_i = w_i \, \partial \ell_i / \partial \boldeta_i$.
The working weights $\bW_i$ here are 
$w_i \Var(\partial \ell_i / \partial \boldeta_i)$
(which, under regularity conditions, is equal to
$-w_i \, E[ \partial^2 \ell_i / (\partial \boldeta_i \,
\partial \boldeta_i^{\top})]$),
giving rise to the Fisher~scoring algorithm.


Let $\bX=(\bix_1,\ldots,\bix_n)^{\top}$ be the usual $n \times p$
(LM) model matrix
obtained from the \texttt{formula} argument of \texttt{vglm()}.
Given $\biz_i$, $\bW_i$ and~$\bX{}$ at the current IRLS iteration,
a weighted multivariate regression is performed.
To do this, a \textit{vector linear model} (VLM) model matrix 
$\bX_{\sVLM}$ is formed from~$\bX{}$ and~$\bH_k$
(see Section~\ref{sec:wffc.appendixa.vgams}).
This is has $nM$~rows, and if there are no constraints then $Mp$~columns.
Then $\left(\biz_1^{\top},\ldots,\biz_n^{\top}\right)^{\top}$ is regressed
upon $\bX_{\sVLM}$
with variance-covariance matrix $\diag(\bW_1^{-1},\ldots,\bW_n^{-1})$.
This system of linear equations is converted to one large
WLS fit by premultiplication of the output of
a Cholesky decomposition of the~$\bW_i$.


Fisher~scoring usually has good numerical stability
because the~$\bW_i$ are positive-definite over a larger
region of parameter space than Newton-Raphson. 
For the categorical models in this article the expected
information matrices are simpler than the observed
information matrices, and are easily derived,
therefore all the families in Table~\ref{tab:cat.quantities}
implement Fisher~scoring.



\subsection{VGAMs and constraint matrices}
\label{sec:wffc.appendixa.vgams}


VGAMs provide additive-model extensions to VGLMs, that is,
(\ref{gammod2})~is generalized to
\begin{equation}
\eta_j(\bix) ~=~ \beta_{(j)1} +
\sum_{k=2}^p \; f_{(j)k}(x_k), ~~~~ j = 1,\ldots, M,
\label{addmod}
\end{equation}
a sum of smooth functions of the individual covariates, just as
with ordinary GAMs \citep{hast:tibs:1990}. The $\bif_k =
(f_{(1)k}(x_k),\ldots,f_{(M)k}(x_k))^{\top}$ are centered for uniqueness,
and are estimated simultaneously using \textit{vector smoothers}.
VGAMs are thus a visual data-driven method that is well suited to
exploring data, and they retain the simplicity of interpretation that
GAMs possess.



An important concept, especially for CDA, is the idea of
`constraints-on-the functions'.
In practice we often wish to constrain the effect of a covariate to
be the same for some of the~$\eta_j$ and to have no effect for others.
We shall see below that this constraints idea is important
for several categorical models because of a popular parallelism assumption.
As a specific example, for VGAMs we may wish to take
\begin{eqnarray*}
\eta_1 & = & \beta_{(1)1} + f_{(1)2}(x_2) + f_{(1)3}(x_3), \\
\eta_2 & = & \beta_{(2)1} + f_{(1)2}(x_2),
\end{eqnarray*}
so that $f_{(1)2} \equiv f_{(2)2}$ and $f_{(2)3} \equiv 0$.
For VGAMs, we can represent these models using
\begin{eqnarray}
\boldeta(\bix) & = & \bbeta_{(1)} + \sum_{k=2}^p \, \bif_k(x_k)
\ =\ \bH_1 \, \bbeta_{(1)}^* + \sum_{k=2}^p \, \bH_k \, \bif_k^*(x_k)
\label{eqn:constraints.VGAM}
\end{eqnarray}
where $\bH_1,\bH_2,\ldots,\bH_p$ are known full-column rank
\textit{constraint matrices}, $\bif_k^*$ is a vector containing a
possibly reduced set of component functions and $\bbeta_{(1)}^*$ is a
vector of unknown intercepts. With no constraints at all, $\bH_1 =
\bH_2 = \cdots = \bH_p = \bI_M$ and $\bbeta_{(1)}^* = \bbeta_{(1)}$.
Like the $\bif_k$, the~$\bif_k^*$ are centered for uniqueness.
For VGLMs, the~$\bif_k$ are linear so that
\begin{eqnarray}
{\bB}^{\top} &=&
\left(
\bH_1 \bbeta_{(1)}^*
 \;
\Bigg|
 \;
\bH_2 \bbeta_{(2)}^*
 \;
\Bigg|
 \;
\cdots
 \;
\Bigg|
 \;
\bH_p \bbeta_{(p)}^*
\right) 
\label{eqn:lin.coefs4}
\end{eqnarray}
for some vectors
$\bbeta_{(1)}^*,\ldots,\bbeta_{(p)}^*$.


The
$\bX_{\sVLM}$ matrix is constructed from \bX{} and the $\bH_k$ using
Kronecker product operations.
For example, with trivial constraints,
$\bX_{\sVLM} = \bX \otimes \bI_M$.
More generally,
\begin{eqnarray}
\bX_{\sVLM} &=& 
\left(
\left( \bX \, \bie_{1} \right) \otimes \bH_1
 \;
\Bigg|
 \;
\left( \bX \, \bie_{2} \right) \otimes \bH_2
 \;
\Bigg|
 \;
\cdots
 \;
\Bigg|
 \;
\left( \bX \, \bie_{p} \right) \otimes \bH_p
\right)
\label{eqn:X_vlm_Hk}
\end{eqnarray}
($\bie_{k}$ is a vector of zeros except for a one in the $k$th~position)
so that 
$\bX_{\sVLM}$ is $(nM) \times p^*$ where
$p^* = \sum_{k=1}^{p} \mbox{\textrm{ncol}}(\bH_k)$ is the total number
of columns of all the constraint matrices.
Note that $\bX_{\sVLM}$ and \bX{} can be obtained by
\texttt{model.matrix(vglmObject, type = "vlm")}
and
\texttt{model.matrix(vglmObject, type = "lm")}
respectively.
Equation~\ref{eqn:lin.coefs4} focusses on the rows of~\bB{} whereas
\ref{eq:lin.pred}~is on the columns.


VGAMs are estimated by applying a modified vector backfitting algorithm
\citep[cf.][]{buja:hast:tibs:1989} to the $\biz_i$.



\subsection{Vector splines and penalized likelihood}
\label{sec:ex.vspline}

If~(\ref{eqn:constraints.VGAM}) is estimated using a vector spline (a
natural extension of the cubic smoothing spline to vector responses)
then it can be shown that the resulting solution maximizes a penalized
likelihood; some details are sketched in~\cite{yee:step:2007}. In fact,
knot selection for vector spline follows the same idea as O-splines
\citep[see][]{wand:orme:2008} in order to lower the computational cost.


The usage of \texttt{vgam()} with smoothing is very similar
to~\texttt{gam()} \citep{gam:pack:2009}, e.g.,
to fit a nonparametric proportional odds model
\citep[cf.~p.179 of][]{mccu:neld:1989}
to the pneumoconiosis data one could try
<<eval=T>>=
pneumo <- transform(pneumo, let = log(exposure.time))
fit <- vgam(cbind(normal, mild, severe) ~ s(let, df = 2),
            cumulative(reverse = TRUE, parallel = TRUE), pneumo)
@
Here, setting \texttt{df = 1} means a linear fit so that
\texttt{df = 2} affords a little nonlinearity.






% ----------------------------------------------------------------------
\section[VGAM family functions]{\pkg{VGAM} family functions}
\label{sec:jsscat.vgamff}



This section summarizes and comments on the~\VGAM{} family functions
of Table~\ref{tab:cat.quantities} for a categorical response variable
taking values $Y=1,2,\ldots,M+1$. In its most basic invokation, the usage
entails a trivial change compared to \texttt{glm()}: use \texttt{vglm()}
instead and assign the \texttt{family} argument a \VGAM{}~family function.
The use of a \VGAM{}~family function to fit a specific model is far
simpler than having a different modeling function for each model.
Options specific to that model appear as arguments of that \VGAM{}~family
function.





While writing \texttt{cratio()} it was found that various authors defined
the quantity ``continuation ratio'' differently, therefore it became
necessary to define a ``stopping ratio''. Table~\ref{tab:cat.quantities}
defines these quantities for \VGAM{}.




The multinomial logit model is usually described by choosing the first or
last level of the factor to be baseline. \VGAM{}~chooses the last level
(Table~\ref{tab:cat.quantities}) by default, however that can be changed
to any other level by use of the \texttt{refLevel} argument.




If the proportional odds assumption is inadequate then one strategy is
to try use a different link function (see Section~\ref{sec:jsscat.links}
for a selection). Another alternative is to add extra terms such as
interaction terms into the linear predictor
\citep[available in the \proglang{S}~language;][]{cham:hast:1993}.
Another is to fit the so-called \textit{partial}
proportional odds model \citep{pete:harr:1990}
which \VGAM{} can fit via constraint matrices.



In the terminology of~\cite{agre:2002},
\texttt{cumulative()} fits the class of \textit{cumulative link models},
e.g.,
\texttt{cumulative(link = probit)} is a cumulative probit model.
For \texttt{cumulative()}
it was difficult to decide whether
\texttt{parallel = TRUE}
or
\texttt{parallel = FALSE}
should be the default.
In fact, the latter is (for now?).
Users need to set
\texttt{cumulative(parallel = TRUE)} explicitly to
fit a proportional odds model---hopefully this will alert
them to the fact that they are making
the proportional odds assumption and
check its validity (\cite{pete:1990}; e.g., through a deviance or
likelihood ratio test). However the default means numerical problems
can occur with far greater likelihood.
Thus there is tension between the two options.
As a compromise there is now a \VGAM{} family function
called \texttt{propodds(reverse = TRUE)} which is equivalent to
\texttt{cumulative(parallel = TRUE, reverse = reverse, link = "logit")}.



By the way, note that arguments such as 
\texttt{parallel}
can handle a slightly more complex syntax.
A call such as
\code{parallel = TRUE ~ x2 + x5 - 1} means the parallelism assumption
is only applied to~$X_2$ and~$X_5$.
This might be equivalent to something like
\code{parallel = FALSE ~ x3 + x4}, i.e., to the remaining
explanatory variables.





% ----------------------------------------------------------------------
\section{Other models}
\label{sec:jsscat.othermodels}


Given the VGLM/VGAM framework of Section~\ref{sec:jsscat.VGLMVGAMoverview}
it is found that natural extensions are readily proposed in several
directions. This section describes some such extensions.




\subsection{Reduced-rank VGLMs}
\label{sec:jsscat.RRVGLMs}


Consider a multinomial logit model where $p$ and $M$ are both large.
A (not-too-convincing) example might be the data frame \texttt{vowel.test}
in the package \pkg{ElemStatLearn} \citep[see][]{hast:tibs:buja:1994}.
The vowel recognition data set involves $q=11$ symbols produced from
8~speakers with 6~replications of each. The training data comprises
$10$~input features (not including the intercept) based on digitized
utterances. A multinomial logit model fitted to these data would
have $\widehat{\bB}$ comprising of $p \times (q-1) = 110$ regression
coefficients for $n=8\times 6\times 11 = 528$ observations. The ratio
of $n$~to the number of parameters is small, and it would be good to
introduce some parsimony into the model.



A simple and elegant solution is to represent $\widehat{\bB}$ by
its reduced-rank approximation. To do this, partition $\bix$ into
$(\bix_1^{\top}, \bix_2^{\top})^{\top}$ and $\bB = (\bB_1^{\top} \;
\bB_2^{\top})^{\top}$ so that the reduced-rank regression is applied
to~$\bix_2$. In general, \bB{} is a dense matrix of full rank, i.e., rank
$=\min(M,p)$, and since there are $M \times p$ regression coefficients
to estimate this is `too' large for some models and/or data sets.
If we approximate~$\bB_2$ by a reduced-rank regression \begin{equation}
\label{eq:rrr.BAC} \bB_2 ~=~ \bC{} \, \bA^{\top} \end{equation} and if
the rank~$R$ is kept low then this can cut down the number of regression
coefficients dramatically. If~$R=2$ then the results may be biplotted
(\texttt{biplot()} in \VGAM{}). Here, \bC{} and \bA{} are $p_2 \times R$
and $M \times R$ respectively, and usually they are `thin'.


More generally, the class of \textit{reduced-rank VGLMs} (RR-VGLMs)
is simply a VGLM where~$\bB_2$ is expressed as a product of two thin
estimated matrices (Table~\ref{tab:rrvglam.jss.subset}). Indeed,
\cite{yee:hast:2003} show that RR-VGLMs are VGLMs with constraint
matrices that are unknown and estimated. Computationally, this is
done using an alternating method: in~(\ref{eq:rrr.BAC}) estimate~\bA{}
given the current estimate of~\bC{}, and then estimate~\bC{} given the
current estimate of~\bA{}. This alternating algorithm is repeated until
convergence within each IRLS iteration.


Incidentally, special cases of RR-VGLMs have appeared in the
literature. For example, a RR-multinomial logit model, is known as the
\textit{stereotype} model \citep{ande:1984}. Another is \cite{good:1981}'s
RC~model (see Section~\ref{sec:jsscat.rrr.goodman}) which is reduced-rank
multivariate Poisson model. Note that the parallelism assumption of the
proportional odds model \citep{mccu:neld:1989} can be thought of as a
type of reduced-rank regression where the constraint matrices are thin
($\bone_M$, actually) and known.



The modeling function \texttt{rrvglm()} should work with any \VGAM{}
family function compatible with \texttt{vglm()}. Of course, its
applicability should be restricted to models where a reduced-rank
regression of~$\bB_2$ makes sense.









\subsection[Goodman's R x C association model]{Goodman's $R \times C$ association model}
\label{sec:jsscat.rrr.goodman}





Let~$\bY = [(y_{ij})]$ be a $n \times M$ matrix of counts.
Section~4.2 of~\cite{yee:hast:2003} shows that Goodman's~RC$(R)$ association
model \citep{good:1981} fits within the VGLM framework by setting up
the appropriate indicator variables, structural zeros and constraint
matrices. Goodman's model fits a reduced-rank type model to~\bY{}
by firstly assuming that~$Y_{ij}$ has a Poisson distribution, and that
\begin{eqnarray}
\log \, \mu_{ij} &=& \mu + \alpha_{i} + \gamma_{j} + 
\sum_{k=1}^R a_{ik} \, c_{jk} , 
\ \ \ i=1,\ldots,n;\ \ j=1,\ldots,M,
\label{eqn:goodmanrc}
\end{eqnarray}
where $\mu_{ij} = E(Y_{ij})$ is the mean of the $i$-$j$ cell, and the
rank~$R$ satisfies $R < \min(n,M)$.


The modeling function \texttt{grc()} should work on any two-way
table~\bY{} of counts generated by~(\ref{eqn:goodmanrc}) provided
the number of 0's is not too large. Its usage is quite simple, e.g.,
\texttt{grc(Ymatrix, Rank = 2)} fits a rank-2 model to a matrix of counts.
By default a \texttt{Rank = 1} model is fitted.




\subsection{Bradley-Terry models}
\label{sec:jsscat.brat}

Consider
an experiment consists of $n_{ij}$ judges who compare
pairs of items $T_i$, $i=1,\ldots,M+1$.
They express their preferences between $T_i$ and $T_j$. 
Let $N=\sum \sum_{i<j} n_{ij}$ be the total number of pairwise
comparisons, and assume independence for ratings of the same pair
by different judges and for ratings of different pairs by the same judge.
Let $\pi_i$ be the \textit{worth} of item~$T_i$,
$$
\pr(T_i > T_j) ~=~ p_{i/ij} ~=~ \frac{\pi_i}{\pi_i + \pi_j},
\ ~~~~~i \neq {j},
$$
where ``$T_i>T_j$'' means~$i$ is preferred over~$j$.
Suppose that $\pi_i > 0$.
Let~$Y_{ij}$ be the number of times that $T_i$ is preferred
over~$T_j$ in the~$n_{ij}$ comparisons of the pairs.
Then~$Y_{ij} \sim {\rm Bin}(n_{ij},p_{i/ij})$.
This is a Bradley-Terry model (without ties),
and the \VGAM{} family function is~\texttt{brat()}.


Maximum likelihood estimation of the parameters $\pi_1,\ldots,\pi_{M+1}$
involves maximizing
$$
\prod_{i<j}^{M+1}
\left(
\begin{array}{c}
n_{ij} \\
y_{ij}
\end{array} \right)
\left(
\frac{\pi_i}{\pi_i + \pi_j}
\right)^{y_{ij}}
\left(
\frac{\pi_j}{\pi_i + \pi_j}
\right)^{n_{ij}-y_{ij}} .
$$
By default, $\pi_{M+1} \equiv 1$ is used for identifiability,
however, this can be changed very easily.
Note that one can define 
linear predictors $\eta_{ij}$ of the form
\begin{equation}
\label{eq:bradter.logit}
\logit 
\left(
\frac{\pi_i}{\pi_i + \pi_j}
\right) ~=~ \log 
\left(
\frac{\pi_i}{\pi_j}
\right) ~=~ \lambda_i - \lambda_j .
\end{equation}
The VGAM{} framework can handle the Bradley-Terry model only for
intercept-only models; it has
\begin{equation}
\label{eq:bradter}
\lambda_j ~=~ \eta_j ~=~ \log\, \pi_j = \beta_{(1)j},
\ \ \ \ j=1,\ldots,M.
\end{equation}


As well as having many applications in the field of preferences,
the Bradley-Terry model has many uses in modeling `contests' between
teams~$i$ and~$j$, where only one of the teams can win in each
contest (ties are not allowed under the classical model).
The {packaging} function \texttt{Brat()} can be used to
convert a square matrix into one that has more columns, to
serve as input to \texttt{vglm()}.
For example,
for journal citation data where a citation of article~B
by article~A is a win for article~B and a loss for article~A.
On a specific data set,
<<>>=
journal <- c("Biometrika", "Comm.Statist", "JASA", "JRSS-B")
squaremat <- matrix(c(NA, 33, 320, 284,   730, NA, 813, 276,
                      498, 68, NA, 325,   221, 17, 142, NA), 4, 4)
dimnames(squaremat) <- list(winner = journal, loser = journal)
@
then \texttt{Brat(squaremat)} returns a~$1 \times 12$ matrix.







\subsubsection{Bradley-Terry model with ties}
\label{sec:cat.bratt}


The \VGAM{} family function \texttt{bratt()}
implements
a Bradley-Terry model with ties (no preference), e.g.,
where both $T_i$ and $T_j$ are equally good or bad.
Here we assume
\begin{eqnarray*}
 \pr(T_i > T_j) &=& \frac{\pi_i}{\pi_i + \pi_j + \pi_0},
\ ~~~~~
 \pr(T_i = T_j) ~=~ \frac{\pi_0}{\pi_i + \pi_j + \pi_0},
\end{eqnarray*}
with $\pi_0 > 0$ as an extra parameter.
It has 
$$
\boldeta=(\log \pi_1,\ldots, \log \pi_{M-1}, \log \pi_{0})^{\top}
$$
by default, where there are $M$~competitors and $\pi_M \equiv 1$.
Like \texttt{brat()}, one can choose a different reference group
and reference value.


Other \R{}~packages for the Bradley-Terry model
include \pkg{BradleyTerry2}
by H.~Turner and D.~Firth
\citep[with and without ties;][]{firth:2005,firth:2008}
and \pkg{prefmod} \citep{Hatzinger:2009}.




\begin{table}[tt]
\centering
\begin{tabular}[small]{|l|c|}
\hline
\pkg{VGAM} family function & Independent parameters \\
\hline
\texttt{ABO()} & $p, q$ \\
\texttt{MNSs()} & $m_S, m_s, n_S$ \\
\texttt{AB.Ab.aB.ab()} & $p$ \\
\texttt{AB.Ab.aB.ab2()} & $p$ \\
\texttt{AA.Aa.aa()} & $p_A$ \\
\texttt{G1G2G3()} & $p_1, p_2, f$ \\
\hline
\end{tabular}
\caption{Some genetic models currently implemented
and their unique parameters.
\label{tab:gen.all}
}
\end{table}





\subsection{Genetic models}
\label{sec:jsscat.genetic}


There are quite a number of population genetic models based on the
multinomial distribution,
e.g., \cite{weir:1996}, \cite{lang:2002}.
Table~\ref{tab:gen.all} lists some \pkg{VGAM}~family functions for such.




For example the ABO blood group system
has two independent parameters~$p$ and~$q$, say.
Here,
the blood groups A, B and O~form six possible combinations (genotypes)
consisting of AA, AO, BB, BO, AB, OO
(see Table~\ref{tab:ABO}). A and~B are dominant over
bloodtype~O. Let $p$, $q$ and $r$ be the probabilities
for A, B and~O respectively (so that
$p+q+r=1$) for a given population. 
The log-likelihood function is 
\[
\ell(p,q) \;=\; n_A\, \log(p^2 + 2pr) + n_B\, \log(q^2 + 2qr) + n_{AB}\,
\log(2pq) + 2 n_O\, \log(1-p-q),
\]
where $r = 1 - p -q$, $p \in (\,0,1\,)$,
$q \in (\,0,1\,)$, $p+q<1$.
We let $\boldeta = (g(p), g(r))^{\top}$ where $g$ is the link function.
Any~$g$ from Table~\ref{tab:jsscat.links} appropriate for
a parameter $\theta \in (0,1)$ will do.


A toy example where $p=p_A$ and $q=p_B$ is
<<>>=
abodat <- data.frame(A = 725, B = 258, AB = 72, O = 1073)
fit <- vglm(cbind(A, B, AB, O) ~ 1, ABO, abodat)
coef(fit, matrix = TRUE)
Coef(fit) # Estimated pA and pB
@
The function \texttt{Coef()}, which applies only to intercept-only models,
applies to $g_{j}(\theta_{j})=\eta_{j}$
the inverse link function $g_{j}^{-1}$ to~$\widehat{\eta}_{j}$
to give~$\widehat{\theta}_{j}$.







\begin{table}[tt]
% Same as Table 14.1 of E-J, and Table 2.6 of Weir 1996
\begin{center}
\begin{tabular}{|l|cc|cc|c|c|}
\hline
Genotype   & AA  & AO  & BB  &  BO  & AB  &  OO  \\
Probability&$p^2$&$2pr$&$q^2$&$ 2qr$&$2pq$& $r^2$\\
Blood group&  A  &  A  &  B  &  B   &  AB &  O \\
\hline
\end{tabular}
\end{center}
\caption{Probability table for the ABO blood group system.
Note that $p$~and $q$~are the parameters and $r=1-p-q$.
\label{tab:ABO}
}
\end{table}





\subsection{Three main distributions}
\label{sec:jsscat.3maindist}

\cite{agre:2002} discusses three main distributions for categorical
variables: binomial, multinomial, and Poisson
\citep{thom:2009}.
All these are well-represented in the \VGAM{} package,
accompanied by variant forms.
For example,
there is a
\VGAM{} family function named \texttt{mbinomial()}
which implements a 
matched-binomial (suitable for matched case-control studies),
Poisson ordination (useful in ecology for multi-species-environmental data),
negative binomial families,
positive and zero-altered and zero-inflated variants,
and the bivariate odds ratio model
\citep[\texttt{binom2.or()}; see Section~6.5.6 of][]{mccu:neld:1989}.
The latter has an \texttt{exchangeable} argument to allow for an
exchangeable error structure:
\begin{eqnarray}
\bH_1 ~=~
\left( \begin{array}{cc}
1 & 0 \\
1 & 0 \\
0 & 1 \\
\end{array} \right), ~~~~~
\bH_k ~=~
\left( \begin{array}{cc}
1 \\
1 \\
0 \\
\end{array} \right), ~~k=2,\ldots,p,
\label{eqn:blom.exchangeable}
\end{eqnarray}
since, for data $(Y_1,Y_2,\bix)$,
$\logit \, P\!\left( Y_{j} = 1 \Big{|} \bix \right) = 
\eta_{j}$ for ${j}=1,2$, and
$\log \, \psi = \eta_{3}$
where $\psi$~is the odds ratio,
and so $\eta_{1}=\eta_{2}$.
Here, \texttt{binom2.or(zero = 3)} by default meaning $\psi$~is
modelled as an intercept-only
(in general, \texttt{zero} may be assigned an integer vector
such that the value~$j$ means $\eta_{j} = \beta_{(j)1}$,
i.e., the $j$th~linear/additive predictor is an intercept-only).
See the online help for all of these models.















% ----------------------------------------------------------------------
\section{Some user-oriented topics}
\label{sec:jsscat.userTopics}


Making the most of \VGAM{} requires an understanding of the general
VGLM/VGAM framework described Section~\ref{sec:jsscat.VGLMVGAMoverview}.
In this section we connect elements of that framework with the software.
Before doing so it is noted that
a fitted \VGAM{} categorical model has access to the usual
generic functions, e.g.,
\texttt{coef()} for
$\left(\widehat{\bbeta}_{(1)}^{*T},\ldots,\widehat{\bbeta}_{(p)}^{*T}\right)^{\top}$
(see Equation~\ref{eqn:lin.coefs4}),
\texttt{constraints()} for $\bH_k$,
\texttt{deviance()} for $2\left(\ell_{\mathrm{max}} - \ell\right)$,
\texttt{fitted()} for $\widehat{\bmu}_i$,
\texttt{logLik()} for $\ell$,
\texttt{predict()} for $\widehat{\boldeta}_i$,
\texttt{print()},
\texttt{residuals(..., type = "response")} for $\biy_i - \widehat{\bmu}_i$ etc.,
\texttt{summary()},
\texttt{vcov()} for $\widehat{\Var}(\widehat{\bbeta})$,
etc.
The methods function for the extractor function
\texttt{coef()} has an argument \texttt{matrix}
which, when set \texttt{TRUE}, returns~$\widehat{\bB}$
(see Equation~\ref{gammod}) as a $p \times M$ matrix,
and this is particularly useful for confirming that a fit
has made a parallelism assumption.







\subsection{Common arguments}
\label{sec:jsscat.commonArgs}


The structure of the unified framework given in
Section~\ref{sec:jsscat.VGLMVGAMoverview}
appears clearly through
the pool of common arguments
shared by the
\VGAM{} family functions in Table~\ref{tab:cat.quantities}.
In particular,
\texttt{reverse} and
\texttt{parallel}
are prominent with CDA.
These are merely convenient shortcuts for the argument \texttt{constraints},
which accepts a named list of constraint matrices~$\bH_k$.
For example, setting
\texttt{cumulative(parallel = TRUE)} would constrain the coefficients $\beta_{(j)k}$
in~(\ref{gammod2}) to be equal for all $j=1,\ldots,M$,
each separately for $k=2,\ldots,p$.
That is, $\bH_k = \bone_M$.
The argument~\texttt{reverse} determines the `direction' of
the parameter or quantity.

Another argument not so much used with CDA is~\texttt{zero};
this accepts a vector specifying which~$\eta_j$ is to be modelled as
an intercept-only; assigning a \texttt{NULL} means none.








\subsection{Link functions}
\label{sec:jsscat.links}

Almost all \VGAM{} family functions
(one notable exception is \texttt{multinomial()})
allow, in theory, for any link function to be assigned to each~$\eta_j$.
This provides maximum capability.
If so then there is an extra argument to pass in any known parameter
associated with the link function.
For example, \texttt{link = "logoff", earg = list(offset = 1)}
signifies a log link with a unit offset:
$\eta_{j} = \log(\theta_{j} + 1)$ for some parameter~$\theta_{j}\ (> -1)$.
The name \texttt{earg} stands for ``extra argument''.
Table~\ref{tab:jsscat.links} lists some links relevant to categorical data.
While the default gives a reasonable first choice,
users are encouraged to try different links.
For example, fitting a binary regression model
(\texttt{binomialff()}) to the coal miners data set \texttt{coalminers} with
respect to the response wheeze gives a
nonsignificant regression coefficient for $\beta_{(1)3}$~with probit analysis
but not with a logit link when
$\eta = \beta_{(1)1} + \beta_{(1)2} \, \mathrm{age} + \beta_{(1)3} \, \mathrm{age}^2$.
Developers and serious users are encouraged to write and use
new link functions compatible with~\VGAM.






\begin{table*}[tt]
\centering
\medskip
\begin{tabular}{|l|c|c|}
\hline
Link function & $g(\theta)$ & Range of $\theta$ \\
\hline
\texttt{cauchit()} & $\tan(\pi(\theta-\frac12))$ & $(0,1)$ \\
\texttt{cloglog()} & $\log_e\{-\log_e(1 - \theta)\}$ & $(0,1)$ \\
\texttt{fisherz()} & 
$\frac12\,\log_e\{(1 + \theta)/(1 - \theta)\}$ & $(-1,1)$ \\
\texttt{identity()} & $\theta$ & $(-\infty,\infty)$ \\
\texttt{logc()} & $\log_e(1 - \theta)$ & $(-\infty,1)$ \\
\texttt{loge()} & $\log_e(\theta)$ & $(0,\infty)$ \\
\texttt{logit()} & $\log_e(\theta/(1 - \theta))$ & $(0,1)$ \\
\texttt{logoff()} & $\log_e(\theta + A)$ & $(-A,\infty)$ \\
\texttt{probit()} & $\Phi^{-1}(\theta)$ & $(0,1)$ \\
\texttt{rhobit()} & $\log_e\{(1 + \theta)/(1 - \theta)\}$ & $(-1,1)$ \\
\hline
\end{tabular}
\caption{
Some \VGAM{} link functions pertinent to this article.
\label{tab:jsscat.links}
}
\end{table*}









% ----------------------------------------------------------------------
\section{Examples}
\label{sec:jsscat.eg}

This section illustrates CDA modeling on three
data sets in order to give a flavour of what is available in the package.



\subsection{2008 World Fly Fishing Championships}
\label{sec:jsscat.eg.WFFC}

The World Fly Fishing Championships (WFFC)
is a prestigious catch-and-release competition held annually.
In 2008 it was held in New~Zealand during the month of March.
The data was released and appears in~\VGAM{} as the data frames
\texttt{wffc},
\texttt{wffc.nc},
\texttt{wffc.indiv} and
\texttt{wffc.teams}.
Details about the competition are found
in the online help, as well as~\cite{yee:2010v}.


Briefly, we will model the abundance of fish caught during each
three-hour session amongst the 90 or so competitors (from about
19~countries) who fished all their sessions. There were five~sectors
(locations) labelled I--V for the Whanganui~River, Lake~Otamangakau,
Lake~Rotoaira, Waihou~River and Waimakariri~River, respectively. The
sessions were sequentially labelled 1--6 where odd and even
numbers denote morning and afternoon respectively. There were
three consecutive days of fishing during which each sector experienced
a rest session.



\cite{yee:2010v} fitted Poisson and negative
binomial regressions to the numbers caught at each competitor-session
combination.
The negative
binomial regression had an intercept-only for its
index parameter~$k$ and
$\Var(Y) = \mu(1+\mu / k)$.
Both models had the log-linear relationship
\begin{eqnarray}
\label{eq:wffc.use.loglinear}
\log \, \mu_{adsc} &=&
\eta ~=~
\beta_{(1)1} +
\alpha_{s} +
\beta_{a} +
\gamma_{d} +
\delta_{c}.
\end{eqnarray}
where $\mu = E(Y)$ is the mean number caught,
$\beta_{(1)1}$~is the intercept,
$\alpha_{s}$~are the sector effects for $s=1,\ldots,5$ sectors,
$\delta_{c}$~are the ``competitor effects'' for $c=1,\ldots,91$ competitors
(8~competitors who did not fish all~5 sessions were excluded),
$\beta_{a}$~are the morning ($a=1$) and afternoon ($a=2$) effects,
$\gamma_{d}$~are the day effects for
day $d=1,2,3$.
Recall for factors that the first level is baseline, e.g.,
$\alpha_1=\beta_1=0$ etc.
Not used here is $b=1,\ldots,19$ for which beat/boat was
fished/used (e.g., fixed locations on the river).
We will fit a proportional odds model with essentially the RHS
of~(\ref{eq:wffc.use.loglinear}) as the linear predictor.



Here is a peek at the data frame used.
Each row of~\texttt{wffc.nc} is the number of captures
by each sector-session-beat combination.
<<>>=
head(wffc.nc, 5)
@ 


We first process the data a little: create the regressor
variables and restrict the analysis to anglers who fished all their
sessions.
Here, ``\texttt{nc}'' stands for numbers caught, and
``\texttt{f}'' stands for factor.
<<>>=
fnc <- transform(wffc.nc,
                 finame = factor(iname),
                 fsector = factor(sector),
                 fday = factor(ceiling(session / 2)),
                 mornaft = 1 - (session %% 2),
                 fbeatboat = factor(beatboat))

fnc <- fnc[with(fnc, !is.element(comid, c(99,72,80,93,45,71,97,78))),] 
fnc <- transform(fnc,
                ordnum = ifelse(numbers <= 02, "few",
                         ifelse(numbers <= 10, "more", "most")))
fnc$ordnum <- ordered(fnc$ordnum, levels = c("few", "more", "most"))
@
The variable \texttt{ordnum} is ordinal with 3~levels.
The cut-points chosen here were decided upon by manual
inspection; they gave approximately the same numbers in each level:
<<>>=
with(fnc, table(ordnum))
@



Now we are in a position to fit a proportional odds model
to mimic~(\ref{eq:wffc.use.loglinear}).
<<>>=
fit.pom <- vglm(ordnum ~
          fsector +
          mornaft +
          fday +
          finame,
          family = cumulative(parallel = TRUE, reverse = TRUE),
          data = fnc)
@
Here, we set \texttt{reverse = TRUE} so that the coefficients
have the same direction as a logistic regression.
It means that if a regression coefficient is positive then 
an increasing value of an explanatory variable is associated
with an increasing value of the response.
One could have used \texttt{family = propodds} instead.



Before interpreting some output let's check that the input was alright.
<<>>=
head(fit.pom@y, 3)
colSums(fit.pom@y)
@
The checking indicates no problems with the input.



Now let's look at some output. Note that the Whanganui~River,
Mornings and Day~1 are the baseline levels of the factors. Also, the
variable \texttt{mornaft} is 0~for morning and 1~for afternoons.
Likewise, the factor \texttt{fday} has values \texttt{1}, \texttt{2}
and \texttt{3}.
<<>>=
head(coef(fit.pom, matrix = TRUE), 10)
#head(summary(fit.pom)@coef3, 10) # Old now since 0.7-10 is nicer
@
verifies the parallelism assumption.
Standard errors and Wald statistics may be obtained by
<<>>=
head(coef(summary(fit.pom)), 10)
@
Not surprisingly, these results agree with the Poisson
and negative binomial regressions
\citep[reported in][]{yee:2010v}.
The most glaring qualitative results are as follows.
We use the rough rule of thumb that if
the absolute value of the $t$~{statistic} is greater than~$2$
then it is `statistically significant'.
\begin{itemize}

\item
The two lakes were clearly less productive than the rivers.
However, neither of the other two rivers were significantly different
from the Whanganui~River.



\item
There is a noticeable day effect: the second day is not significantly
different from the opening day but it is for the third day.
The decreasing values of the fitted coefficients show there is an
increasing catch-reduction (fish depletion if it were catch-and-keep)
as the competition progressed. Replacing \texttt{fday} by a
variable~\texttt{day} and entering that linearly gave a $t$~statistic
of~$-4.0$: there is a significant decrease in catch over time.



\item
Mornings were more productive than afternoons. The $p$~value for this
would be close to~5\%. This result is in line with the day effect: fishing
often results in a `hammering' effect over time on fish populations,
especially in small streams. Since the morning and afternoon sessions
were fixed at 9.00am--12.00pm and 2.30--5.30pm daily, there was only
$2\frac12$~hours for the fish to recover until the next angler arrived.


\end{itemize}





Let us check the proportional odds assumption with respect
to the variable~\texttt{mornaft}.
<<>>=
fit.ppom <- vglm(ordnum ~
          fsector +
          mornaft +
          fday +
          finame,
          cumulative(parallel = FALSE ~ 1 + mornaft, reverse = TRUE),
          data = fnc)
head(coef(fit.ppom, matrix = TRUE),  8)
@
As expected, all rows but~\texttt{(Intercept)} and~\texttt{mornaft}
are identical due to the parallelism.
Then
<<>>=
pchisq(deviance(fit.pom) - deviance(fit.ppom),
       df = df.residual(fit.pom) - df.residual(fit.ppom), lower.tail=FALSE)
@
gives a likelihood ratio test $p$~value which is non-significant.
Repeating the testing for each variable separately indicates
that the parallelism assumption seems reasonable here except
with~\texttt{fday} ($p$~value $\approx 0.012$).
For this model
<<>>=
fit2.ppom <- vglm(ordnum ~
          fsector +
          mornaft +
          fday +
          finame,
          family = cumulative(parallel = FALSE ~ 1 + fday, reverse = TRUE),
          data = fnc)
head(coef(fit2.ppom, matrix = TRUE), 8)
@


Some miscellaneous output is as follows.
<<>>=
head(fitted(fit2.ppom), 3)
@
are the fitted probabilities $\widehat{P}(Y={j})$ which sum to unity for
each row.
The $i$th~row of
<<>>=
head(predict(fit2.ppom), 3)
@
is $\widehat{\boldeta}(\bix_i)^{\top}$.
The dimensions of the LM and VLM design matrices are
<<>>=
dim(model.matrix(fit2.ppom, type = "lm"))
dim(model.matrix(fit2.ppom, type = "vlm"))
@
which shows the VLM matrix grows quickly with respect to~$M$.
Lastly,
<<>>=
constraints(fit2.ppom)[c(1, 2, 5, 6)]
@
shows some of the constraint matrices,
$\bH_1=\bI_2$
and
$\bH_2=\bH_5=\bH_6=\bone_2$
(see Equations~\ref{eqn:constraints.VGAM}--\ref{eqn:lin.coefs4}).







\subsection{Marital status data}
\label{sec:jsscat.eg.mstatus}

We fit a nonparametric multinomial logit model to data collected from
a self-administered questionnaire administered in a large New Zealand
workforce observational study conducted during 1992--3.
The data were augmented by a second study consisting of retirees.
For homogeneity, this analysis is restricted
to a subset of 6053 European males with no missing values.
The ages ranged between~16 and 88~years.
The data can be considered a reasonable representation of the white
male New Zealand population in the early 1990s, and
are detailed in~\cite{macm:etal:1995} and~\cite{yee:wild:1996}.
We are interested in exploring how $Y=$ marital status varies as a function
of $x_2=$ age. The nominal response~$Y$ has four levels;
in sorted order, they are divorced or separated, married or partnered,
single and widower.
We will write these levels as $Y=1$, $2$, $3$, $4$, respectively,
and will choose the married/partnered (second level) as the reference group
because the other levels emanate directly from it.

Suppose the data is in a data frame called \texttt{marital.nz}
and looks like
<<>>=
head(marital.nz, 4)
summary(marital.nz)
@
We fit the VGAM
<<>>=
fit.ms <- vgam(mstatus ~ s(age, df = 3), multinomial(refLevel = 2),
               data = marital.nz)
@

Once again let's firstly check the input.
<<>>=
head(fit.ms@y, 4)
colSums(fit.ms@y)
@
This seems ok.




Now the estimated component functions $\widehat{f}_{(s)2}(x_2)$
may be plotted with
<<fig=F>>=
# Plot output
mycol <- c("red","darkgreen","blue")
 par(mfrow=c(2,2))
plot(fit.ms, se=TRUE, scale=12,
         lcol=mycol, scol=mycol)

# Plot output overlayed
#par(mfrow=c(1,1))
plot(fit.ms, se=TRUE, scale=12,
         overlay=TRUE,
         llwd=2,
         lcol=mycol, scol=mycol)
@
to produce Figure~\ref{fig:jsscat.eg.mstatus}.
The \texttt{scale} argument is used here to ensure that the $y$-axes have
a common scale---this makes comparisons between the component functions
less susceptible to misinterpretation.
The first three plots are the (centered) $\widehat{f}_{(s)2}(x_2)$ for
$\eta_1$,
$\eta_2$,
$\eta_3$,
where
\begin{eqnarray}
\label{eq:jsscat.eg.nzms.cf}
\eta_{s} ~=~
\log(\pr(Y={t}) / \pr(Y={2})) ~=~
\beta_{(s)1} + f_{(s)2}(x_2),
\end{eqnarray}
$(s,t) = (1,1), (2,3), (3,4)$,
and~$x_2$ is~\texttt{age}.
The last plot are the smooths overlaid to aid comparison.


It may be seen that the $\pm 2$ standard error bands
about the \texttt{Widowed} group is particularly wide at
young ages because of a paucity of data, and
likewise at old ages amongst the \texttt{Single}s.
The $\widehat{f}_{(s)2}(x_2)$ appear as one would expect.
The log relative risk of
being single relative to being married/partnered drops sharply from
ages~16 to~40.
The fitted function for the~\texttt{Widowed} group increases
with~\texttt{age} and looks reasonably linear.
The $\widehat{f}_{(1)2}(x_2)$
suggests a possible maximum around 50~years old---this
could indicate the greatest marital conflict occurs during
the mid-life crisis years!



\setkeys{Gin}{width=0.9\textwidth} % 0.8 is the current default

\begin{figure}[tt]
\begin{center}
<<fig=TRUE,width=8,height=5.6,echo=FALSE>>=
# Plot output
mycol <- c("red","darkgreen","blue")
 par(mfrow=c(2,2))
 par(mar=c(4.2,4.0,1.2,2.2)+0.1)
plot(fit.ms, se=TRUE, scale=12,
         lcol=mycol, scol=mycol)

# Plot output overlaid
#par(mfrow=c(1,1))
plot(fit.ms, se=TRUE, scale=12,
         overlay=TRUE,
         llwd=2,
         lcol=mycol, scol=mycol)
@
\caption{
Fitted (and centered) component functions
$\widehat{f}_{(s)2}(x_2)$
from the NZ marital status data
(see Equation~\ref{eq:jsscat.eg.nzms.cf}).
The bottom RHS plot are the smooths overlaid.
\label{fig:jsscat.eg.mstatus}
}
\end{center}
\end{figure}

\setkeys{Gin}{width=0.8\textwidth} % 0.8 is the current default



The methods function for~\texttt{plot()} can also plot the
derivatives of the smooths.
The call
<<fig=F>>=
plot(fit.ms, deriv=1, lcol=mycol, scale=0.3)
@
results in Figure~\ref{fig:jsscat.eg.mstatus.cf.deriv}.
Once again the $y$-axis scales are commensurate.

\setkeys{Gin}{width=\textwidth} % 0.8 is the current default

\begin{figure}[tt]
\begin{center}
<<fig=TRUE,width=7.2,height=2.4,echo=FALSE>>=
# Plot output
 par(mfrow=c(1,3))
 par(mar=c(4.5,4.0,0.2,2.2)+0.1)
plot(fit.ms, deriv=1, lcol=mycol, scale=0.3)
@
\caption{
Estimated first derivatives of the component functions,
$\widehat{f'}_{(s)2}(x_2)$,
from the NZ marital status data
(see Equation~\ref{eq:jsscat.eg.nzms.cf}).
\label{fig:jsscat.eg.mstatus.cf.deriv}
}
\end{center}
\end{figure}

\setkeys{Gin}{width=0.8\textwidth} % 0.8 is the current default


The derivative for the \texttt{Divorced/Separated} group appears
linear so that a quadratic component function could be tried.
Not surprisingly the \texttt{Single} group shows the greatest change;
also, $\widehat{f'}_{(2)2}(x_2)$ is approximately linear till~50
and then flat---this suggests one could fit a piecewise quadratic
function to model that component function up to 50~years.
The~\texttt{Widowed} group appears largely flat.
We thus fit the parametric model
<<>>=
foo <- function(x, elbow=50)
    poly(pmin(x, elbow), 2)

clist <- list("(Intercept)" = diag(3),
             "poly(age, 2)" = rbind(1, 0, 0),
             "foo(age)" = rbind(0, 1, 0),
             "age" = rbind(0, 0, 1))
fit2.ms <-
    vglm(mstatus ~ poly(age, 2) + foo(age) + age,
         family = multinomial(refLevel = 2),
         constraints = clist,
         data = marital.nz)
@
Then
<<>>=
coef(fit2.ms, matrix = TRUE)
@
confirms that one term was used for each component function.
The plots from
<<fig=F>>=
par(mfrow=c(2,2))
plotvgam(fit2.ms, se = TRUE, scale = 12,
         lcol = mycol[1], scol = mycol[1], which.term = 1)
plotvgam(fit2.ms, se = TRUE, scale = 12,
         lcol = mycol[2], scol=mycol[2], which.term = 2)
plotvgam(fit2.ms, se = TRUE, scale = 12,
         lcol = mycol[3], scol = mycol[3], which.term = 3)
@
are given in Figure~\ref{fig:jsscat.eg.mstatus.vglm}
and appear like
Figure~\ref{fig:jsscat.eg.mstatus}.


\setkeys{Gin}{width=0.9\textwidth} % 0.8 is the current default

\begin{figure}[tt]
\begin{center}
<<fig=TRUE,width=8,height=5.6,echo=FALSE>>=
# Plot output
par(mfrow=c(2,2))
 par(mar=c(4.5,4.0,1.2,2.2)+0.1)
plotvgam(fit2.ms, se = TRUE, scale = 12,
         lcol = mycol[1], scol = mycol[1], which.term = 1)
plotvgam(fit2.ms, se = TRUE, scale = 12,
         lcol = mycol[2], scol = mycol[2], which.term = 2)
plotvgam(fit2.ms, se = TRUE, scale = 12,
         lcol = mycol[3], scol = mycol[3], which.term = 3)
@
\caption{
Parametric version of~\texttt{fit.ms}: \texttt{fit2.ms}.
The component functions are now quadratic, piecewise quadratic/zero,
or linear.
\label{fig:jsscat.eg.mstatus.vglm}
}
\end{center}
\end{figure}

\setkeys{Gin}{width=0.8\textwidth} % 0.8 is the current default




It is possible to perform very crude inference based on heuristic theory
of a deviance test:
<<>>=
deviance(fit.ms) - deviance(fit2.ms)
@
is small, so it seems the parametric model is quite reasonable
against the original nonparametric model.
Specifically,
the difference in the number of `parameters' is approximately
<<>>=
(dfdiff <- df.residual(fit2.ms) - df.residual(fit.ms))
@
which gives an approximate $p$~value of
<<>>=
1-pchisq(deviance(fit.ms) - deviance(fit2.ms), df=dfdiff)
@
Thus \texttt{fit2.ms} appears quite reasonable.








The estimated probabilities of the original fit can be plotted
against~\texttt{age} using
<<fig=F>>=
ooo <- with(marital.nz, order(age))
with(marital.nz, matplot(age[ooo], fitted(fit.ms)[ooo,],
     type="l", las=1, lwd=2, ylim=0:1,
     ylab="Fitted probabilities",
     xlab="Age", # main="Marital status amongst NZ Male Europeans",
     col=c(mycol[1], "black", mycol[-1])))
legend(x=52.5, y=0.62, # x="topright",
       col=c(mycol[1], "black", mycol[-1]),
       lty=1:4,
       legend=colnames(fit.ms@y), lwd=2)
abline(v=seq(10,90,by=5), h=seq(0,1,by=0.1), col="gray", lty="dashed")
@
which gives Figure~\ref{fig:jsscat.eg.mstatus.fitted}.
This shows that between 80--90\%~of NZ white males
aged between their early~30s to mid-70s
were married/partnered.
The proportion widowed
started to rise steeply from 70~years onwards but remained below~0.5
since males die younger than females on average.


\setkeys{Gin}{width=0.8\textwidth} % 0.8 is the current default

\begin{figure}[tt]
\begin{center}
<<fig=TRUE,width=8,height=4.8,echo=FALSE>>=
 par(mfrow=c(1,1))
 par(mar=c(4.5,4.0,0.2,0.2)+0.1)
ooo <- with(marital.nz, order(age))
with(marital.nz, matplot(age[ooo], fitted(fit.ms)[ooo,],
     type="l", las=1, lwd=2, ylim=0:1,
     ylab="Fitted probabilities",
     xlab="Age",
     col=c(mycol[1], "black", mycol[-1])))
legend(x=52.5, y=0.62,
       col=c(mycol[1], "black", mycol[-1]),
       lty=1:4,
       legend=colnames(fit.ms@y), lwd=2.1)
abline(v=seq(10,90,by=5), h=seq(0,1,by=0.1), col="gray", lty="dashed")
@
\caption{
Fitted probabilities for each class for the
NZ male European
marital status data
(from Equation~\ref{eq:jsscat.eg.nzms.cf}).
\label{fig:jsscat.eg.mstatus.fitted}
}
\end{center}
\end{figure}

\setkeys{Gin}{width=0.8\textwidth} % 0.8 is the current default







\subsection{Stereotype model}
\label{sec:jsscat.eg.grc.stereotype}

We reproduce some of the analyses of \cite{ande:1984} regarding the
progress of 101~patients with back pain
using the data frame \texttt{backPain} from \pkg{gnm}
\citep{Rnews:Turner+Firth:2007,Turner+Firth:2009}.
The three prognostic variables are
length of previous attack ($x_1=1,2$),
pain change ($x_2=1,2,3$) 
and lordosis ($x_3=1,2$).
Like him, we treat these as numerical and standardize and negate them.
%
The output
<<>>=
# Scale the variables? Yes; the Anderson (1984) paper did (see his Table 6).
head(backPain, 4)
summary(backPain)
backPain <- transform(backPain, sx1 = -scale(x1), sx2 = -scale(x2), sx3 = -scale(x3))
@
displays the six ordered categories.
Now a rank-1 stereotype model can be fitted with
<<>>=
bp.rrmlm1 <- rrvglm(pain ~ sx1 + sx2 + sx3, multinomial, backPain)
@
Then
<<>>=
Coef(bp.rrmlm1)
@
are the fitted \bA, \bC{} and~$\bB_1$ (see Equation~\ref{eq:rrr.BAC}) and
Table~\ref{tab:rrvglam.jss.subset}) which agrees with his Table~6.
Here, what is known as ``corner constraints'' is used
($(1,1)$ element of \bA{} $\equiv 1$),
and only the intercepts are not subject to any reduced-rank regression
by default.
The maximized log-likelihood from \textsl{\texttt{logLik(bp.rrmlm1)}}
is $\Sexpr{round(logLik(bp.rrmlm1), 2)}$.
The standard errors of each parameter can be obtained by
\textsl{\texttt{summary(bp.rrmlm1)}}.
The negative elements of~$\widehat{\bC}$ imply the
latent variable~$\widehat{\nu}$ decreases in value with increasing
\textsl{\texttt{sx1}},
\textsl{\texttt{sx2}} and
\textsl{\texttt{sx3}}.
The elements of~$\widehat{\bA}$ tend to decrease so it suggests
patients get worse as $\nu$~increases,
i.e., get better as \textsl{\texttt{sx1}},
\textsl{\texttt{sx2}} and
\textsl{\texttt{sx3}} increase.






<<echo=FALSE>>=
set.seed(123)
@
A rank-2 model fitted \textit{with a different normalization}
<<>>=
bp.rrmlm2 <- rrvglm(pain ~ sx1 + sx2 + sx3, multinomial, backPain, Rank = 2,
                   Corner = FALSE, Uncor = TRUE)
@
produces uncorrelated $\widehat{\bnu}_i = \widehat{\bC}^{\top} \bix_{2i}$.
In fact \textsl{\texttt{var(lv(bp.rrmlm2))}} equals $\bI_2$
so that the latent variables are also scaled to have unit variance.
The fit was biplotted
(rows of $\widehat{\bC}$ plotted as arrow;
 rows of $\widehat{\bA}$ plotted as labels) using
<<figure=F>>=
biplot(bp.rrmlm2, Acol="blue", Ccol="darkgreen", scores=TRUE,
#      xlim=c(-1,6), ylim=c(-1.2,4), # Use this if not scaled
       xlim=c(-4.5,2.2), ylim=c(-2.2, 2.2), # Use this if scaled
       chull=TRUE, clty=2, ccol="blue")
@
to give Figure~\ref{fig:jsscat.eg.rrmlm2.backPain}.
It is interpreted via inner products due to~(\ref{eq:rrr.BAC}).
The different normalization means that the interpretation of~$\nu_1$
and~$\nu_2$ has changed, e.g., increasing
\textsl{\texttt{sx1}},
\textsl{\texttt{sx2}} and
\textsl{\texttt{sx3}} results in increasing $\widehat{\nu}_1$ and
patients improve more.
Many of the latent variable points $\widehat{\bnu}_i$ are coincidental
due to discrete nature of the~$\bix_i$. The rows of~$\widehat{\bA}$
are centered on the blue labels (rather cluttered unfortunately) and
do not seem to vary much as a function of~$\nu_2$.
In fact this is confirmed by~\cite{ande:1984} who showed a rank-1
model is to be preferred.



This example demonstrates the ability to obtain a low dimensional view
of higher dimensional data. The package's website has additional
documentation including more detailed Goodman's~RC and stereotype
examples.





\setkeys{Gin}{width=0.8\textwidth} % 0.8 is the current default

\begin{figure}[tt]
\begin{center}
<<fig=TRUE,width=8,height=5.3,echo=FALSE>>=
# Plot output
 par(mfrow=c(1,1))
 par(mar=c(4.5,4.0,0.2,2.2)+0.1)

biplot(bp.rrmlm2, Acol="blue", Ccol="darkgreen", scores=TRUE,
#      xlim=c(-1,6), ylim=c(-1.2,4),  # Use this if not scaled
       xlim=c(-4.5,2.2), ylim=c(-2.2, 2.2),  # Use this if scaled
       chull=TRUE, clty=2, ccol="blue")
@
\caption{
Biplot of a rank-2 reduced-rank multinomial logit (stereotype) model
fitted to the back pain data.
A convex hull surrounds the latent variable scores
$\widehat{\bnu}_i$
(whose observation numbers are obscured because of their discrete nature).
The position of the $j$th~row of~$\widehat{\bA}$
is the center of the label ``\texttt{log(mu[,j])/mu[,6])}''.
\label{fig:jsscat.eg.rrmlm2.backPain}
}
\end{center}
\end{figure}

\setkeys{Gin}{width=0.8\textwidth} % 0.8 is the current default










% ----------------------------------------------------------------------
\section{Some implementation details}
\label{sec:jsscat.implementDetails}

This section describes some implementation details of~\VGAM{}
which will be more of interest to the developer than to the casual user.



\subsection{Common code}
\label{sec:jsscat.implementDetails.code}

It is good programming practice to write reusable code where possible.
All the \VGAM{}~family functions in Table~\ref{tab:cat.quantities}
process the response in the same way because the same segment of code
is executed. This offers a degree of uniformity in terms of how input is
handled, and also for software maintenance
(\cite{altm:jack:2010} enumerates good programming techniques and references).
As well, the default initial values are computed in the same manner
based on sample proportions of each level of~$Y$.





\subsection[Matrix-band format of wz]{Matrix-band format of \texttt{wz}}
\label{sec:jsscat.implementDetails.mbformat}

The working weight matrices $\bW_i$ may become large for categorical
regression models. In general, we have to evaluate the~$\bW_i$
for~$i=1,\ldots,n$, and naively, this could be held in an \texttt{array} of
dimension~\texttt{c(M, M, n)}. However, since the~$\bW_i$ are symmetric
positive-definite it suffices to only store the upper or lower half of
the matrix.



The variable~\texttt{wz} in \texttt{vglm.fit()}
stores the working weight matrices $\bW_i$ in 
a special format called the \textit{matrix-band} format. This
format comprises a $n \times M^*$ matrix where
$$
M^* ~=~ \sum_{i=1}^{\footnotesize \textit{hbw}} \;
\left(M-i+1\right) ~=~ 
\frac12 \, \textit{hbw}\, \left(2\,M - \textit{hbw} +1\right)
$$
is the number of columns. Here, \textit{hbw} refers to the
\textit{half-bandwidth} of the matrix, which is an integer
between~1 and~$M$ inclusive. A diagonal matrix has
unit half-bandwidth, a tridiagonal matrix has half-bandwidth~2, etc.


Suppose $M=4$. Then \texttt{wz} will have up to $M^*=10$ columns
enumerating the unique elements of~$\bW_i$ as follows:
\begin{eqnarray}
\bW_i ~=~ 
\left( \begin{array}{rrrr}
1 & 5 & 8 & 10 \\
  & 2 & 6 & 9 \\
  &   & 3 & 7 \\
  &   &   & 4 
\end{array} \right).
\label{eqn:hbw.eg}
\end{eqnarray}
That is, the order is firstly the diagonal, then the band above that,
followed by the second band above the diagonal etc.
Why is such a format adopted? 
For this example, if $\bW_i$ is diagonal then only the first 4 columns
of \texttt{wz} are needed. If $\bW_i$ is tridiagonal then only the
first~7 columns of \texttt{wz} are needed. 
If $\bW_i$ \textit{is} banded then \texttt{wz} needs not have
all $\frac12 M(M+1)$ columns; only~$M^*$ columns suffice, and the
rest of the elements of~$\bW_i$ are implicitly zero.
As well as reducing the size of \texttt{wz} itself in most cases, the
matrix-band format often makes the computation of \texttt{wz} very
simple and efficient. Furthermore, a Cholesky decomposition of a
banded matrix will be banded. A final reason is that sometimes we
want to input~$\bW_i$ into \VGAM: if \texttt{wz} is $M \times M \times
n$ then \texttt{vglm(\ldots, weights = wz)} will result in an error
whereas it will work if \texttt{wz} is an $n \times M^*$ matrix.



To facilitate the use of the matrix-band format,
a few auxiliary functions have been written.
In particular, there is \texttt{iam()} which gives the indices
for an array-to-matrix.
In the $4\times 4$ example above,
<<>>=
iam(NA, NA, M = 4, both = TRUE, diag = TRUE)
@
returns the indices for the respective array coordinates for
successive columns of matrix-band format
(see Equation~\ref{eqn:hbw.eg}).
If \texttt{diag = FALSE} then the first~4 elements in each vector
are omitted. Note that the first two arguments of 
\texttt{iam()} are not used here and have been assigned
\texttt{NA}s for simplicity.
For its use on the multinomial logit model, where
$(\bW_i)_{jj} = w_i\,\mu_{ij} (1-\mu_{ij}),\ j=1,\ldots,M$, and 
$(\bW_i)_{jk} = -w_i\,\mu_{ij} \mu_{ik},\ j\neq k$,
this can be programmed succinctly like
\begin{Code}
wz <- mu[, 1:M] * (1 - mu[, 1:M])
if (M > 1) {
    index <- iam(NA, NA, M = M, both = TRUE, diag = FALSE)
    wz <- cbind(wz, -mu[, index$row] * mu[, index$col])
}
wz <- w * wz
\end{Code}
(the actual code is slightly more complicated).
In general, \VGAM{}~family functions can be remarkably compact,
e.g.,
\texttt{acat()},
\texttt{cratio()}
and
\texttt{multinomial()} are all less than 120~lines of code each.










% ----------------------------------------------------------------------
\section{Extensions and utilities}
\label{sec:jsscat.extnUtil}

This section describes some useful utilities/extensions of the above.



\subsection{Marginal effects}
\label{sec:jsscat.extnUtil.margeff}


Models such as the multinomial logit and cumulative link models
model the posterior probability $p_{j} = \pr(Y=j|\bix)$ directly.
In some applications, knowing the derivative of~$p_{j}$
with respect to some of the~$x_k$ is useful;
in fact, often just knowing the sign is important.
The function \texttt{margeff()} computes the derivatives and
returns them as a $p \times (M+1) \times n$ array.
For the multinomial logit model it is easy to show
\begin{eqnarray}
\frac{\partial \, p_{j}(\bix_i)}{\partial \,
\bix_{i}}
&=&
p_{j}(\bix_i)
\left\{
 \bbeta_{j} -
\sum_{s=1}^{M+1}
p_{s}(\bix_i)
\,
 \bbeta_{s}
\right\},
\label{eqn:multinomial.marginalEffects}
\end{eqnarray}
while for
\texttt{cumulative(reverse = FALSE)}
we have
$p_{j} = \gamma_{j} - \gamma_{j-1} = h(\eta_{j}) - h(\eta_{j-1})$
where $h=g^{-1}$ is the inverse of the link function
(cf.~Table~\ref{tab:cat.quantities})
so that
\begin{eqnarray}
\frac{\partial \, p_{j}(\bix_{})}{\partial \,
\bix}
&=&
h'(\eta_{j}) \, \bbeta_{j} -
h'(\eta_{j-1}) \, \bbeta_{j-1} .
\label{eqn:cumulative.marginalEffects}
\end{eqnarray}




The function \texttt{margeff()} returns an array with these
derivatives and should handle any value of
\texttt{reverse} and \texttt{parallel}.








% ----------------------------------------------------------------------
\subsection[The xij argument]{The \texttt{xij} argument}
\label{sec:jsscat.extnUtil.xij}

There are many models, including those for categorical data,
where the value of an explanatory variable~$x_k$ differs depending
on which linear/additive predictor~$\eta_{j}$.
Here is a well-known example from {consumer choice} modeling.
Suppose an econometrician is interested in peoples'
choice of transport for travelling to work
and that there are four choices:
$Y=1$ for ``bus'',
$Y=2$ ``train'', 
$Y=3$ ``car'' and
$Y=4$ means ``walking''.
Assume that people only choose one means to go to work.
Suppose there are three covariates:
$X_2=$ cost,
$X_3=$ journey time, and
$X_4=$ distance.
Of the covariates only~$X_4$ (and the intercept~$X_1$)
is the same for all transport choices;
the cost and journey time differ according to the means chosen.
Suppose a random sample of~$n$ people is collected
from some population, and that each person has
access to all these transport modes.
For such data, a natural regression model would be a 
multinomial logit model with~$M=3$:
for $j=1,\ldots,M$, we have
$\eta_{j} =$
\begin{eqnarray}
\log \frac{\pr(Y=j)}{\pr(Y=M+1)}
&=&
\beta_{(j)1}^{*} +
\beta_{(1)2}^{*} \, (x_{i2j}-x_{i24}) +
\beta_{(1)3}^{*} \, (x_{i3j}-x_{i34}) +
\beta_{(1)4}^{*} \, x_{i4},
\label{eqn:xij.eg.gotowork}
\end{eqnarray}
where, for the~$i$th person,
$x_{i2j}$ is the cost for the~$j$th transport means, and
$x_{i3j}$ is the journey time of the~$j$th transport means.
The distance to get to work is $x_{i4}$; it has the same value
regardless of the transport means.


Equation~\ref{eqn:xij.eg.gotowork}
implies $\bH_1=\bI_3$ and $\bH_2=\bH_3=\bH_4=\bone_3$.
Note
also that if the last response category is used as the baseline or
reference group (the default of \texttt{multinomial()}) then $x_{ik,M+1}$
can be subtracted from $x_{ikj}$ for~$j=1,\ldots,M$---this
is the natural way $x_{ik,M+1}$ enters into the model.




Recall from~(\ref{gammod2}) that we had
\begin{equation}
\eta_j(\bix_i) ~=~ \bbeta_j^{\top} \bix_i ~=~
\sum_{k=1}^{p} \, x_{ik} \, \beta_{(j)k} .
\label{eqn:xij0}
\end{equation}
Importantly, this can be generalized to
\begin{equation}
\eta_j(\bix_{ij}) ~=~ \bbeta_j^{\top} \bix_{ij} ~=~
\sum_{k=1}^{p} \, x_{ikj} \, \beta_{(j)k} ,
\label{eqn:xij}
\end{equation}
or writing this another way (as a mixture or hybrid),
\begin{equation}
\eta_j(\bix_{i}^{*},\bix_{ij}^{*}) ~=~ 
\bbeta_{j}^{*T} \bix_{i}^{*} + \bbeta_{j}^{**T} \bix_{ij}^{*} .
\label{eqn:xij2}
\end{equation}
Often $\bbeta_{j}^{**} = \bbeta_{}^{**}$, say.
In~(\ref{eqn:xij2}) the variables in~$\bix_{i}^{*}$ are common to
all~$\eta_{j}$, and the variables in~$\bix_{ij}^{*}$ have
different values for differing~$\eta_{j}$.
This allows for covariate values that are specific to each~$\eta_j$,
a facility which is very important in many applications.


The use of the \texttt{xij} argument with the \VGAM{} family function
\texttt{multinomial()} has very important applications in economics.
In that field the term ``multinomial logit model'' includes a variety of
models such as the ``generalized logit model'' where (\ref{eqn:xij0})
holds, the ``conditional logit model'' where~(\ref{eqn:xij}) holds,
and the ``mixed logit model,'' which is a combination of the two,
where~(\ref{eqn:xij2}) holds.
The generalized logit model focusses on the individual as the unit of
analysis, and uses individual characteristics as explanatory variables,
e.g., age of the person in the transport example.
The conditional logit model assumes different values for each
alternative and the impact of a unit of~$x_k$ is assumed to be constant
across alternatives, e.g., journey time in the choice of transport mode.
Unfortunately, there is confusion in the literature for the terminology
of the models. Some authors call \texttt{multinomial()}
with~(\ref{eqn:xij0}) the ``generalized logit model''.
Others call the mixed
logit model the ``multinomial logit model'' and view the generalized
logit and conditional logit models as special cases.
In~\VGAM{} terminology there is no need to give different names to
all these slightly differing special cases. They are all still called
multinomial logit models, although it may be added that there are
some covariate-specific linear/additive predictors.
The important thing is that the framework accommodates~$\bix_{ij}$,
so one tries to avoid making life unnecessarily complicated.
And~\texttt{xij} can apply in theory to any VGLM and not just to the
multinomial logit model.
\cite{imai:king:lau:2008} present another perspective on the
$\bix_{ij}$ problem with illustrations from \pkg{Zelig}
\citep{Zelig:2009}.





\subsubsection[Using the xij argument]{Using the \texttt{xij} argument}
\label{sec:xij.sub}

\VGAM{} handles variables whose values depend on $\eta_{j}$,
(\ref{eqn:xij2}), using the \texttt{xij} argument.
It is assigned an~S formula or a list of \proglang{S}~formulas.
Each formula, which must have~$M$ \textit{different} terms,
forms a matrix that premultiplies a constraint matrix.
In detail, (\ref{eqn:xij0})~can be written in vector form as
\begin{equation}
\boldeta(\bix_i) ~=~ \bB^{\top} \bix_i ~=~
\sum_{k=1}^{p} \, \bH_{k} \, \bbeta_{k}^{*} \, x_{ik},
\label{eqn:xij0.vector}
\end{equation}
where
$\bbeta_{k}^{*} =
\left( \beta_{(1)k}^{*},\ldots,\beta_{(r_k)k}^{*} \right)^{\top}$
is to be estimated.
This may be written
\begin{eqnarray}
\boldeta(\bix_{i})
&=&
\sum_{k=1}^{p} \, \diag(x_{ik},\ldots,x_{ik}) \,
\bH_k \, \bbeta_{k}^{*}.
\label{eqn:xij.d.vector}
\end{eqnarray}
To handle~(\ref{eqn:xij})--(\ref{eqn:xij2})
we can generalize~(\ref{eqn:xij.d.vector}) to
\begin{eqnarray}
\boldeta_i
&=&
\sum_{k=1}^{p} \, \diag(x_{ik1},\ldots,x_{ikM}) \;
\bH_k \, \bbeta_{k}^{*}
\ \ \ \ \left(=
\sum_{k=1}^{p} \, \bX_{(ik)}^{*} \,
\bH_k \, \bbeta_{k}^{*} ,
\mathrm{\ say} \right).
\label{eqn:xij.vector}
\end{eqnarray}
Each component of the list \texttt{xij} is a formula having~$M$ terms
(ignoring the intercept) which
specifies the successive diagonal elements of the matrix~$\bX_{(ik)}^{*}$.
Thus each row of the constraint matrix may be multiplied by a different
vector of values.
The constraint matrices themselves are not affected by the
\texttt{xij} argument.





How can one fit such models in \VGAM{}?
Let us fit~(\ref{eqn:xij.eg.gotowork}).
Suppose the journey cost and time variables have had the
cost and time of walking subtracted from them.
Then,
using ``\texttt{.trn}'' to denote train,
\begin{Code}
fit2 <- vglm(cbind(bus, train, car, walk) ~ Cost + Time + Distance,
             fam = multinomial(parallel = TRUE ~ Cost + Time + Distance - 1),
             xij = list(Cost ~ Cost.bus + Cost.trn + Cost.car,
                        Time ~ Time.bus + Time.trn + Time.car),
             form2 = ~  Cost.bus + Cost.trn + Cost.car +
                        Time.bus + Time.trn + Time.car +
                        Cost + Time + Distance,
             data = gotowork)
\end{Code}
should do the job.
Here, the argument \texttt{form2} is assigned a second \proglang{S}~formula which
is used in some special circumstances or by certain types
of~\VGAM{} family functions.
The model has $\bH_{1} = \bI_{3}$ and $\bH_{2} = \bH_{3} = \bH_{4} = \bone_{3}$
because the lack of parallelism only applies to the intercept.
However, unless \texttt{Cost} is the same as \texttt{Cost.bus} and
\texttt{Time} is the same as \texttt{Time.bus},
this model should not be plotted with \texttt{plotvgam()};
see the author's homepage for further documentation.


By the way,
suppose 
$\beta_{(1)4}^{*}$
in~(\ref{eqn:xij.eg.gotowork})
is replaced by~$\beta_{(j)4}^{*}$.
Then the above code but with
\begin{Code}
  fam = multinomial(parallel = FALSE ~ 1 + Distance),
\end{Code}
should fit this model.
Equivalently,
\begin{Code}
  fam = multinomial(parallel = TRUE ~ Cost + Time - 1),
\end{Code}






\subsubsection{A more complicated example}
\label{sec:xij.complicated}

The above example is straightforward because the
variables were entered linearly. However, things
become more tricky if data-dependent functions are used in
any \texttt{xij} terms, e.g., \texttt{bs()}, \texttt{ns()} or \texttt{poly()}.
In particular, regression splines such as \texttt{bs()} and \texttt{ns()}
can be used to estimate a general smooth function~$f(x_{ij})$, which is
very useful for exploratory data analysis.



Suppose we wish to fit the variable \texttt{Cost} with a smoother.
This is possible with regression splines and using a trick.
Firstly note that
\begin{Code}
fit3 <- vglm(cbind(bus, train, car, walk) ~ ns(Cost) + Time + Distance,
             multinomial(parallel = TRUE ~ ns(Cost) + Time + Distance - 1),
             xij = list(ns(Cost) ~ ns(Cost.bus) + ns(Cost.trn) + ns(Cost.car),
                        Time ~ Time.bus + Time.trn + Time.car),
             form2 = ~  ns(Cost.bus) + ns(Cost.trn) + ns(Cost.car) +
                        Time.bus + Time.trn + Time.car +
                        ns(Cost) + Cost + Time + Distance,
             data = gotowork)
\end{Code}
will \textit{not} work because the basis functions for
\texttt{ns(Cost.bus)}, \texttt{ns(Cost.trn)} and \texttt{ns(Cost.car)}
are not identical since the knots differ.
Consequently, they represent different functions despite
having common regression coefficients.


Fortunately, it is possible to force the~\texttt{ns()} terms
to have identical basis functions by using a trick:
combine the vectors temporarily.
To do this, one can let
\begin{Code}
NS <- function(x, ..., df = 3)
      ns(c(x, ...), df = df)[1:length(x), , drop = FALSE]
\end{Code}
This computes a natural cubic B-spline evaluated at~\texttt{x} but it uses the
other arguments as well to form an overall vector from which to obtain
the (common) knots.
Then the usage of \texttt{NS()} can be something like
\begin{Code}
fit4 <- vglm(cbind(bus, train, car, walk) ~ NS(Cost.bus, Cost.trn, Cost.car)
                                          + Time + Distance,
             multinomial(parallel = TRUE ~  NS(Cost.bus, Cost.trn, Cost.car)
                                          + Time + Distance - 1),
             xij = list(NS(Cost.bus, Cost.trn, Cost.car) ~
                        NS(Cost.bus, Cost.trn, Cost.car) +
                        NS(Cost.trn, Cost.car, Cost.bus) +
                        NS(Cost.car, Cost.bus, Cost.trn),
                        Time ~ Time.bus + Time.trn + Time.car),
             form2 = ~  NS(Cost.bus, Cost.trn, Cost.car) +
                        NS(Cost.trn, Cost.car, Cost.bus) +
                        NS(Cost.car, Cost.bus, Cost.trn) +
                        Time.bus + Time.trn + Time.car +
                        Cost.bus + Cost.trn + Cost.car +
                        Time + Distance,
             data = gotowork)
\end{Code}
So \texttt{NS(Cost.bus, Cost.trn, Cost.car)}
is the smooth term for
\texttt{Cost.bus}, etc.
Furthermore, \texttt{plotvgam()} may be applied to
\texttt{fit4}, in which case the fitted regression spline is plotted
against its first inner argument, viz.~\texttt{Cost.bus}.


One of the reasons why it will predict correctly, too,
is due to ``smart prediction''
\citep{Rnews:Yee:2008}.



\subsubsection{Implementation details} 
\label{sec:jss.xij.implementationDetails} 

The~\texttt{xij} argument operates \textit{after} the
ordinary $\bX_{\sVLM}$ matrix is created. Then selected columns
of~$\bX_{\sVLM}$ are modified from the constraint matrices, \texttt{xij}
and~\texttt{form2} arguments. That is, from \texttt{form2}'s model
matrix $\bX_{\sformtwo}$, and the~$\bH_k$. This whole operation
is possible because $\bX_{\sVLM}$ remains structurally the same.
The crucial equation is~(\ref{eqn:xij.vector}).


Other \texttt{xij} examples are given in the online help of
\texttt{fill()} and \texttt{vglm.control()},
as well as at the package's webpage.











% ----------------------------------------------------------------------
\section{Discussion}
\label{sec:jsscat.discussion}


This article has sought to convey how VGLMs/VGAMs are well suited for
fitting regression models for categorical data. Its primary strength
is its simple and unified framework, and when reflected in software,
makes practical CDA more understandable and efficient. Furthermore,
there are natural extensions such as a reduced-rank variant and
covariate-specific~$\eta_{j}$. The \VGAM{}~package potentially offers
a wide selection of models and utilities.


There is much future work to do.
Some useful additions to the package include:
\begin{enumerate}

\item
Bias-reduction \citep{firt:1993} is a method for removing the~$O(n^{-1})$
bias from a maximum likelihood estimate. For a substantial class of
models including GLMs it can be formulated in terms of a minor adjustment
of the score vector within an IRLS algorithm \citep{kosm:firt:2009}.
One by-product, for logistic regression, is that while the maximum
likelihood estimate (MLE) can be infinite, the adjustment leads to
estimates that are always finite. At present the \R{}~package \pkg{brglm}
\citep{Kosmidis:2008} implements bias-reduction for a number of models.
Bias-reduction might be implemented by adding an argument
\texttt{bred = FALSE}, say, to some existing \VGAM{} family functions.


\item
Nested logit models were developed to overcome a fundamental shortcoming
related to the multinomial logit model, viz.~the independence of
irrelevant alternatives~(IIA) assumption. Roughly, the multinomial logit
model assumes the ratio of the choice probabilities of two alternatives
is not dependent on the presence or absence of other alternatives in
the model. This presents problems that are often illustrated by the
famed red bus-blue bus problem.




\item
The generalized estimating equations (GEE) methodology is largely
amenable to IRLS and this should be added to the package in the future
\citep{wild:yee:1996}.


\item
For logistic regression \proglang{SAS}'s \code{proc logistic} gives
a warning if the data is {completely separate} or {quasi-completely
separate}. Its effects are that some regression coefficients tend to~$\pm
\infty$. With such data, all (to my knowledge) \R{}~implementations
give warnings that are vague, if any at all, and this is rather
unacceptable \citep{alli:2004}. The \pkg{safeBinaryRegression} package
\citep{Konis:2009} overloads \code{glm()} so that a check for the
existence of the MLE is made before fitting a binary response GLM.


\end{enumerate}


In closing, the \pkg{VGAM} package is continually being developed,
therefore some future changes in the implementation details and usage
may occur. These may include non-backward-compatible changes (see the
\code{NEWS} file.) Further documentation and updates are available at
the author's homepage whose URL is given in the \code{DESCRIPTION} file.



% ----------------------------------------------------------------------
\section*{Acknowledgments}

The author thanks Micah Altman, David Firth and Bill Venables for helpful
conversations, and Ioannis Kosmidis for a reprint.
Thanks also to The Institute for Quantitative Social Science at Harvard
University for their hospitality while this document was written during a
sabbatical visit.





\bibliography{categoricalVGAMbib}

\end{document}




